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

HAL
open science

# Starlet higher order statistics for galaxy clustering and weak lensing

Virginia Ajani, Joachim Harnois-Déraps, Valeria Pettorino, Jean-Luc Starck

► To cite this version:

Virginia Ajani, Joachim Harnois-Déraps, Valeria Pettorino, Jean-Luc Starck. Starlet higher order statistics for galaxy clustering and weak lensing. Astronomy and Astrophysics - A&A, 2023, 672, pp.L10. 10.1051/0004-6361/202245510. hal-03892296

HAL Id: hal-03892296

https://hal.science/hal-03892296

Submitted on 13 Oct 2023

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
LETTER TO THE EDITOR

# Starlet higher order statistics for galaxy clustering and weak lensing

Virginia Ajani<fnref n="1" />,<fnref n="2" />, Joachim Harnois-Déraps<fnref n="3" />, Valeria Pettorino<fnref n="1" />, and Jean-Luc Starck<fnref n="1" />

Received 18 November 2022 / Accepted 15 March 2023

## ABSTRACT

We present a first application to photometric galaxy clustering and weak lensing of wavelet-based multi-scale (beyond two points) summary statistics: starlet peak counts and starlet $\ell_1$-norm. Peak counts are the local maxima in the map, and $\ell_1$-norm is computed via the sum of the absolute values of the starlet (wavelet) decomposition coefficients of a map, providing a fast multi-scale calculation of the pixel distribution, encoding the information of all pixels in the map. We employ the cosmo-SLICS simulations sources and lens catalogues, and we compute wavelet-based non-Gaussian statistics in the context of combined probes and their potential when applied to the weak-lensing convergence maps and galaxy maps. We obtain forecasts on the matter density parameter $\Omega_m$, the reduced Hubble constant $h$, the matter fluctuation amplitude $\sigma_8$, and the dark energy equation of state parameter $w_0$. In our setting for this first application, we consider the two probes to be independent. We find that the starlet peaks and the $\ell_1$-norm represent interesting summary statistics that can improve the constraints with respect to the power spectrum, even in the case of photometric galaxy clustering and when the two probes are combined.

**Key words.** large-scale structure of Universe – methods: statistical – cosmological parameters

## 1. Introduction

Past, present, and future cosmological surveys, such as the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS) (Heymans et al. 2012), the Kilo-Degree Survey (KiDS) (Heymans et al. 2021), the Dark Energy Survey (DES) (Abbott et al. 2019; To et al. 2021), Hyper SuprimeCam (HSC) Mandelbaum & Hyper Suprime-Cam (HSC) Collaboration 2017, Euclid (Laureijs et al. 2011), and the Rubin Observatory (Abell 2009), use weak gravitational lensing by large-scale structure and galaxy clustering as the main physical probes to explore the unknown properties of dark energy and dark matter. The correlation between the positions of foreground galaxies with the shapes of the source galaxies, also known as 3 × 2 pt analysis, is very sensitive to the amplitude of matter clustering in the low-redshift universe and represents a powerful test of consistency between the growth of structure and the expansion history of the Universe (Porredon et al. 2021; Abbott et al. 2022). Since they bring complementary cosmological information, the combination of cosmic shear, galaxy clustering, and galaxy-galaxy lensing has proved to be very powerful to test cosmological models and modified gravity models and in enhancing the constraining power (Zhang et al. 2007; Eriksen & Gaztañaga 2018; Heymans et al. 2021; Euclid Collaboration 2020; Tutusaus et al. 2020 and references therein) and for dealing with systematic effects (Mandelbaum et al. 2013; Camera et al. 2016; Harnois-Déraps et al. 2018; Joachimi et al. 2021; Krause et al. 2021) as they are typically distinct and uncorrelated in dif-

ferent probes. Furthermore, the increasing precision reached with next-generation surveys will enable us to access the non-Gaussian part of cosmological signals, induced by the non-linear evolution of structure on small scales and low redshifts, which is not captured with second-order summary statistics alone. Specifically for weak lensing, a rich literature proposing several non-Gaussian statistics (Euclid Collaboration 2023), such as Minkowski functionals (e.g., Kratochvil et al. 2012 and Parroni et al. 2020), higher order moments (e.g., Petri et al. 2016 and Gatti et al. 2020), bispectrum (Takada & Jain 2004; Coulton et al. 2019), peak counts (Kruse & Schneider 1999; Dietrich & Hartlap 2010; Liu et al. 2015; Lin & Kilbinger 2015; Peel et al. 2017; Martinet et al. 2017; Li et al. 2019; Ajani et al. 2020; Zücher et al. 2022b; Ayçoberry et al. 2023), Betti numbers (Parroni et al. 2021), the scattering transform (Cheng et al. 2020), wavelet phase harmonic statistics (Allys et al. 2020), and machine learning-based methods (e.g., Fluri et al. 2018 and Shirasaki et al. 2021), is catching the attention of the community. The $\ell_1$-norm of wavelet coefficients of weak-lensing convergence maps has been proposed (Ajani et al. 2021) as a new summary statistics for weak lensing as it provides a unified framework to perform a multi-scale analysis that takes into account the information encoded in all pixels of the map.

