Starlet higher order statistics for galaxy clustering and weak lensing

We present a ﬁrst application to photometric galaxy clustering and weak lensing of wavelet-based multi-scale (beyond two points) summary statistics: starlet peak counts and starlet (cid:96) 1 -norm. Peak counts are the local maxima in the map, and (cid:96) 1 -norm is computed via the sum of the absolute values of the starlet (wavelet) decomposition coe ﬃ cients 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 Ω m , the reduced Hubble constant h , the matter ﬂuctuation amplitude σ 8 , and the dark energy equation of state parameter w 0 . In our setting for this ﬁrst application, we consider the two probes to be independent. We ﬁnd that the starlet peaks and the (cid:96) 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


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. 2020), 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 (LSST Science Collaboration et al. 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 × 2pt 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. 2020;Euclid Collaboration et al. 2020;Tutusaus, I. 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, B. et al. 2021;Krause et al. 2021) as they are typically distinct and uncorrelated in different 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 et al. 2023), such as Minkowski functionals (e.g.Kratochvil et al. 2012 andParroni et al. 2020), higher order moments (e.g.Petri et al. 2016 andGatti 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. 2017Martinet et al. 2017, Li et al. 2019, Ajani et al. 2020Zürcher et al. 2022b, Ayçoberry et al. 2023), Betti numbers (Parroni, Carolina 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 andShirasaki et al. 2021), is catching the attention of the community.The 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 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), Ajani et al. (2021), andZürcher 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 Ω m , the reduced Hubble constant h, the matter fluctuation amplitude σ 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.

Weak lensing
The distortions caused by gravitational lensing in the original source image can be summarised in the distortion matrix: Here γ = (γ 1 , γ 2 ) is the shear, also defined in terms of the lensing potential ψ as , which describes the anisotropic distortion of the image; κ is the convergence field, also defined as κ = 1 2 ∇ 2 ψ, which describes the isotropic dilation of the source; and g is the reduced shear.The observed ellipticity obs is related with the intrinsic ellipticity of the galaxy s and the reduced shear through the relation (Seitz & Schneider 1997) where g is the complex conjugate of g.In the weak-lensing regime, where |γ|, κ 1, this relation can be approximated as obs ≈ s + γ.Moreover, if the intrinsic ellipticity of galaxies has no preferred orientation, its expectation value vanishes, s = 0, and the observed ellipticity can be employed as an unbiased estimator of the reduced shear, as obs ≈ γ (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 γ and κ to the lensing potential along with the fact that κ 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) where

Photometric galaxy clustering
As second-order statistics for photometric galaxy clustering, we employed the angular power spectra in harmonic space, C GG ( ), measured on the galaxy maps.Given the redshift distribution of the galaxies n G i (z) in each tomographic bin i of a photometric survey, the average galaxy density in this redshift bin is given by (Euclid Collaboration et al. 2020) Then the radial weight function for a given tomographic bin i for galaxy clustering can be defined as The observable spectrum is theoretically given by the expression where P photo gg ( , z) is the galaxy-galaxy power spectrum, which, under the assumption of a linear galaxy bias, is linked to the matter power spectrum P δδ ( , z) through the bias b g : 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.

Mock data
We modelled our non-Gaussian statistics with mock data constructed from the cosmo-SLICS (Harnois-Déraps, J. ).These were then ray-traced under the Born approximation to construct the sources and lenses catalogues of 100 deg 2 .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.

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 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 (γ 1 , γ 2 ) and the values for the components of the observed ellipticities ( obs 1 , obs 2 ).

Building convergence maps
To get the convergence map, we used the KS inversion: we projected the shear onto a 600 2 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.We provide an example of our validation against the theoretical prediction in Appendix C.

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. 2021, 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 ∆χ = 257.5 Mpc/h.The catalogues are generated assuming linear galaxy bias, whose values and uncertainties are from Vakili et al. (2020) and summarised in Table B.1.

Building galaxy maps
Starting from the lens catalogue described in Sect.3.2, we built galaxy maps for each light cone by assigning the galaxies onto a 600 2 grid.We validated our maps by comparing our measured power spectra against the theoretical prediction obtained with the Python library pyccl.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 max = 780.An example of the comparison between our measurement on the simulated galaxy map and the corresponding theoretical prediction can be found in

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 θ ker = 1 arcmin.For weak lensing we use 33 bins in the range = [180, 2040], which lies within the region where we are confident that the power spectra are in agreement with the theoretical predictions.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 where W(θ ker ) is the Starlet filter, m is the noisy map, and * indicates the convolution operation.The quantity σ filt n 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 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 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.

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: 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, µ is the model prediction as a function of the cosmological parameters θ, and α 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 θ = [Ω m , σ 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 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).
scikit-learn Python package.For more details, the emulator used in this paper is the same as the one employed in Li et al. (2019) and in Ayçoberry et al. (2023). 1

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 1 The code for the GP emulator is publicly available at https:// github.com/CosmoStat/shear-pipe-peaksare computed as where N = 800 is the number of observations, x r i is the value of the i-th data element for a given realisation r, and µ 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 , where n bins is the number of bins and C * the covariance matrix computed for the power spectrum, the peak counts, and the 1 -norm, whose elements are given by Equation ( 9).We area-rescale the covariance matrix to model the statistical accuracy of the KiDS-1000 survey, with A KiDS−1000 = 1000 deg 2 .

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 Equation ( 8).We use 120 walkers initialised in a tiny Gaussian ball of radius 10 −3 around the fiducial cosmology [Ω m , σ 8 , h, w 0 ] = [0.2905,0.8364, 0.6898, −1].

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 (Ω m , σ 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 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 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 reported in Table B.1.The impact is major, where a 1σ shift in galaxy bias offsets the inferred σ 8 by 2σ.To assess then the constraining power of the starlet peaks and the 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 1 -norm both significantly improve the constraints, also in the case of galaxy clustering, with the 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 1 -norm for galaxy clustering and weak lensing.We show these results in Fig. 4 for the weaklensing starlet 1 -norm, the galaxy clustering starlet 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.Fig. 3. Forecast contours at confidence levels of 68% and 95% using the power spectrum (orange), starlet peaks (light red), and the starlet 1 -norm(green) computed on galaxy density maps.Fig. 4. Forecast contours at confidence levels of 68% and 95% using the starlet 1 -norm measured on convergence maps (blue), on galaxy density maps (green), and their combination as independent probes (red).

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 cosmo-SLICS 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 starletbased statistics, such as starlet peak counts and the starlet 1norm 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.  .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'.
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σ, +1σ, −2σ, and −1σ 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 ±1σ uncertainties in galaxy bias already induce a 2σ offset in the inferred σ 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.
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 Ω m in the range [0.10, 0.55], the reduced Hubble constant h in the range [0.60, 0.82], the combination between Ω m and the matter fluctuation amplitude σ 8 , S 8 = σ 8 √ Ω 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 3 particles in a box of comoving side of 505 h −1 Mpc with the Poisson solver cubep 3 m (Harnois-Déraps et al. 2013 Fig. C.1.

Fig. 1 .
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.

Fig
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 TableB.1; the light pink dashed contours using the columns '−1σ' and '−2σ'; and the orange continuous contours refer to the column 'best'.
Ajani et al. (2020)15) 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),Ajani  et al. ( For galaxy clustering, we use 14 bins in the range = [120, 780], following the pessimistic setting of Euclid Collaboration et al. (2020).These regions in are shown in Appendix C. For the non-Gaussian statistics, we filter the noisy maps with a starlet filter (