AMICO galaxy clusters in KiDS-DR3: Cosmological constraints from angular power spectrum and correlation function

We study the tomographic clustering properties of the photometric cluster catalogue derived from the Third Data Release of the Kilo Degree Survey, focusing on the angular correlation function and its spherical harmonic counterpart, the angular power spectrum. We measure the angular correlation function and power spectrum from a sample of 5162 clusters, with an intrinsic richness $\lambda^*\geq 15$, in the photometric redshift range $z\in [0.1, 0.6]$, comparing our measurements with theoretical models, in the framework of the $\Lambda$-Cold Dark Matter cosmology. We perform a Monte Carlo Markov Chain analysis to constrain the cosmological parameters $\Omega_{\mathrm{m}}$, $\sigma_8$ and the structure growth parameter $S_8\equiv\sigma_8 \sqrt{\Omega_{\mathrm{m}}/0.3}$. We adopt Gaussian priors on the parameters of the mass-richness relation, based on the posterior distributions derived from a previous joint analysis of cluster counts and weak lensing mass measurements carried out with the same catalogue. From the angular correlation function, we obtain $\Omega_{\mathrm{m}}=0.32^{+0.05}_{-0.04}$, $\sigma_8=0.77^{+0.13}_{-0.09}$ and $S_8=0.80^{+0.08}_{-0.06}$, in agreement, within $1\sigma$, with 3D clustering result based on the same cluster sample and with existing complementary studies on other datasets. For the angular power spectrum, we derive statistically consistent results, in particular $\Omega_{\mathrm{m}}=0.24^{+0.05}_{-0.04}$ and $S_8=0.93^{+0.11}_{-0.12}$, while the constraint on $\sigma_8$ alone is weaker with respect to the one provided by the angular correlation function, $\sigma_8=1.01^{+0.25}_{-0.17}$. Our results show that the 2D clustering from photometric cluster surveys can provide competitive cosmological constraints with respect to the full 3D clustering statistics, and can be successfully applied to ongoing and forthcoming spectro/photometric surveys.


Introduction
The spatial properties of the large-scale structure (LSS) of the Universe have been recognized as key cosmological probes.According to the Λ−cold dark matter (ΛCDM) model, galaxy clusters are the largest gravitationally bound systems that emerge from the cosmic web of LSS (e.g.Kaiser 1984).They trace peaks in the large-scale matter density field, produced by gravitational infall and hierarchical merging of dark matter haloes (Bardeen et al. 1986;Tormen 1998;Despali et al. 2016).Since their growth is related to the expansion rate of the Universe and to the underlying distribution of matter, cluster statistics represent a powerful tool to understand the structure formation process, to constrain the neutrino mass (e.g.Marulli et al. 2011;Villaescusa-Navarro et al. 2014;Roncarelli et al. 2015), to in-vestigate the nature of dark matter and dark energy (e.g.Mantz et al. 2008;Vikhlinin et al. 2009;Marulli et al. 2012;Sartoris et al. 2016;Costanzi et al. 2019;Moresco et al. 2021;Lesci et al. 2022a,b) and of gravity itself (e.g.Marulli et al. 2021).
Despite the fact that cluster catalogues contain usually a lower number of objects with respect to galaxy catalogues, cosmology with clusters presents a series of key advantages.Galaxy clusters are hosted by the most massive virialised haloes, so they are highly biased tracers, i.e. more clustered than galaxies (e.g.Mo & White 1996;Moscardini et al. 2001;Sheth et al. 2001;Hütsi 2010;Allen et al. 2011;Moresco et al. 2021).Furthermore, thanks to their lower peculiar velocities, galaxy clusters are relatively less affected by nonlinear dynamics at small scales, in particular by the effect of incoherent motions within virialised structures, that generates the so-called fingers-of-God effect (e.g.Veropalumbo et al. 2014;Sereno et al. 2015;Marulli et al. 2017).The impact of redshift-space distorsion (RSD) is thus reduced, allowing us to simplify theoretical assumptions in the modelling of their clustering signal.
During the past years, cluster catalogues have been constructed from observations at several wavelengths (Allen et al. 2011), for example by exploiting the X-ray emission from the diffuse intracluster medium (ICM; Rosati et al. 2002;Böhringer et al. 2004;Pacaud et al. 2016), the millimeter Sunyaev-Zel'dovich effect produced by the inverse Compton scattering between the hot ICM electrons and the cosmic microwave background (CMB) photons (Vanderlinde et al. 2010;Planck Collaboration et al. 2011), and the optical and near-infrared (IR) starlight emission from galaxies (Eisenhardt et al. 2008;Bellagamba et al. 2018).The importance of this multiwavelength approach relies on the possibility to relate different observables, accessible with spectro-photometric observations, to the total mass of clusters, mostly composed of dark matter.The existence of a so-called mass-observable scaling relation (Okabe et al. 2010;Allen et al. 2011;Giodini et al. 2013;Sereno & Ettori 2015;Bellagamba et al. 2019;Sereno et al. 2020;Giocoli et al. 2021) represents a useful link between the theoretical mass function and the distribution of clusters in the space of survey observables, and gives us the opportunity to predict the effective bias of the cluster sample (e.g.Branchini et al. 2017;Lesci et al. 2022b).
In the last decades, the cosmic distribution of the LSS has been investigated in a progressively more accurate and precise way.Typically, measurements of clustering are based on some cosmological assumption for the redshift-distance relation, and require an appropriate reconstruction of the position of cosmic structures, which can be provided by spectroscopic redshift surveys, like the Sloan Digital Sky Survey (SDSS; see Tegmark et al. 2004) or, more recently, the Baryon Oscillation Spectroscopic Survey (BOSS; see Tojeiro et al. 2012) and the Dark Energy Spectroscopic Instrument Legacy Survey (DESI; see Hang et al. 2021).However, spectroscopic surveys are time-consuming thus, in a given amount of observational time, they have a series of limitations in terms of sky coverage and numbers of detected objects.On the other hand, ongoing and future photometric surveys, like the Kilo Degree Survey (KiDS; see de Jong et al. 2017;Kuijken et al. 2019), the Dark Energy Survey (DES; see Dark Energy Survey Collaboration et al. 2016), the Hyper Suprime-Cam (HSC) Subaru Strategic Program (HSC-SSP; see Aihara et al. 2018), the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; see LSST Dark Energy Science Collaboration 2012) and the Euclid mission (Laureijs et al. 2011;Scaramella et al. 2014;Amendola et al. 2018;Euclid Collaboration et al. 2022), will allow us to cover a wider area, imaging also faint sources at high z.
One of the most powerful tools of modern cosmology is the analysis of the two-point correlation function.The simplest and historically first used point-process statistics are the two-point angular correlation function, in configuration space, and its harmonic-space counterpart, the angular power spectrum (Hauser & Peebles 1973;Peebles 1973).In principle, these two statistics bring the same cosmological information, though in practice they have different sensitivities to different scales, due to the finite sizes of real catalogues, and thus the limited range of scales that can be probed.One of their fundamental advantages with respect to the full 3D study is that we can measure the clustering signal from the angular position alone, without any cosmological assumption in converting redshifts to distances (Asorey et al. 2012;Salazar-Albornoz et al. 2014).
The aim of this work is to perform a cosmological analysis based on the catalogue of galaxy clusters identified by the Adaptive Matched Identifier of Clustered Objects (AMICO; see Bellagamba et al. 2018) algorithm from the third data release of the Kilo Degree Survey (KIDS-DR3), presented in Maturi et al. (2019).Here, the availability of photometric redshift measurements allows us to divide the catalogue in shells and to perform a tomographic study, which can provide independent constraints relative to the 3D reconstruction.
The current analysis has been performed with the Cos-moBolognaLib (Marulli et al. 2016) 1 , a set of free software C++ and Python libraries that we used to manage cluster catalogues, to measure their statistical quantities and to perform the Bayesian inferences.
This work is part of a series of papers which aims at exploiting distant clusters in KiDS for both cosmological (Bellagamba et al. 2019;Giocoli et al. 2021;Ingoglia et al. 2022;Lesci et al. 2022a,b;Busillo et al. 2023) and astrophysical studies (Radovich et al. 2020;Puddu et al. 2021).
The paper is organised as follows.In Sect. 2 we present the AMICO KiDS cluster catalogue.In Sects. 3 and 4 we describe the methods used to measure and model the cluster angular correlation function and power spectrum, respectively.In Sect. 5 we discuss the results of the cosmological analysis.Finally, in Sect.6 we draw our conclusions.