Moreover, Grewal et al. (2022) have explored the addition of DES-like and LSST-like mock clustering maps to an analysis of Minkowski functionals; they find a significant improvement with respect to a weak lensing-only analysis. Kacprzak & Fluri (2022) show that a deep learning analysis of combined weak

<footnote n="1">¹ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">² Institute for Particle Physics and Astrophysics, Department of Physics, ETH Zürich, Wolfgang Pauli Strasse 27, 8093 Zürich, Switzerland
e-mail: vajani@phys.ethz.ch</footnote>

<footnote n="3">³ School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, NE1 7RU Newcastle-upon-Tyne, UK</footnote>
lensing and galaxy clustering at the map level can help break the degeneracy between cosmological and astrophysical parameters. As the wavelet-based statistics presented in Ajani et al. (2020, 2021), and Zücher et al. (2022a) were limited to weak lensing, we propose in this work a first application of such statistics in the context of joint weak lensing and galaxy clustering with the final goal of showing the potential of wavelet-based non-Gaussian statistics in the framework of probe combination. The scope of this study is to extend the application of starlet peaks and $ℓ_1$-norm to photometric galaxy clustering and its combination with weak lensing. Specifically, we employ lensing and clustering data from the cosmo-SLICS simulations to build convergence and density maps from which we obtain various weak lensing and clustering predictions. We perform a likelihood analysis using the power spectrum, the starlet peaks, and $ℓ_1$-norm, and constrain the matter density parameter $\Omega_m$, the reduced Hubble constant $h$, the matter fluctuation amplitude $\sigma_8$, and the dark energy equation of state parameter $w_0$.

The Letter is structured as follows. In Sect. 2 we provide some background definitions for weak lensing and photometric galaxy clustering. In Sect. 3 we present the mock data, the catalogues, and our map-making procedure. We describe our analysis in Sect. 4. In Sect. 5 we present and discuss the results, limitations, and perspectives of this work. We conclude in Sect. 6.

## 2. Background

### 2.1. Weak lensing

The distortions caused by gravitational lensing in the original source image can be summarised in the distortion matrix:

$$ \mathcal{A} = \begin{pmatrix} 1 - \kappa - \gamma_1 & -\gamma_2 \\ -\gamma_2 & 1 - \kappa + \gamma_1 \end{pmatrix} = (1 - \kappa) \begin{pmatrix} 1 - g_1 & -g_2 \\ -g_2 & 1 + g_1 \end{pmatrix}. $$

Here $\gamma = (\gamma_1, \gamma_2)$ is the shear, also defined in terms of the lensing potential $\psi$ as $\gamma_1 = \frac{1}{2}(\partial_1\partial_1 - \partial_2\partial_2)\psi$ and $\gamma_2 = \partial_1\partial_2\psi$, which describes the anisotropic distortion of the image; $\kappa$ is the convergence field, also defined as $\kappa = \frac{1}{2}\nabla^2\psi$, which describes the isotropic dilation of the source; and $g$ is the reduced shear. The observed ellipticity $\epsilon^{\text{obs}}$ is related with the intrinsic ellipticity of the galaxy $\epsilon^s$ and the reduced shear through the relation (Seitz & Schneider 1997)

$$ \epsilon^{\text{obs}} = \frac{\epsilon^s + g}{1 + g^* \epsilon^s}, \qquad (1) $$

where $g^*$ is the complex conjugate of $g$. In the weak-lensing regime, where $|\gamma|, \kappa \ll 1$, this relation can be approximated as $\epsilon^{\text{obs}} \approx \epsilon^s + \gamma$. Moreover, if the intrinsic ellipticity of galaxies has no preferred orientation, its expectation value vanishes, $\langle \epsilon^s \rangle = 0$, and the observed ellipticity can be employed as an unbiased estimator of the reduced shear, as $\langle \epsilon^{\text{obs}} \rangle \approx \gamma$ (Kilbinger 2015). In this study, we work under these assumptions; including the effect of the intrinsic alignments goes beyond the purposes of this work. The shear field can be related to the underlying matter density field using the relations connecting $\gamma$ and $\kappa$ to the lensing potential along with the fact that $\kappa$ is a line-of-sight integral of the matter density field. Working in Fourier space, we obtain an estimator of the convergence field starting from the shear field through the Kaiser-Squires inversion (KS) (Kaiser & Squires 1993) as

$$ \hat{k} = \frac{k_1^2 - k_2^2}{k^2} \hat{\gamma}_1 + \frac{2k_1 k_2}{k^2} \hat{\gamma}_2, \qquad (2) $$

where $k^2 = k_1^2 + k_2^2$.

### 2.2. Photometric galaxy clustering

As second-order statistics for photometric galaxy clustering, we employed the angular power spectra in harmonic space, $C^{\text{GG}}(\ell)$, measured on the galaxy maps. Given the redshift distribution of the galaxies $n_i^G(z)$ in each tomographic bin $i$ of a photometric survey, the average galaxy density in this redshift bin is given by (Euclid Collaboration 2020)

$$ \bar{n}_i^G = \int_{z_{\text{min}}}^{z_{\text{max}}} dz n_i^G(z). \qquad (3) $$

Then the radial weight function for a given tomographic bin $i$ for galaxy clustering can be defined as

$$ W_i^G(z) = \frac{n_i^G(z) H(z)}{\bar{n}_i^G c}. \qquad (4) $$

The observable spectrum is theoretically given by the expression

$$ C_{ij}^{\text{GG}}(\ell) = \int dz \frac{W_i^G(z) W_j^G(z)}{H(z) r^2(z)} P_{\text{gg}}^{\text{photo}} \left[ \frac{\ell + 1/2}{r(z)}, z \right], \qquad (5) $$

where $P_{\text{gg}}^{\text{photo}}(\ell, z)$ is the galaxy–galaxy power spectrum, which, under the assumption of a linear galaxy bias, is linked to the matter power spectrum $P_{\delta\delta}(\ell, z)$ through the bias $b_g$:

$$ P_{\text{gg}}^{\text{photo}}(\ell, z) = [b_{\text{g}}^{\text{photo}}(z)]^2 P_{\delta\delta}(\ell, z). \qquad (6) $$

In this first application we consider only the auto-spectra ($i = j$) for our power spectrum measurements, and that the galaxy bias $b_g$ is fixed in the simulations.

## 3. Mock data

We modelled our non-Gaussian statistics with mock data constructed from the cosmo-SLICS (Harnois-Déraps et al. 2019), a suite of high-resolution N-body simulations consisting of 25 different cosmologies organised in a Latin hypercube obtained by varying the matter density parameter $\Omega_m$ in the range [0.10, 0.55], the reduced Hubble constant $h$ in the range [0.60, 0.82], the combination between $\Omega_m$ and the matter fluctuation amplitude $\sigma_8$, $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$ in the range [0.60, 0.90], and the dark energy equation of state parameter $w_0$ between [-2.0, -0.5]. Each gravity-only calculation evolves 1536³ particles in a box of comoving side of 505 $h^{-1}$ Mpc with the Poisson solver CUBEP³M (Harnois-Déraps et al. 2013). These were then ray-traced under the Born approximation to construct the sources and lenses catalogues of 100 deg². To reduce the impact of cosmic variance, two simulations were provided for each cosmology with different initial conditions, each re-sampled to create 25 pseudo-independent light cones.

### 3.1. Source catalogue

For the weak lensing analysis we used the KiDS-1000-like sources catalogue described in Harnois-Déraps et al. (2021), which mimics the survey properties of the KiDS-1000 data described in Giblin et al. (2018) and Hildebrandt et al. (2021). These catalogues are available for five different tomographic redshift bins, with a corresponding galaxy number density of $n_{\text{gal}}$ = [0.62, 1.18, 1.85, 1.26, 1.31]. For each source galaxy, the catalogue contains the positions $(x, y)$ in arcmin and in redshift $z$, the values for the components of the true shear ($\gamma_1, \gamma_2$) and the values for the components of the observed ellipticities ($\epsilon_1^{\text{obs}}, \epsilon_2^{\text{obs}}$).
## 3.2. Building convergence maps

To get the convergence map, we used the KS inversion: we projected the shear onto a 600² Cartesian grid (resolution of 1 arcmin per pixel), which we then smoothed with a Gaussian filter of width 0.7 arcmin following Giblin et al. (2018) to reduce the impact of potential empty pixels. This was implemented using the Python package lenspack<fnref n="1" />. We provide an example of our validation against the theoretical prediction in Appendix C.

## 3.3. Lens catalogues

To perform the galaxy clustering analysis, we employed the KiDS-1000-like lens catalogues of luminous red galaxies (LRGs) provided by the cosmo-SLICS simulations. This simulated foreground LRG sample reproduces the characteristics described in Vakili et al. (2020) and is split into four tomographic bins (see Burger et al. 2022, for a recent study employing the same catalogues). Each catalogue contains the position and redshift of each galaxy. The LRGs trace the underlying projected mass sheets, hence their redshift distribution is heavily discretised in bins of $\Delta\chi = 257.5$ Mpc h⁻<fnref n="1" />. The catalogues are generated assuming linear galaxy bias, whose values and uncertainties are from Vakili et al. (2020) and summarised in Table B.1.

## 3.4. Building galaxy maps