Data: the AMICO KiDS-DR3 catalogue
KiDS is a European Southern Observatory (ESO) public optical imaging survey, obtained with the OmegaCAM wide-field imager (Kuijken 2011) mounted on the 2.6m Very Large Telescope (VLT) Survey Telescope (VST), at the Paranal Observatory.This work is focused on its third release, KiDS-DR3 (de Jong et al. 2017).The DR3 covers a total area of 438 deg2 over two fields, one equatorial (KiDS-N) and the other near the South Galactic Pole (KiDS-S), with aperture photometry in u, g, r, and i bands down to the limiting magnitudes of 24.3, 25.1, 24.9, and 23.8, respectively.A final effective area of 377 deg 2 (Maturi et al. 2019), displayed by the pixelated KiDS-DR3 footprint binary mask presented in Fig. 1 (Hildebrandt et al. 2017), which can be found on the KiDS website 2 , was obtained after rejecting all regions affected by satellite tracks, within haloes produced by bright stars and within secondary/tertiary halo masks used for the weak lensing analysis (de Jong et al. 2015;Kuijken et al. 2015).
As discussed in Maturi et al. (2019), clusters were detected thanks to the application of the AMICO algorithm (Bellagamba et al. 2018), which identifies galaxy overdensities by exploiting a linear matched optimal filter.Cluster detection relies only on angular coordinates, magnitudes, and photometric redshifts of galaxies, avoiding a colour-based selection related to the redsequence of clusters.
The complete sample contains 7988 objects, with a signalto-noise ratio (S /N) > 3.5, in the redshift range z ∈ [0.10, 0.80].We limit the current study to z ∈ [0.10, 0.60] because our model is based on the mass-richness scaling relation estimated from the stacked weak-lensing analysis presented in Bellagamba et al. (2019) and Lesci et al. (2022a), which has been calibrated in Fig. 1.The KiDS-DR3 footprint binary mask (Hildebrandt et al. 2017), pixelated with Healpix (Górski et al. 2005), with a resolution given by N side = 512.The top and the middle panels cover the extension of KiDS-N, while the bottom panel refers to KiDS-S.Clusters are distributed within the unmasked area (yellow).Small holes and irregularities reflect the presence of severe incompleteness and substantial photometric degradation, due to satellite tracks, stars etc. this photo-z range.The cluster detection algorithm returns an unbiased redshift estimate with respect to the input photometric catalogue, but it is sensitive to the photo-z bias of the underlying galaxy population, discussed in de Jong et al. ( 2017) for the KiDS survey.Thus, as suggested by Maturi et al. (2019), we correct the estimated cluster redshifts with the relation z = z est − 0.02(1 + z est ).
The mass proxy for the scaling relation used in this paper is the intrinsic richness provided directly by AMICO, defined as the sum of the membership probabilities: where P i ( j) is the probability assigned by AMICO to the i-th galaxy of being a member of a given detection j, r i ( j) is the distance of the i-th galaxy from the centre of the j-th cluster, R 200 (z j ) is the sphere radius in which the mean density is 200 times the critical density at redshift z j , and m * is a function of redshift representing the typical magnitude used in the Schechter function of the cluster model employed by the AMICO algorithm (Maturi et al. 2019).Thus for a given detection the intrinsic richness represents the expected number of visible galaxies, under the condition expressed in Eq. ( 1).According to this definition, it is a nearly redshift-independent quantity because the threshold m * + 1.5 is well below the magnitude limit of the galactic sample, in the considered redshift interval.For our analysis we select clusters with λ * ≥ 15, which ensures a purity higher than 97% and a completeness higher than 50% over the whole sample (Maturi et al. 2019).Finally, to perform a tomographic analysis, we split our catalogue in three different redshift bins.Thinner shells preserve clustering information along the line of sight and are closer to a full 3D study (Asorey et al. 2012;Salazar-Albornoz et al. 2014; Balaguera-Antolínez et al.

The angular correlation function of AMICO KiDS-DR3 catalogue
In this section we describe the methods used to measure and model the angular correlation function.

The angular correlation function estimator
The joint probability of finding two clusters in the solid angle elements δΩ 1 and δΩ 2 , at a distance θ is given by: where n Ω is the mean number of clusters per unit solid angle.Thus, the angular correlation function w(θ) represents the excess probability of finding a pair of objects separated by the angular distance θ, relative to that expected from a random distribution.We measured the observed angular correlation function using the Landy & Szalay (1993) (LS) estimator: where DD(θ), RR(θ), and DR(θ) are the number of data-data, random-random, and data-random pairs in the angular bin θ ± ∆θ/2, respectively.The measurement is performed in eight logarithmically-spaced bins, between 10-20 and 200 arcmin, with a conservative angular separation which takes into account, for the lower limit, the maximum virial cluster size in every redshift bin and, for the upper limit, the angular scale of the survey.
The results are shown for each redshift bin in Fig. 2, where they are also compared to the model presented in Sect.3.2 and to the result of the cosmological analysis of Sect. 5. We construct the random catalogue by randomly extracting the angular (RA, Dec) cluster coordinates within the survey observational tiles, using the same masks adopted in Maturi et al. (2019).To limit shot noise effects, our random catalogue is 30 times larger than the original one.The covariance matrix is estimated through the jackknife method (Norberg et al. 2009).Specifically, for each redshift slice, we project our catalogue onto the celestial sphere, using the equal-area Healpix pixelisation scheme (Górski et al. 2005, see Sect.4.1), with a low-resolution N side = 128, i.e. corresponding to a pixel side of 27 arcmin.Clusters belonging to the same pixel are considered part of a unique region.Therefore, the exact number of regions depends on the quantity of clusters available in each redshift bin and is of the order of 1000.This allows us to estimate the covariance matrix with different N JK measurements of the angular correlation function, obtained after removing one region at a time.

The angular correlation function model
On linear scales, the cluster density field, δ cl (x), is related to the dark matter density field, δ DM (x), through a scale-independent bias, b cl (z): which implies that in Fourier space, P(k) = b 2 cl P DM (k), where n cl (x) is the cluster density, ncl is its mean value and P DM (k) is the linear dark matter power spectrum.We employed the fitting formulae provided by Eisenstein & Hu (1998), which in the angular range of our w(θ) and C ℓ analyses produce consistent results to those provided by CAMB (Lewis et al. 2000) and CLASS (Lesgourgues 2011; Blas et al. 2011).Given the normalised selection function ϕ i (z) in the i-th photometric redshift bin (see Sect. 3.3), we can project the density field onto the celestial sphere, in a given direction n on the sky: The angular correlation function at a given separation θ is the projection of the 3D spatial correlation function, ξ(s): where s = r 2 (z 1 ) + r 2 (z 2 ) − 2r(z 1 )r(z 2 ) cos θ and r(z) is the comoving distance at redshift z.In this work, we assume the planeparallel approximation and we parameterise the linear power spectrum in redshift space as: where f ≡ d ln D d ln a is the linear growth rate, µ is the cosine of the angle between k and the line of sight and b eff is the effective bias, i.e. the halo bias weighted with the halo mass function (see Sect. 3.4).The Fourier transform of the power spectrum gives us the 3D correlation function, which can be expressed in terms of multipoles ξ l (s) and Legendre polynomials P l (µ) (Hamilton 1992): We keep only the monopole because it includes most of the information (Salazar-Albornoz et al. 2014; García-Farieta et al.

Redshift selection function model
Photometric redshifts have larger uncertainties than spectroscopic ones.Because of photo-z errors, cluster photometric redshift distributions can be different from the true redshift distributions, thus we need to account for the conditional probability of having a cluster at the true redshift, z, given the observed redshift, z phot .In Eqs. ( 5) and ( 6) we consider the radial selection function as the normalised cluster distribution in a given redshift bin ∆z i .In other words, it represents the probability of including a cluster in the corresponding photometric shell, depending on the selection characteristics of our study, namely on the binning strategy (Asorey et al. 2012).Our photometric volume-limited survey is selected by the top-hat window function: where z min and z max represent the lower and the upper limits of our photometric redshift interval, respectively.Including objects into redshift shells of a given redshift width ∆z i allows us to 'integrate out' the effect of photo-z errors (Bykov et al. 2023).The conditional probability of detecting a cluster in a sample selected with this window function, i.e. our normalised redshift distribution, is obtained with the following convolution (Budavári et al. 2003;Crocce et al. 2011;Hütsi et al. 2014): where ϕ(z) is the true underlying full redshift distribution: and P(z phot |z) is a Gaussian distribution, as derived by Lesci et al. (2022a,b) from the mock catalogues described in Maturi et al. (2019), whose mean is z, while the standard deviation is equal to: with σ z,0 = 0.02.An estimate of the true redshift distribution can be then obtained as follows: where Ω sky is the survey area in rad 2 .Here, the cosmological dependence is provided by dV dzdΩ , the derivative of the comoving volume, by dn(M,z)  dM , the halo mass function modelled with the functional form by Tinker et al. (2008), and by P(λ * |M, z), where λ * is the intrinsic richness.The latter integral convolves the theoretical mass function, taking into account the shape of the mass-observable scaling relation of the cluster sample (Lesci et al. 2022a): where P(λ * |z) is a cosmology-independent power-law with an exponential cut-off calibrated from mock catalogues (Lesci et al. 2022a), P(M|z) is a normalisation factor computed as the integral of the numerator, P(M|λ * , z) is modelled as a log-normal distribution, in which the mean is given by the mass-observable scaling relation and the root mean square is a free parameter of the model: Here, where E(z) ≡ H(z)/H 0 .We set λ piv = 30 and z piv = 0.35, which represent the central values of intrinsic richness and redshift in the ranges covered by the whole sample, as found by Bellagamba et al. (2019).The intrinsic scatter is modelled with two free parameters, σ intr,0 and σ intr,λ * , as follows: Finally, we need to compute the cluster redshift distribution in a given redshift bin, accounting also for the probability of measuring λ * obs given the true intrinsic richness λ * , and modelling the selection effects and the incompleteness of the sample.This requires a further convolution in the intrinsic richess bin ∆λ * i = ∆λ * (∆z i ), with a Gaussian P(λ * obs |λ * ).We considered an uncertainty of 17% on λ * obs , as in Lesci et al. (2022a,b).In Fig. 3 we show the KiDS-DR3 photometric redshift distribution and the complete number counts model, expressed by:

Effective bias of the cluster sample
Under the assumption that clusters trace the locations of the large-scale dark matter haloes, we adopt a simple linear relation between dark matter and cluster power spectra: b 2 cl = P(k)/P DM (k).While in general the bias is expected to be a function of scale, at large linear scales it can be considered as a scale-independent quantity (e.g.Estrada et al. 2009;Manera & Gaztañaga 2011;Sawangwit et al. 2011).On the other hand, the evolution of both dark matter and cluster clustering leads to a bias which is a function of halo mass and redshift (Estrada et al. 2009).
Following the approach of Lesci et al. (2022b), we assume a constant bias, neglecting any redshift evolution within the broad photometric shell.The effective bias is derived theoretically, as the average over the number of clusters N i in the photometric sample: where b(M, z) is the halo-bias according to the model presented in Tinker et al. (2010).

The angular power spectrum of AMICO KiDS-DR3 catalogue
In this section we describe the methods used to measure and model the angular power spectrum.

The pixelated density map
For every different redshift bin, we generate a cluster density map by projecting the catalogue objects onto the celestial sphere, using the Healpix pixelisation (Górski et al. 2005), with a resolution N side = 512, which ensures a pixel size of approximately 7 arcmin, comparable with the minimum angular scale exploited in the w(θ) analysis.The density contrast δ cl,I in each pixel I is given by: where n cl,I is the cluster number density in the I-th pixel and ncl is its mean, computed in the unmasked area of the survey.The pixelisation procedure smooths information on scales smaller than the pixel size, i.e. ℓ > 1500, which is, in our case, well above the maximum value of ℓ used in our analysis.Different observational effects, like stellar density, air-mass, sky flux and reddening might introduce biases in the galaxy photometry, if not taken into account (Loureiro et al. 2019).This could lead to biases in the cluster detection, thus in the clustering statistics.Anyway, regions affected by these effects are already excluded, since we estimate the cluster density in the unmasked field only, given by the KiDS-DR3 footprint binary mask.Furthermore, as Maturi et al. (2019) imposed the strict magnitude cut at r = 24 (note that the limit in the r band for KiDS-DR3 is 24.9, see Sect. 2), corresponding to the depth of the shallowest tile, the mean cluster count does not depend on the sky position (Lesci et al. 2022b).For these reasons, following for example Branchini et al. (2017), clusters are counted as single objects, thus no weighting scheme has been applied to account for their selection effects.In Fig. 1 we show the KiDS-DR3 footprint mask, degraded from a resolution of N side = 2048 to N side = 512, keeping the 377 deg 2 survey area, which yields to a final unmasked sky fraction of f sky ≈ 0.9%.The irregular geometry of KiDS-DR3 survey is reflected in the large number of small holes contained in the mask.Thus, we avoid any apodisation of the mask, which would lead to a significant modification in the shape and to a non-trivial loss of area (White et al. 2022).

The angular power spectrum estimator
The angular power spectrum of clusters in a given redshift bin i can be measured from the harmonic decomposition of the observed density field.The pixelated density contrast, being defined on a 2D sphere, can be expanded in a series of spherical harmonics: where Y ℓm are the spherical harmonics, computed at the direction on the sky n ≡ (θ, φ), a ℓm are the harmonic coefficients, defined by: the symbol * indicates the complex conjugation operator, while ∆Ω I is the area of the I-th pixel.In this analysis, we used the angular power spectrum estimator introduced by Peebles (1973) and Hauser & Peebles (1973).For a partial sky survey, the masked density contrast is related to the full-sky field through a binary mask function, δcl ( n) = M( n)δ cl ( n), thus the measured pseudo-C ℓ in Fig. 4, named K i j ℓ , is corrected for the sky fraction and defined as follows: where w ℓ is the Healpix pixel window function, which removes the effect of the pixelisation, depending on the parameter N side .
The case i = j refers to the auto power spectrum, while i j to the cross power spectrum between different redshift bins.We removed the shot noise contribution from the measured power spectrum, which accounts for the unclustered part of the measure, given by a discrete distribution of point-like sources.In first approximation, it depends only on the ratio ∆Ω N cl δ i j K , where δ i j K is the Kronecker delta, equal to zero in the cross-correlation case.Possible deviations from this term are further investigated in Sect.4.2.1.Due to the limited size of our cluster catalogues, the shot noise becomes the dominant part of the total signal at As the KiDS-DR3 catalogue does not cover the full-sky, spherical harmonics no longer provide a complete orthonormal basis to expand the angular overdensity field (Camacho et al. 2019).Thus, the measured power spectrum at multipole ℓ depends on an underlying range of multipoles ℓ ′ (Blake et al. 2007).This mode-mode coupling determines a power transfer between different multipoles and it is summarised in R ℓℓ ′ , the socalled mixing matrix, which depends only on the geometry of the angular mask.The ensemble average of the measured power spectrum is related to the theoretical one through (Balaguera-Antolínez et al. 2018): The mixing matrix is equal to an identity matrix in full-sky surveys, where f sky = 1.Starting from the survey window function, R ℓℓ ′ can be expressed in terms of the Wigner 3 j symbols: where: and I ℓm represents the spherical harmonic coefficient of the mask, given by: We estimate the mixing matrix using the publicly available code NaMaster (Alonso et al. 2019), which provides a general framework for the pseudo-C ℓ analysis.The convolution function R ℓ * ℓ ′ is shown in Fig. 5.It is peaked at ℓ * , while it has a drop at different multipoles, depending on the survey area and geometry.In this sense, the mixing matrix gives us an indication about the size of ∆ℓ bands used to bin the measurements.Indeed, with ∆ℓ = 25 we can keep most of the clustering signal, reducing both the effect of the window function, the size of the covariance matrix and the correlation between different bands (Blake et al. 2007;Balaguera-Antolínez et al. 2018;Loureiro et al. 2019).For each bin we compute the weighted average: After the convolution with the mixing matrix, the theoretical C ℓ is averaged with the same weights.As shown in Fig. 4, we restrict our analysis in the range 10 < ℓ < 175.The lower limit reflects the validity of the Limber approximation, while the upper is a conservative value that accounts for the impact of the shot noise (Balaguera-Antolínez et al. 2018), beyond which the clustering signal can be considered negligible.
The analytical error estimation adopted in Blake et al. ( 2007 2019) contains contributions from the cosmic variance and the shot noise, with a boost factor f sky (Blake et al. 2007;Thomas et al. 2011).It is based on the assumption that the coefficients a ℓm are Gaussian distributed, while the effect of the angular window function is modelled by the parameter f sky (Blake et al. 2007;Balaguera-Antolínez et al. 2018), resulting in a diagonal covariance matrix.These approximations do not hold in our case, due to the irregular shape of the KiDS-DR3 survey, which introduces a non-negligible correlation between different multipoles, resulting in an error leak towards other ℓ modes and in a reduction of the diagonal errors (Crocce et al. 2011).Thus we estimate random errors directly from the dataset, using jackknife resampling.In practice, we divided our binary mask in N JK = 400 non-overlapping regions, removing one of them at a time.The cluster density field and the survey area used for the angular power spectrum measurements are then updated, for every realisation, by multiplying the original Healpix map by the new mask.

Shot noise correction
The quantity measured with Eq. ( 24) is the sum of two contributions: the signal and the shot noise.The latter represents the unclustered part of the power spectrum caused by the discreteness of the cluster distribution.The Poisson sampling of pointlike sources contributes to the auto-correlation at null separation, in real space, which brings to a constant power spectrum in harmonic space (Paech et al. 2017), equal to ∆Ω N cl δ i j K .However, when dealing with real data, this simple relation may not always hold.We verify the validity of the analytical relation for the shot noise in two steps.First, we check the power spectrum of 100 density maps derived from random cluster positions within the unmasked regions, considering the same number of cluster as in the real AMICO-KiDS catalogue.We find that the coupling with irregular mask merely increases the dispersion of the shot noise power spectrum around its theoretical prediction, with respect to full-sky surveys, without altering its mean value.Second, we expect deviations from the Poissonian shot noise due to halo exclusion and nonlinear effects (Giocoli et al. 2010;Baldauf et al. 2013).The exclusion simply consists in the fact that clusters cannot overlap, i.e. the distance between clusters cannot be smaller than the sum of their radii, R (Baldauf et al. 2013;Paech et al. 2017).In other words, the cluster-cluster real-space correlation function ξ cc (r) = −1, for r < R, which means that the probability to find another cluster is zero (Baldauf et al. 2013).The exclusion introduces a mass-dependent deviation from the Poissonian shot noise, since objects with different mass occupy different fractions of the sampling volume, while the nonlinear clustering effectively increases the noise term (Paech et al. 2017).
To test the robustness of the Poissonian shot noise hypothesis at the angular scales of our interest, we need a pratical way to disentangle signal and noise.Following Ando et al. (2018), Makiya et al. (2018) and Ibitoye et al. (2022), we first randomly divide the catalogue into two submaps, δ 1,cl and δ 2,cl , both of which contain roughly the same number of clusters.Then we build the half-sum, HS = 1 2 (δ 1,cl + δ 2,cl ), and the half-difference density fluctuation maps, HD = 1 2 (δ 1,cl − δ 2,cl ).By construction, the former contains both signal and noise, while in the latter the signal cancels out, leaving only the shot noise contribution.Since the division of the catalogue into subsets is a random process, the estimated HD map and its power spectrum C HD ℓ slightly change for different realisations (Ando et al. 2018).Thus we averaged over 100 realisations, finding that the shot noise approximation holds in every redshift bin for all the multipoles considered in our analysis.In Fig. 6 we show the ensamble average of C HS ℓ and C HD ℓ power spectra, and the comparison with the Poissonian shot noise.Due to the very low number of clusters in our redshift bins, following Loureiro et al. (2019) we include an extra shot noise term, such that C th,i ∆ℓ → C th,i ∆ℓ + S i .The S i nuisance parameters are forward modelled at the likelihood level, allowing Fig. 6.Shot noise estimation over 100 realisations.Blue dots, black squares and red triangles represent the angular power spectrum of the half-sum (HS) maps, which contain the contribution of both signal and noise, in our three redshift bins (panels from left to right).The purple pentagons show instead the angular power spectrum of the half-difference (HD) maps, which provides a direct estimation of the shot noise, with their average (purple dashed line) and their standard deviation (purple shaded band) in agreement with the theoretical Poissonian value (cyan solid line).
them to vary within a Gaussian prior given by a mean equal to zero and the same standard deviation of C HD ℓ , as shown in Table 1.

The angular power spectrum model
The angular power spectrum is modelled from the spatial power spectrum through a projection kernel, which takes into account redshift evolution and radial selection effects.Its exact computation is given by (Padmanabhan et al. 2007;Thomas et al. 2011;Asorey et al. 2012;Camacho et al. 2019): The kernel function Ψ i ℓ (k) describes the mapping of k to ℓ in real space and is defined as: where ϕ i (z) and b(z) are computed with Eqs. ( 11) and ( 20), respectively, j ℓ (x) is the spherical Bessel function and D(z) is the linear growth factor, normalised such that D(0) = 1 (Camacho et al. 2019).The second term incorporates the linear Kaiser effect (Padmanabhan et al. 2007;Asorey et al. 2012), i.e. the enhancement in 3D power spectrum due to cluster peculiar velocities: Cluster peculiar velocities are smaller with respect to galaxy ones, thus they have a minor effect in our broad redshift selection function, and RSD are erased by the radial projection (Padmanabhan et al. 2007).In particular, for ℓ ≫ 0, the term Ψ i,r tends to zero, so that the total window function reduces to Eq. ( 31) (Padmanabhan et al. 2007;Thomas et al. 2011).However, since the evaluation of spherical Bessel functions is still quite computationally demanding, we make use of the Limber approximation (Limber 1953): where H(z) is the Hubble parameter and the photo-z effects are included through the radial selection function, ϕ(z) (Asorey et al. 2012), see Sect.3.3.Here we underline that P DM (k, z) = P DM (k)D(z) 2 is strictly valid only in linear theory (Blake et al. 2007).
Finally, there are two ways to consider the mode-mode coupling induced by the mask.One can solve the linear system in Eq. ( 25).This requires to bin the pseudo-C ℓ into even more larger bandpowers, since it is computationally expensive and unstable to deconvolve the effect of the mixing matrix from our noisy data (Andrade-Oliveira et al. 2021).For this reason, we decide to include the angular selection directly at the level of the likelihood analysis, choosing the forward modelling (Balaguera-Antolínez et al. 2018;Loureiro et al. 2019;Xavier et al. 2019).
In Fig. 7 we show how redshift-space distortions and partial sky convolution can alter the shape of the angular power spectrum.In particular, the effect of the Limber approximation can be noted only for ℓ ≲ 10, i.e. for multipoles already excluded from our analysis.On the other hand, modifications due to the mixing matrix affect much smaller scales (ℓ ≲ 150), thus we need to properly include them in our model.31) and ( 32) with a dashed line.The convolution with the mixing matrix is shown, respectively, with a solid line (Limber approximation, i.e. the model used in our analysis), and with a dot-dashed line (exact computation).After the mode-mode coupling, the Limber approximation deviates only for ℓ ≲ 10, i.e. for the angular range already excluded from our analysis (gray shaded region).

Likelihood
The analyses of the angular correlation function and power spectrum are performed through Bayesian statistics.We estimate the set of cosmological parameters in Table 1, which enter the model m, by adopting a Gaussian likelihood: in the k-th redshift bin, where: where N is the number of angular bins, µ is the correlation statistic involved, i.e. angular correlation function or power spectrum, and the superscript d and m refer to the quantities obtained from the data and computed with the model, respectively.C −1 i, j,(k) is the inverse of the covariance matrix in the k-th redshift bin, estimated directly from the data for N JK jackknife resamplings (for example Norberg et al. 2009): with expectation values µ i = N res ℓ=1 µ ℓ i /N JK .We neglect correlations between different redshift slices, so that the total likelihood is simply given by the product of the individual likelihoods, L k .

Cosmological results
The Bayesian analysis is performed by adopting uniform priors on both Ω m and σ 8 , while we assume Gaussian priors around the mean values from Planck Collaboration et al. (2020, Table 2, TT, TE and EE+lowE) for the baryon density, Ω b , the primordial spectral index, n s and the normalised Hubble constant, h.They are summarised in Table 1.We use Gaussian priors also for α, β, γ, σ intr,0 and σ intr,λ * , the parameters of the mass-richness scaling relation, whose medians and standard deviations are taken from the posterior distribution of the cluster counts and weak lensing joint analysis, as derived by Lesci et al. (2022a).
In Fig. 8 we show the results of the Monte Carlo Markov Chain (MCMC) cosmological analysis, in the Ω m − σ 8 plane.For the angular correlation function we find Ω m = 0.32 +0.05 −0.04 and σ 8 = 0.77 +0.13  −0.09 as the medians, the 16th and the 84th percentiles of the marginalised 1D posterior distributions.The constraint on the structure growth parameter, S 8 = 0.80 +0.08 −0.06 , is presented in Fig. 9 as our main outcome, and compared with several previous studies available in the literature.In particular, we find a 1σ agreement with cosmological constraints from Wilkinson Microwave Anisotropy Probe (WMAP, Hinshaw et al. 2013, Table 3, WMAP-only, Year 9) and Planck (Planck Collaboration et al. 2020, Table 2, TT, TE, and EE+lowE).An equivalent level of agreement is found also with the number counts analysis presented in Lesci et al. (2022a) using the same AMICO KiDS-DR3 cluster sample, in Costanzi et al. (2019) based on SDSS-DR8 cluster data, in Bocquet et al. (2019) with the 2500 deg 2 South Pole Telescope -Sunyaev-Zel'dovich (SPT-SZ) survey data and with constraints from cosmic shear in DES Year 3 (Amon et al. 2022;Secco et al. 2022), HSC Year 3 (Dalal et al. 2023;Li et al. 2023) and KiDS-DR4 (Asgari et al. 2021).Instead, we note a 1σ tension with S 8 from cluster abundances and weak lensing in DES Year 1 data (Abbott et al. 2020), probably related to richness-dependent effects, since it significantly reduces when their sample is limited to clusters with λ * ≥ 30.
Concerning cluster clustering, our main achievement shows the competitiveness of the 2D correlation function with respect to its 3D counterpart.Indeed, the current study provides more constraining power than the 3D correlation case discussed in Lesci et al. (2022b), where the measurement of ξ(r) is performed within two redshift bins, 0.1 ≤ z ≤ 0.3 and 0.35 ≤ z ≤ 0.6.This is partially due to the slightly larger sample considered, with 228 more clusters, and highlights the importance of the tomographic strategy adopted in photometric redshift surveys, for which several parameters, such as the bin width, the photometric redshift error and the number density of detections, need to be balanced, as discussed in Sect. 2. Our tighter constraints are also confirmed by repeating the full MCMC analysis over 1000 bootstrap resampling, with replacement, which yields to Ω m = 0.33 +0.04 −0.04 , σ 8 = 0.75 +0.11  −0.08 and S 8 = 0.79 +0.08 −0.05 , in excellent agreement with our findings.
The angular power spectrum analysis includes an extra shot noise parameter for each redshift bin, with a Gaussian prior derived with the methodology described in Sect.4.2.1.As shown in Fig. 8 the angular correlation function and the angular power spectrum produce statistically consistent results, although the latter exhibits a much lower constraining power, with wider constraints in particular on σ 8 .This is probably due to the relative importance of the shot noise, which equals the signal contribution even at sub-degree scales and prevents us from extending our angular range to ℓ ≳ 175.We found Ω m = 0.24 +0.05 −0.04 , in excellent agreement with Ω m = 0.24 +0.03 −0.04 presented in Lesci et al. (2022a) with number counts, using the same cluster catalogue, and σ 8 = 1.01 +0.25 −0.17 , which yield to S 8 = 0.93 +0.11 −0.12 , consistent within 1σ with both WMAP and Planck constraints.As shown in Sartoris et al. (2016) and Garrel et al. (2022), the combination of cluster clustering with the more constraining cluster number counts is particularly important since it can highly improve the parameter estimation accuracy both in the ΛCDM Ω m − σ 8 plane and in the dynamical dark energy w 0 − w a plane.The full list of parameters, with prior intervals, the medians, the 16th and the 84th percentiles of the marginalised posterior distributions are shown in Table 1.
Finally, we perform some tests to confirm the robustness of our results.First, we verify that our findings are stable if we adopt the halo mass function parameters provided by Despali et al. (2016) in the model of the cluster redshift distribution.Second, we focus on the inspection of Figs. 2 and 4, which reveals an excess of clustering in the redshift range 0.45 < z ≤ 0.60, common both to w(θ) and C ℓ , with respect to the model median prediction.This tension is marginal at 2σ and has a negligible impact on the overall conclusion.We check our constraints by repeating the entire MCMC analysis, in the first and second redshift bins only.We found consistent results, but with a larger uncertainty due to the lower statistics, since we excluded approximately the 40% of the clusters.In particular, the exclusion of the third bin causes a general broadening of the posterior distributions, with a slight shift for Ω m , to values of 0.25 +0.07 −0.04 and 0.34 +0.07 −0.05 , for angular power spectrum and correlation function, respectively.These results are fully consistent with our expectations.We underline that the excess of clustering appears to be independent of the selection effects, since it remains even if we select more massive clusters with λ * ≥ 20, for which both the purity and the completeness are higher (Maturi et al. 2019).The underestimation in the model might stem from a mismatch between the recovered and the true cluster masses derived from weak lensing calibration, which propagates in the effective bias of the cluster sample.Alternatively, inaccuracies in the cluster selection function, especially at higher redshifts where the cluster redshift distribution deviates more from the theoretical one, could also contribute.However, these systematics are not fundamentally inconsistent with our findings, as we show repeating the analysis without the last redshift bin, and with previous cosmological results obtained in Lesci et al. (2022a,b), using number counts and 3D clustering.We aim at better investigating this aspect by exploiting the larger statistics offered by the KiDS-DR4 (Kuijken et al. 2019) data.

Conclusions
In this paper we presented the cosmological constraints derived from the angular clustering properties of the KiDS-DR3 cluster catalogue.The sample of clusters, which has been constructed with the AMICO algorithm, consists of 5162 galaxy clusters with intrinsic richness λ * ≥ 15.Using a tomographic approach, we measured the angular correlation function and power spectrum in three different photometric redshift bins, z ∈ (0.10, 0.30], z ∈ (0.30, 0.45] and z ∈ (0.45, 0.60], whose widths were selected in order to balance the statistics and the photometric errors.For the angular power spectrum, we verified that the Poissonian shot noise approximation holds in every redshift bin and for all the multipoles considered in our study. We modelled the clustering signal by taking into account the effects of the photometric errors on the redshift selection function, and considering the mass-richness scaling relation from the weak lensing analysis by Lesci et al. (2022a), in order to estimate the effective bias and the redshift distribution of the cluster sample.For the first time, we found cosmological constraints from the angular correlation function and power spectrum of a photometric-redshift cluster catalogue.From the MCMC analyses, we obtained Ω m = 0.32 +0.05 −0.04 , σ 8 = 0.77 +0.13 −0.09 and S 8 = 0.80 +0.08 −0.06 for w(θ), and Ω m = 0.24 +0.05 −0.04 , σ 8 = 1.01 +0.25 −0.17 and S 8 = 0.93 +0.11  −0.12 for C ℓ .Both exhibit a 1σ agreement with the literature results reported in Fig. 8, which includes different cosmological probes like CMB, cluster number counts and cluster clustering.From a comparison with Lesci et al. (2022b), our work has shown that the 2D clustering from a photometric-redshift survey can provide competitive constraints with respect to the full 3D clustering, with the advantage that our findings are cosmological independent, since they rely on the cluster angular positions alone, without any cosmological assumption in converting redshifts to distances.Indeed, from the angular correlation function we derived tighter uncertainties based on the same AM-ICO KiDS-DR3 cluster catalogue, but with a slightly larger sample.This fact reveals the importance of the tomographic strategy adopted in order to fully exploit the cosmological information contained in the cluster catalogue.On the other hand, the angular power spectrum yielded to a wider posterior, in particular with regard to the parameter σ 8 , and does not allow to cover the full angular range explored with w(θ), due to the relative importance of the shot noise.
We tested the the robustness of our study with respect to the parameterisation of the halo mass function presented in Despali et al. (2016).Moreover, we detected an excess of clustering in the redshift range 0.45 < z ≤ 0.60, which does not depend on the selection in richness, since it remains even if we consider only clusters with λ * ≥ 20.However, this tension is marginal at 2σ and does not affect our results, which are stable even if we exclude the third redshift bin, restricting our redshift range to (0.10, 0.45]. We expect more stringent constraints from the analyses of KiDS-DR4 (Kuijken et al. 2019), which covers an area of approximately 1000 deg 2 and includes the photometry of the VISTA Kilo-degree INfrared Galaxy survey (VIKING; see Edge et al. 2013), and of the final KiDS-DR5 (Wright et al. 2023), which will contain data from the full 1350 deg 2 of the KiDS/VIKING footprint.They will allow us to better investigate the differences, the behaviour and the benefits of these two complementary statistics, as well as their combination with number counts and other independent cosmological probes.In the future, we expect an extensive use of the angular clustering of galaxy clusters within the next-generation photometric redshift surveys, like for example Euclid (Laureijs et al. 2011;Scaramella et al. 2014;Amendola et al. 2018;Euclid Collaboration et al. 2022), which will allow us to constrain the parameters of the Dark Energy equation of state, leading to significant advances in the field of the observational cosmology.
2018), but they reduce the projected number density of clusters, and thus the accuracy and the precision of the measurements (Salazar-Albornoz et al. 2014).The redshift bin width is chosen to be five times larger than the maximum photometric error, while we increase the amplitude of the first redshift shell to improve the statistics.The final sample contains 5162 clusters, 1019 in z ∈ (0.10, 0.30], 2072 in z ∈ (0.30, 0.45] and 2071 in z ∈ (0.45, 0.60].