Starting from the lens catalogue described in Sect. 3.3, we built galaxy maps for each light cone by assigning the galaxies onto a 600² grid. We validated our maps by comparing our measured power spectra against the theoretical prediction obtained with the Python library pyccl<fnref n="2" />. The shot noise of the galaxy maps was approximated as Gaussian noise with a variance equal to the mean of each galaxy map.

As the map is intrinsically noisy (due to the presence of shot noise), while the theoretical prediction is obtained for a noiseless case, we computed the power spectrum of the noise map and then subtracted it from the power spectrum of the map. Since we did not have a model including non-linear bias, we started with a conservative analysis and included scales up to $\ell_{max} = 780$. An example of the comparison between our measurement on the simulated galaxy map and the corresponding theoretical prediction can be found in Fig. C.1.

# 4. Analysis

## 4.1. Summary statistics

We compute our summary statistics on convergence maps for weak lensing and galaxy maps for galaxy clustering. When computing the power spectra, we smooth the noisy convergence and galaxy maps with a Gaussian kernel of size $\theta_{ker} = 1$ arcmin. For weak lensing we use 33 bins in the range $\ell = [180, 2040]$, which lies within the region where we are confident that the power spectra are in agreement with the theoretical predictions. For galaxy clustering, we use 14 bins in the range $\ell = [120, 780]$, following the pessimistic setting of Euclid Collaboration (2020). These regions in $\ell$ are shown in Appendix C. For the non-Gaussian statistics, we filter the noisy maps with a starlet filter (Starck et al. 2007), a wavelet transform that allows for a fast multi-scale image analysis (see Appendix A). It has proved to be

powerful for extracting cosmological information in the context of non-Gaussian statistics (e.g., Lin & Kilbinger 2015; Lin et al. 2016; Peel et al. 2017; Ajani et al. 2020, 2021). A visualisation of this decomposition is provided for both probes in Fig. 1. Following Ajani et al. (2020), we then computed peak counts on the starlet decomposed map, finding pixels with values larger than their eight neighbours. They can be quantified as a function of height, namely the pixel’s value, or the signal-to-noise ratio (S/N) (Fan et al. 2010; Maturi et al. 2010; Shan et al. 2017). We compute the S/N as

$$S/N \equiv \frac{(W * m)(\theta_{ker})}{\sigma_{n}^{\text{filt}}(\theta_{ker})}, \qquad (7)$$

where $W(\theta_{ker})$ is the Starlet filter, $m$ is the noisy map, and $*$ indicates the convolution operation. The quantity $\sigma_{n}^{\text{filt}}$ is the standard deviation of the smoothed noise maps: for each weak-lensing map we estimate the shape noise as the standard deviation of the observed ellipticities from the source catalogue and for the galaxy clustering maps as the mean of each galaxy density maps. We estimate the noise level at each filter scale with automatic noise estimation from the multi-resolution support (Starck & Murtagh 1998). On the same maps, we also compute the starlet $\ell_1$-norm, which is the absolute value of the sum of all wavelet coefficients for each scale of the decomposition within a chosen range of amplitude (Ajani et al. 2021). For both starlet peaks and $\ell_1$-norm we use for inference four scales corresponding to [2', 4', 8', 16', coarse], with 29 linearly spaced bins for each scale between the minimum and maximum values of each S/N map.

## 4.2. Likelihood

To perform Bayesian inference and obtain the probability distributions of the cosmological parameters we want to constrain, we use a Gaussian likelihood for a cosmology-independent covariance:

$$\log \mathcal{L}(\theta) = \frac{1}{2} \sum_{\alpha} (d - \mu(\theta))_{\alpha}^{T} C_{\alpha}^{-1} (d - \mu(\theta))_{\alpha}. \qquad (8)$$

Here $d$ is the data array given by the average of the summary statistics over the different realisations at fiducial model, $C$ is the covariance matrix of the observable defined in the next section, $\mu$ is the model prediction as a function of the cosmological parameters $\theta$, and $\alpha$ denotes the physical probe (i.e., weak lensing and galaxy clustering). At second order, the cross-correlation between the galaxy clustering and weak-lensing mock catalogues from two different surveys can be neglected if the overlap between the two surveys is sufficiently small (Heymans et al. 2021). We are aware that in our case the galaxy density and convergence maps can in reality be correlated, as shown in Fig. 1, where the structures of the galaxy map are traced by the lensing convergence map. For this first application, we ignore these correlations that could actually bring additional information. For map level studies it has been shown that cross-probes can help in breaking degeneracies and in increasing the constraining power (Kacprzak & Fluri 2022; Grewal et al. 2022). We leave for a future study the inclusion of their cross-correlation. The cosmological parameters are the ones that vary in the simulations, namely $\theta = [\Omega_m, \sigma_8, h, w_0]$. To model the summary statistics at arbitrary cosmologies, we employ an interpolation with Gaussian processes regression (Rasmussen & Williams 2005), using the scikit-learn<fnref n="3" /> Python package. For more details, the emulator