Fig. 2 .
Fig. 2. The angular correlation function measured in three redshift bins: z ∈ (0.10, 0.30] (blue circles), z ∈ (0.30, 0.45] (black squares), z ∈ (0.45, 0.60] (red triangles).Error bars are estimated as the diagonal terms of the jackknife covariance matrix.The dashed lines represent the model computed with cosmological parameters by Planck Collaboration et al. (2020, Table2, TT, TE and EE+lowE).The solid lines show the median of the model distribution computed from the combined posterior of our cosmological analysis, while the shaded regions highlight the 68% confidence levels.

Fig. 3 .
Fig. 3. Normalised redshift distributions.The histograms represent ϕ(z phot ), from the photometric cluster catalogue, while ϕ(z true |W) is predicted from Eq. (11), using the theoretical halo mass function model by Tinker et al. (2008), with cosmological parameters provided by Planck Collaboration et al. (2020, Table 2, TT, TE and EE+lowE), and convolved with our photometric window function, W. The shaded areas indicate the limits of our photometric redshift bins.

Fig. 4 .
Fig. 4. The angular power spectrum measured in three redshift bins: z ∈ (0.10, 0.30] (blue circles), z ∈ (0.30, 0.45] (black squares), z ∈ (0.45, 0.60] (red triangles).Error bars are estimated as the diagonal terms of the jackknife covariance matrix.The dashed lines represent the model computed with cosmological parameters by Planck Collaboration et al. (2020, Table2, TT, TE and EE+lowE).The solid lines show the median of the model distribution computed from the combined posterior of our cosmological analysis, while the shaded regions highlight the 68% confidence levels.