<footnote n="1">¹ https://github.com/CosmoStat/lenspack</footnote>

<footnote n="2">² https://pypi.org/project/pyccl/</footnote>

<footnote n="3">³ https://scikit-learn.org/</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Starlet decomposition for a convergence map (top panel) and a galaxy map (bottom panel). The value in arcmin corresponds to the resolution of the corresponding starlet scale. The resolution of the original map is 1 arcmin.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Forecast contours at confidence levels of 68% and 95% using the weak-lensing power spectrum (light blue) compared to galaxy clustering power spectrum (orange) and their combination as independent probes (violet).</figcaption></figure>

used in this Letter is the same as the one employed in Li et al. (2019) and in Ayçoberry et al. (2023)<fnref n="4" />.

### 4.3. Covariance matrix

To compute the covariance matrix we employ the SLICS simulations (Harnois-Déraps et al. 2018). These are obtained as described in Sect. 3 but are specifically designed for the estimation of covariance matrices: they consist of 800 fully independent ΛCDM runs in which the cosmological parameters are fixed to the cosmology [Ω_m, σ_8, h, n_s, Ω_b] = [0.2905, 0.826, 0.6898, 0.969, 0.0473] and the random seeds in the initial conditions are varied. The covariance matrix elements are computed as

$$C_{ij} = \sum_{r=1}^{N} \frac{(x_i^r - \mu_i)(x_j^r - \mu_j)}{N-1}, \quad (9)$$

where N = 800 is the number of observations, $x_i^r$ is the value of the ith data element for a given realisation r, and $\mu_i$ is the mean over all realisations. We take into account the loss of information due to the finite number of bins and realisations by adopting the estimator introduced by Hartlap et al. (2007) for the inverse of the covariance matrix $C^{-1} = \frac{N-n_{bins}-2}{N-1}C_*^{-1}$, where $n_{bins}$ is the number of bins and $C_*$ the covariance matrix computed for the power spectrum, the peak counts, and the $\ell_1$-norm, whose elements are given by Eq. (9). We area-rescale the covariance matrix to model the statistical accuracy of the KiDS-1000 survey, with $A_{KiDS-1000} = 1000 \text{ deg}^2$.

### 4.4. MCMC simulations and posterior distributions

We explore and constrain the parameter space with the emcee package, the affine invariant ensemble sampler for Markov chain Monte Carlo (MCMC) introduced by Foreman-Mackey et al. (2013). We assume a flat prior for all parameters over the range modelled by the cosmo-SLICS, with the Gaussian likelihood defined in Eq. (8). We use 120 walkers initialised in a tiny Gaussian ball of radius $10^{-3}$ around the fiducial cosmology [$\Omega_m, \sigma_8, h, w_0$] = [0.2905, 0.8364, 0.6898, -1].

## 5. Results

We first derive constraints with the power spectrum for both weak lensing and galaxy clustering, shown in Fig. 2. We see how, in our setting, the galaxy clustering alone outperforms the weak-lensing contours, breaking the degeneracies especially in the ($\Omega_m, \sigma_8$) plane. At the same time, combining them brings additional information, improving the constraining power and the degeneracy break even more, underlying the complementarity of cosmological information encoded in the two probes.

In Fig. 3 we show the comparison among the constraints obtained using the galaxy clustering power spectrum, the starlet peak counts measured on galaxy density maps, and the galaxy clustering $\ell_1$-norm to explore this statistics in the context of galaxy clustering. We see how, as for weak lensing, the starlet peaks improve the constraining power with respect to the power spectrum and the $\ell_1$-norm outperforms the other two. We are aware that the massive improvement in constraining power for galaxy clustering is very likely overestimated due to the fact that the galaxy bias in these simulations is kept fixed and assumed to be perfectly linear. In Appendix B we perform an analysis using the power spectrum with four different values of the galaxy bias parameters sampled from the measured uncertainty on these, as

<footnote n="4">⁴ The code for the GP emulator is publicly available at https://github.com/CosmoStat/shear-pipe-peaks</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Forecast contours at confidence levels of 68% and 95% using the power spectrum (orange), starlet peaks (light red), and the starlet $l_1$-norm(green) computed on galaxy density maps.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Forecast contours at confidence levels of 68% and 95% using the starlet $l_1$-norm measured on convergence maps (blue), on galaxy density maps (green), and their combination as independent probes (red).</figcaption></figure>

reported in Table B.1. The impact is major, where a $1\sigma$ shift in galaxy bias offsets the inferred $\sigma_8$ by $2\sigma$. To assess then the constraining power of the starlet peaks and the $l_1$-norm with respect to the power spectrum within our setting in a fairer comparison, we show the results for galaxy clustering alone, where the setting for all statistics in terms of maps employed and assumptions on the galaxy bias is the same. As in this setting the bias is fixed for all three statistics, we are more confident in claiming that starlet peaks and the starlet $l_1$-norm both significantly improve the constraints, also in the case of galaxy clustering, with the $l_1$-norm outperforming the peaks. Given this result, to fully exploit the potential of multi-scale non-Gaussian statistics in the context of map-level probe combination, we get the constraints using the combination of starlet $l_1$-norm for galaxy clustering and weak lensing. We show these results in Fig. 4 for the weak-lensing starlet $l_1$-norm, the galaxy clustering starlet $l_1$-norm, and their combination. We see how, also in this context, the combination of the two probes might help in tightening the constraints, suggesting that this beyond-two-point statistics computed at the map level could be potentially very powerful for probe combination. We are aware that the relative importance of clustering will degrade once the bias is allowed to vary. From the results of Appendix B we see how it is crucial to take this into account to avoid shifts on the inference contours for a real data application. We leave to a future study a deeper investigation of how we might include this uncertainty, for example by emulating the bias values across the uncertainties reported in Table B.1 and treating it as nuisance parameter.

## 6. Conclusions

In this Letter we proposed using for the first time starlet-based non-Gaussian statistics on weak-lensing convergence maps and on galaxy maps to constrain cosmological parameters. We employed a mock source catalogue for weak lensing and a mock lens catalogue for galaxy clustering available from the cosmological

simulations provided with a galaxy bias estimated from the LRG set in the fourth data release of the Kilo-Degree Survey. We performed a single probe and a combined non-Gaussian statistics analysis at the map level. We find that using the starlet-based statistics, such as starlet peak counts and the starlet $l_1$-norm also in the context of galaxy clustering, can potentially help in improving the constraints, and that this gain can be also obtained when combining weak lensing and photometric galaxy clustering. We are aware that there are many important aspects to explore before a real data application, such as the validity of linear bias approximation in the context of such higher order statistics. In Appendix B we show the results we obtained considering a different level of uncertainty on the bias and its impact on the constraints for the power spectrum. The next step was to build the framework to include galaxy–galaxy lensing analysis as well as cross-correlations between the two probes in the covariance matrix. We expect these to help in breaking degeneracies and improving the constraining power. However, to exploit this step for a real data application, several aspects need to be considered, such as observing conditions, masks, baryonic effects, and source clustering, among others. The goal of this first application represents a proof of concept to explore the applicability of these multi-scale map-based statistics already proven powerful for weak lensing also to galaxy clustering as a starting point for an analysis that accounts for the inclusion of systematic effects for both probes.

**Acknowledgements.** VA wishes to thank François Lanusse, for useful discussions and help during the first phase of development of the project. VA wishes to thank Tomasz Kacprzak and Alexandre Refregier for useful discussions. VA acknowledges support by the Centre National d'Etudes Spatiales and the project Initiative d'Excellence (IdEx) of Université de Paris (ANR-18-IDEX-0001) at the start and during the preparation of this work. VA acknowledges support by the grant Cosmological Weak Lensing and Neutral Hydrogen, 200021_192243 of Prof. Alexandre Refregier. JHD acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). The N-body simulations were enabled by the Digital Research Alliance of Canada (https://alliancecan.ca/).
References

Abbott, T. M. C., Abdalla, F. B., Avila, S., et al. 2019, Phys. Rev. D, 99, 123505

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520

Ajani, V., Peel, A., Pettorino, V., et al. 2020, Phys. Rev. D, 102, 103531

Ajani, V., Starck, J.-L., & Pettorino, V. 2021, A&A, 645, L11

Allys, E., Marchand, T., Cardoso, J.-F., et al. 2020, Phys. Rev. D, 102, 103506

Ayçoberry, E., Ajani, V., Guinot, A., et al. 2023, A&A, 671, A17

Burger, P., Friedrich, O., Harnois-Déraps, J., & Schneider, P. 2022, A&A, 661, A137

Camera, S., Harrison, I., Bonaldi, A., & Brown, M. L. 2016, MNRAS, 464, 4747

Cheng, S., Ting, Y.-S., Ménard, B., & Bruna, J. 2020, MNRAS, 499, 5902

Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2

Coulton, W. R., Liu, J., Madhavacheril, M. S., Böhm, V., & Spergel, D. N. 2019, J. Cosmol. Astropart. Phys., 2019, 043

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

Eriksen, M., & Gaztañaga, E. 2018, MNRAS, 480, 5226

Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191

Euclid Collaboration (Ajani, V., et al.) 2023, A&A, submitted [arXiv:2301.12890]

Fan, Z., Shan, H., & Liu, J. 2010, ApJ, 719, 1408

Fluri, J., Kacprzak, T., Refregier, A., et al. 2018, Phys. Rev. D, 98, 123518

Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306

Gatti, M., Chang, C., Friedrich, O., et al. 2020, MNRAS, 498, 4060

Giblin, B., Heymans, C., Harnois-Déraps, J., et al. 2018, MNRAS, 480, 5529

Grewal, N., Zuntz, J., Tröster, T., Amon, A., et al. 2022, Minkowski Functionals in Joint Galaxy Clustering & Weak Lensing Analyses

Harnois-Déraps, J., Pen, U.-L., Iliev, I. T., Merz, H., et al. 2013, MNRAS, 436, 540

Harnois-Déraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337

Harnois-Déraps, J., Giblin, B., & Joachimi, B. 2019, A&A, 631, A160

Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623

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

Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146

Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Joachimi, B., Lin, C.-A., Asgari, M., et al. 2021, A&A, 646, A129

Kacprzak, T., & Fluri, J. 2022, Phys. Rev. X, 12, 031029

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

Kilbinger, M. 2015, Rep. Progr. Phys., 78, 086901

Kratochvil, J. M., Lim, E. A., Wang, S., et al. 2012, Phys. Rev. D, 85, 103513

Krause, E., Fang, X., Pandey, S., et al. 2021, ArXiv e-prints [arXiv:2105.13548]

Kruse, G., & Schneider, P. 1999, MNRAS, 302, 821

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]