Fig. 5 .
Fig. 5. Normalised elements of the mixing matrix R ℓℓ ′ , centred in three different multipoles ℓ * = 100, 150, 200.The functions decrease as we move away from ℓ * .Their value give us a quantitative amount of the correlation between different multipoles, induced by the mask.The gray shaded regions indicate the bin width ∆ℓ = 25 of our analysis.

Fig. 7 .
Fig. 7. Theoretical angular power spectrum for our three redshift bins.The Limber approximation is represented with a dotted line, while the exact computation presented in Eqs.(31) and (32) with a dashed line.The convolution with the mixing matrix is shown, respectively, with a solid line (Limber approximation, i.e. the model used in our analysis), and with a dot-dashed line (exact computation).After the mode-mode coupling, the Limber approximation deviates only for ℓ ≲ 10, i.e. for the angular range already excluded from our analysis (gray shaded region).

Fig. 8 .
Fig. 8. Cosmological constraints in the Ω m −σ 8 plane, with 68% and 95% confidence intervals, obtained considering angular power spectrum (green contours) and correlation function (purple contours).Our findings are compared to the cosmological constraints derived from Planck Collaboration et al. (2020) (top left) and from the same KiDS-DR3 cluster catalogue (bottom left), in particular to the number counts and 3D correlation function analyses presented inLesci et al. (2022a,b).Right: summary plot containing the 1D marginalised posteriors for the structure growth parameter, S 8 .

Table 1 .
Parameters, prior and posterior mean and percentiles of the cosmological analyses.