Li, Z., Liu, J., Zorrilla Matilla, J. M., & Coulton, W. R. 2019, Phys. Rev. D, 99, 063527

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

Lin, C.-A., Kilbinger, M., & Pires, S. 2016, A&A, 593, A88

Liu, J., Petri, A., Haiman, Z., et al. 2015, Phys. Rev. D, 91, 063507P

LSST Science Collaboration (Abell, P. A., et al.) 2009, ArXiv e-prints [arXiv:0912.0201]

Mandelbaum, R., & Hyper Suprime-Cam (HSC) Collaboration 2017, Am. Astron. Soc. Meet. Abstr., 229, 226.02

Mandelbaum, R., Slosar, A., Baldauf, T., et al. 2013, MNRAS, 432, 1544

Martinet, N., Schneider, P., Hildebrandt, H., et al. 2017, MNRAS, 474, 712

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

Parroni, C., Cardone, V. F., Maoli, R., & Scaramella, R. 2020, A&A, 633, A71

Parroni, C., Tollet, É., Cardone, V. F., Maoli, R., & Scaramella, R. 2021, A&A, 645, A123

Peel, A., Lin, C.-A., Lanusse, F., et al. 2017, A&A, 599, A79

Petri, A., May, M., & Haiman, Z. 2016, Phys. Rev. D, 94, 063534

Porredon, A., Crocce, M., Fosalba, P., et al. 2021, Phys. Rev. D, 103, 043503

Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press)

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

Shan, H., Liu, X., Hildebrandt, H., et al. 2017, MNRAS, 474, 1116

Shirasaki, M., Moriwaki, K., Oogi, T., et al. 2021, MNRAS, 504, 1825

Starck, J.-L., & Murtagh, F. 1998, PASP, 110, 193

Starck, J.-L., Fadili, J., & Murtagh, F. 2007, IEEE Signal Proces. Soc., 16, 297

Takada, M., & Jain, B. 2004, MNRAS, 348, 897

To, C., Krause, E., Rozo, E., et al. 2021, Phys. Rev. Lett., 126, 141301

Tutusaus, I., Martinelli, M., Cardone, V. F., et al. 2020, A&A, 643, A70

Vakili, M., Hoekstra, H., Bilicki, M., et al. 2020, A&A, submitted [arXiv:2008.13154]

Zhang, P., Liguori, M., Bean, R., & Dodelson, S. 2007, Phys. Rev. Lett., 99, 141302

Zücher, D., Fluri, J., Ajani, V., et al. 2022a, Towards a Full wCDM Map-based Analysis for Weak Lensing Surveys

Zücher, D., Fluri, J., Sgier, R., et al. 2022b, MNRAS, 511, 2075
## Appendix A: Starlet transform

We provide here a brief introduction to the starlet filter that we employed in this study to filter our input noisy maps. The starlet is a wavelet transform, namely a function that decomposes a signal into a family of scaled and translated functions. More specifically, wavelets are highly localised functions with vanishing mean since they satisfy the admissibility condition $\psi(\int_0^{+\infty} |\psi(k)|^2 \frac{dk}{k} < +\infty)$, which implies that $\int_0^{+\infty} \psi(x)dx = 0$ (i.e. they integrate to zero over their domain). We can define the starlet $\psi$ in relation to its scaling function, which is a B-spline function $\phi$ of order 3,

$$ \psi(t_1, t_2) = 4\phi(2t_1, 2t_2) - \phi(t_1, t_2), \quad (A.1) $$

with

$$ \phi(t) = \frac{1}{12}(|t-2|^3 - 4|t-1|^3 + 6|t|^3 - 4|t+1|^3 + |t+2|^3) \quad (A.2) $$

and $\phi(t, t') = \phi(t)\phi(t')$ (for a complete description and derivation of the starlet transform algorithm, see Starck et al. 2007). We show its 1D and 2D profiles in Fig. A.1. When applying this transform in practice, an original $N \times N$ map $I$ is decomposed into a coarse version of $c_I$ plus several images of the same size at different resolution scales $j$:

$$ I(x, y) = c_I(x, y) + \sum_{j=1}^{j_{\max}} w_j(x, y). \quad (A.3) $$

Here the images $w_j$ represent the details of the original image at dyadic (powers of two) scales, corresponding to a spatial size of $2^j$ pixels and $J = j_{\max} + 1$. Its shape emphasises round features, making it very efficient when dealing with peaks.

<figure><img src="image_7.png" /><figcaption>Fig. A.1. One-dimensional profile of the starlet function.</figcaption></figure>

## Appendix B: Uncertainty on the bias

As the mock catalogues used in the galaxy clustering analysis are produced with a fixed galaxy bias, we provide here an estimate of the impact of the uncertainty of the galaxy bias on the parameter contours through MCMC forecasts. Specifically, we consider the best values and corresponding $\pm 1\sigma$ and $\pm 2\sigma$ found by Vakili et al. (2020) for the sample of luminous red sequence

<table><caption>**Table B.1.** Model constraints and uncertainties on the galaxy bias for the luminous red galaxies from Vakili et al. (2020) derived from the median and the 68% confidence intervals of the marginalised posterior distributions.</caption><thead><tr><th>bin</th><th>-2σ</th><th>-1σ</th><th>best</th><th>+1σ</th><th>+2σ</th></tr></thead><tbody><tr><td>1</td><td>1.36</td><td>1.53</td><td>1.70</td><td>1.88</td><td>2.06</td></tr><tr><td>2</td><td>1.46</td><td>1.59</td><td>1.72</td><td>1.86</td><td>2.00</td></tr><tr><td>3</td><td>1.62</td><td>1.68</td><td>1.74</td><td>1.80</td><td>1.86</td></tr><tr><td>4</td><td>1.85</td><td>1.93</td><td>2.01</td><td>2.09</td><td>2.17</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. B.1.** Forecast contours at confidence levels of 68% and 95% using the power spectrum for four different uncertainties on the galaxy bias (dashed contours) compared to when the galaxy bias is fixed (continuous contours). The light blue dashed contours are obtained using the columns ‘+1σ’ and ‘+2σ’ in Table B.1; the light pink dashed contours using the columns ‘-1σ’ and ‘-2σ’; and the orange continuous contours refer to the column ‘best’.</figcaption></figure>

galaxies for the fourth data release of the Kilo-Degree Survey. We provide in Table B.1 the values corresponding to Table 5 of Vakili et al. (2020).

In order to do so, we perform the inference analysis using the power spectrum computed for four different cases of uncertainty on the galaxy bias. We consider the bias unknown using as mock data the power spectrum computed on simulations with respectively $+2\sigma$, $+1\sigma$, $-2\sigma$, and $-1\sigma$ uncertainty on its value and for inference what was found as best fit (see ‘best’ in Table B.1) for the KiDS-1000 LRG set. We compare the results for these four settings with the case in which the galaxy bias is instead fixed (orange contours in Fig. B.1), namely the same in both mock data and simulations used for inference. We see how the $\pm 1\sigma$ uncertainties in galaxy bias already induce a $2\sigma$ offset in the inferred $\sigma_8$, highlighting the importance of taking into account this uncertainty for real data application. In the future, we plan to perform a deeper investigation of how the uncertainty on the bias impacts the results, for example by emulating and marginalising over the effect.
## Appendix C: Map validation

We validate the convergence and galaxy maps for the two probes that we built from the simulations by comparing our measured two-point statistics with the theoretical prediction. In order to do so we measure the power spectrum respectively on the convergence maps from the source catalogue and the galaxy density map from the lens catalogues built as described in Sect. 3. We show this in the top panel of Fig. C.1 for weak lensing and the bottom panel for galaxy clustering, where the black dots represent the angular power spectrum we measure from the galaxy map we built from the simulations averaged here across the 50 light cones, and the dashed line is the theoretical prediction using pyccl (Chisari et al. 2019). The grey areas indicate the scales excluded in the analysis. Concerning galaxy clustering, as the map is intrinsically noisy (due to the presence of Poisson noise), while the theoretical prediction is obtained for a noiseless case, we compute the power spectrum of the noise and then subtract it from the power spectrum of the map. At high $\ell$ a noisy behaviour can be seen. The scales where this happens are excluded from the analysis, as we include multipoles up to $\ell_{max} = 780$.

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

**Fig. C.1.** Comparison between the power spectrum measured from maps built from the simulations and the corresponding theoretical prediction. **Top panel:** Weak-lensing power spectrum measured from convergence maps built from the simulations and the corresponding theoretical prediction obtained using pyccl. An example of this comparison is shown for the fiducial cosmology of the *cosmo*-SLICS simulations of the third tomographic bin. **Bottom panel:** Same, but for the galaxy density map used for clustering.