© The Authors 2025

1. Introduction

The growth of cosmic structures is governed by long-range gravitational interactions between a large number of systems that behave as tracers of the total matter distribution. Among them, according to the hierarchical model for the cosmic structure formation, galaxy clusters are the latest objects produced by the gravitational infall and the merging of dark matter haloes (Bardeen et al. 1986; Tormen 1998). Reaching masses of up to 10¹⁵ M_⊙ and radii of the order of a few megaparsecs (McNamara & Nulsen 2007), galaxy clusters are the largest gravitationally bound systems ever formed, located at the nodes of the filamentary cosmic web of large-scale structures (LSSs). These astrophysical objects are multi-component systems, made mainly of dark matter and, to a lesser extent, of baryonic matter in different phases (Voit 2005; Allen et al. 2011). This, together with the variety and number of physical processes involved, including both the non-linear and non-local character of the gravitational force, makes a fully analytical treatment of their formation and evolution not accurate enough, considering the precision of current and future datasets.

For this reason, numerical simulations of the LSS are widely used. Since the most relevant effect in the evolution of large-scale density fluctuations is the gravitational interaction, and since the dominant component is the collisionless dark matter, which does not involve computationally expensive hydrodynamical effects, observations of galaxy clusters can be reproduced, within certain limits, by N-body simulations following only the dark matter component. In fact, the baryonic component has only a negligible effect on the formation and evolution of galaxy clusters, and hence on their clustering properties (Springel et al. 2018; Hernández-Aguayo et al. 2023). In particular, N-body simulations rely on a fiducial cosmology, defined by a certain set of cosmological parameters, and evolve the positions of particles undergoing gravitational interaction over cosmic time (Springel et al. 2001; Springel 2005). Therefore, they require a considerable computational effort to realise a sufficient number of mocks necessary for an adequate statistical treatment of the results.

Despite the advances in computational techniques in recent years, there are still some serious limitations to the generation of cosmological simulations (Garaldi et al. 2020; Vogelsberger et al. 2020; Hernández-Aguayo et al. 2023). An accurate description of the evolution of dark matter perturbations is fundamental for modelling the clustering properties of galaxy clusters. First of all, the study of clustering requires an accurate sampling of the baryon acoustic oscillation (BAO) scale, at around 100 Mpc h⁻¹, which is only possible for large simulations (Monaco 2016). Typical requirements for such mock catalogues imply a volume of about 1 Gpc³ and a mass resolution below 10¹⁰ M_⊙ h⁻¹, and become even more prohibitive if we want to produce a wide past light cone (PLC) without multiple replications along the same line of sight (Monaco et al. 2013; Monaco 2016). Furthermore, a large number of simulations corresponding to N independent realisations of the Universe are needed to evaluate the associated errors and covariance matrices (Manera et al. 2013; Euclid Collaboration 2021, 2024a). A smaller number of realisations would make the estimate very noisy (Giocoli et al. 2017, 2020), with relevant consequences for the cosmological constraints, also because of the problems related to the inversion of the covariance matrix (Manera et al. 2013; Munari et al. 2017). Thus, approximate methods are often used to obtain independent realisations in a much faster way, by evolving the density fluctuations with the linear theory, and computing higher-order corrections through a perturbative approach (see e.g. Monaco 2016, for a review).

Even though cluster catalogues usually contain sparse statistics, they can be used extensively in cosmological analyses. In particular, cluster clustering has a strong dependence on the amplitude of the mass power spectrum and the matter density parameter of the Universe (see e.g. Sartoris et al. 2016), so it can be successfully exploited to constrain these cosmological parameters (Veropalumbo et al. 2014, 2016; Sereno et al. 2015; Marulli et al. 2018, 2021; Lindholm et al. 2021; Lesci et al. 2022b; Fumagalli et al. 2024; Romanello et al. 2024). It can also significantly improve the constraining power of other cosmological probes, such as weak gravitational lensing of galaxy clusters (Sereno et al. 2015) and number counts (Sartoris et al. 2016), when analysed in combination with them.

Currently, cluster clustering is less studied than the traditional galaxy clustering; nevertheless, it has a series of remarkable advantages. Because they are hosted by the most massive virialised haloes, galaxy clusters are more clustered than galaxies, meaning that they are more biased tracers (e.g. Mo & White 1996; Moscardini et al. 2001; Sheth et al. 2001; Hütsi 2010; Allen et al. 2011). This reduces the effect of redshift-space distortions (RSDs), allowing us to simplify the theoretical model for the clustering signal. In contrast to the galaxy bias, which is usually considered a nuisance parameter in cosmological analysis because it is difficult to model correctly, especially at small comoving separation, the effective bias of galaxy cluster samples can be predicted theoretically (Branchini et al. 2017; Paech et al. 2017; Lesci et al. 2022b) by means of fitting functions that depend on the mass and redshift of the host haloes (e.g. Tinker et al. 2010). Furthermore, the BAO peak in the correlation function of galaxy clusters exhibits a low non-linear damping and is sharper compared to other tracers, in other words it is closer to the predictions of linear theory (Veropalumbo et al. 2016).

Clustering information is often investigated with the three-dimensional (3D) two-point correlation function, ξ(r), measured from the 3D comoving spatial coordinates of the tracers; however, they are not directly accessible and require a cosmological assumption to convert the observed redshifts and angular positions to distances. Since the fiducial cosmology assumed for the measurements might be different from the true one, the measured two-point correlation function might be geometrically distorted by the Alcock-Paczynski (AP) effect (Alcock & Paczynski 1979). A way to avoid this is to extract the clustering signal from the angular positions alone, using the two-point angular correlation function, w(θ), in configuration space, and its harmonic-space counterpart, the angular power spectrum, C_ℓ (Hauser & Peebles 1973; Peebles 1973), which do not require any AP correction (Asorey et al. 2012; Salazar-Albornoz et al. 2014; Camera et al. 2018). In principle, these two statistics carry the same cosmological information when the entire spectrum of spherical harmonics and angular scales is considered. However, in practice the limited range of wavelengths that can be probed in the reduced volume of real catalogues affects the correlation function and the power spectrum differently (Ata et al. 2018; Avila et al. 2018). Furthermore, we expect a different response to noise and to potential uncorrected observational or systematic effects, which may lead to complementary though highly correlated estimates of cosmological parameters (Sánchez et al. 2017; Camacho et al. 2019; Romanello et al. 2024).

In this paper we make use of a set of dark matter-only simulations to investigate the clustering properties of massive haloes that reproduce some observational conditions of Stage III surveys, such as 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 2016), the Hyper Suprime-Cam³ (HSC) Subaru Strategic Program (HSC-SSP; see Aihara et al. 2018), and the number of clusters expected in Stage IV surveys, for example Euclid⁴ (Laureijs et al. 2011; Euclid Collaboration 2020, 2022b, 2025). In particular, our aim is to explore the best strategies for maximising the cosmological information that can be extracted from angular cluster clustering, in both its variants (i.e. w(θ) and C_ℓ) to verify if it can provide competitive cosmological constraints with respect to the full 3D study. This tomographic approach examines the clustering properties of individual thickness shells at different depths and requires management with an appropriate binning strategy and handling observational effects, such as RSDs and photo-z errors, in order to quantify their impact on the modelling of the clustering signal.

The paper is organised as follows. In Sect. 2 we describe the dark matter-only simulations, and the adopted angular and mass selections. In Sects. 3, 4, and 5 we present the methods used to measure and model ξ(r), w(θ), and C_ℓ. In Sect. 6 we show the result of the Bayesian cosmological analysis. Finally, in Sect. 7 we draw our conclusions and present future perspectives. In Appendix A we quantify the impact of reducing the sky coverage on the angular clustering; in Appendix B we test some simple analytic models for the covariance matrices of the angular correlation function and power spectrum; and in Appendix C we discuss the role of other observational systematics in the estimation of cosmological parameters.

The management of halo catalogues, the measurements of the relevant statistical quantities, the cosmological modelling, and the Bayesian inference presented in this study were performed with COSMOBOLOGNALIB⁵ (Marulli et al. 2016), a set of free software C++ and Python libraries.

2. Data: PINOCCHIO simulations

The PINpointing Orbit-Crossing Collapsed HIerarchical Objects algorithm (PINOCCHIO; Monaco et al. 2002, 2013; Monaco 2016; Munari et al. 2017) is a fast approximated code for the production of dark matter halo catalogues, which works as follows. A Gaussian density field is generated on a cubic grid and then smoothed on a set of scales R. Following the extended Press & Schechter approach (Press & Schechter 1974; Bond et al. 1991), the ellipsoidal collapse model is then used to estimate the collapse time for each particle on the grid, considering all the smoothing scales. The ellipsoid evolution is approximated by the third-order Lagrangian Perturbation Theory (LPT), which is valid until an orbit crossing event occurs or the ellipsoid collapses on its shortest axis. LPT is also used to calculate halo displacements (Monaco 2016; Munari et al. 2017). When the collapse occurs, the particle is expected to be assigned and incorporated into dark matter haloes or in the filaments between them. This fragmentation is done according to the two main processes that characterise the hierarchical clustering of dark matter haloes, namely matter accretion and merging (Monaco et al. 2002, 2013; Munari et al. 2017).

In this way, the code is able to reproduce the mass and the accretion histories of dark matter haloes (Monaco 2016), with an accuracy for the halo mass function and the linear halo bias of 5% and 10%, respectively, compared to the full N-body simulations (Euclid Collaboration 2021, 2024a; Fumagalli et al. 2022). Even an accuracy of 10% in the power spectrum can be achieved up to z ≈ 1. This means that the LPT is able to reproduce the clustering of haloes, placing them in the right position even at relatively low redshifts (Monaco 2016), with a big advantage in terms of computational time. For these reasons PINOCCHIO has been used for several Euclid preparatory science papers (see e.g. Euclid Collaboration 2021, 2024a, 2022a).

2.1. Dark matter halo catalogues: angular and redshift selections

The PINOCCHIO algorithm simulates cubic boxes with side of L = 3870 Mpc h⁻¹ and periodic boundary conditions, from which replication PLCs are produced, selecting only those dark matter haloes that are causally connected to an observer at the present time (Euclid Collaboration 2021, 2024a). We make use of a set of 1000 PLCs, with an aperture of approximately 60 deg, for a final effective area of 10 313 deg² (i.e. one-quarter of the sky). To study the clustering properties of haloes, we focus only on one mock, while the full set of PLCs is used to build the random catalogue (see Sect. 2.2) and to evaluate the numerical covariance matrices. The simulations have been run assuming Planck Collaboration XVI (2014) cosmology: Ω_m = 0.30711, Ω_b = 0.048254, h = 0.6777, n_s = 0.96, σ₈ = 0.8288 and, unless otherwise stated, this is assumed to compute the clustering models that are compared with the measured quantities and presented in the next figures.

The full catalogue contains haloes in the redshift range 0 < z < 2.5, with virial masses M_vir > 6.7 × 10¹³ M_⊙ h⁻¹, sampled with more than 50 particles. However, we focus only on the interval 0.2 < z < 1.0. In real catalogues, a lower limit of 0.1 − 0.2 is expected, due to the reduced lensing power that does not allow a robust lensing analysis, necessary for the mass calibration and the following cosmological exploitation (Bellagamba et al. 2019; Giocoli et al. 2021; Euclid Collaboration 2024b). The upper limit is quite conservative and depends on several considerations. First, it is because our work is in preparation for the analysis of the cluster catalogues form the next KiDS data releases (KiDS-DR4, Kuijken et al. 2019; KiDS-DR5, Wright et al. 2024), where the expected number of clusters at z > 1 is very low and certainly not sufficient for a mass calibration with weak lensing, as already verified in KiDS-DR3, which is extended only up to 0.6 (see e.g. Bellagamba et al. 2019; Giocoli et al. 2021; Lesci et al. 2022a,b; Giocoli et al. 2024; Romanello et al. 2024). Furthermore, the low number of clusters at z > 1 implies that the shot noise term dominates the signal even at the largest angular scales, with non-negligible consequences for the angular power spectrum measurements. Finally, at higher redshifts, we expect a general reduction in the completeness and purity of the cluster sample affecting the robustness of the cosmological analysis.

As explained in Euclid Collaboration (2021), we use a version of the catalogue in which the halo masses have been rescaled in order to match the Despali et al. (2016) halo mass function model, preserving all the fluctuations due to sample variance and shot noise. We also select only haloes with a virial mass larger than 10¹⁴ M_⊙ h⁻¹, which in good approximation corresponds to the predicted limits of the Stage IV photometric survey cluster selection function (i.e. the limiting cluster mass threshold at z < 2, see e.g. Sartoris et al. 2016). The PLCs are available in both real and redshift space. From the latter, we have also produced mock catalogues by introducing a Gaussian redshift perturbation at the redshift z of each halo, which we will hereafter refer to as photo-z space,

$$\begin{array}{c}\hfill P(z_{\mathrm{phot}}|z)=\frac{1}{\sqrt{2\pi }{\sigma }_z}\exp [-\frac{\delta z^2}{2{\sigma }_z^2}],\end{array}$$(1)

where δz = z_phot-z, so that there is no systematic bias in the photometric distribution of the tracers. On the other hand, the standard deviation follows the typical shape for photometric redshift surveys:

$$\begin{array}{c}\hfill {\sigma }_z={\sigma }_{z,0}(1+z).\end{array}$$(2)

The choice of the σ_z, 0 mimics the characteristics of the observations. In particular, we are interested in broad-band photometric data, for which we set σ_z, 0 = 0.05, corresponding to the accuracy required for the Stage IV galaxy photometric surveys (see e.g. Sartoris et al. 2016; Euclid Collaboration 2019). In general, the redshift errors associated with cluster detection are smaller than the corresponding errors in galaxy photometry, because they are computed using the information coming from a set of galaxies. For this reason, we should consider the previous value as a worst-case scenario and also adopt a value of σ_z, 0 = 0.02, which corresponds to the photo-z uncertainty of the KiDS-DR3 cluster catalogue (Maturi et al. 2019). We divide our catalogue either in four redshift bins of binwidth 0.2, or two redshift bins of binwidth 0.4. We limit our analysis to broad redshift bins, several times larger than the photo-z errors, while we do not evaluate the narrow bin case, where Δz ≲ σ_z, because the redshift subsampling dramatically increases the relative importance of the shot noise over the signal, given the limited number of massive haloes. The chosen binwidths give us a minimum of Δz ≈ 2σ_z (for Δz = 0.2, σ_z, 0 = 0.05 and in 0.8 < z < 1.0) and a maximum of Δz ≈ 14σ_z (for Δz = 0.4, σ_z, 0 = 0.02 and in 0.2 < z < 0.6).

2.2. Random catalogue

The construction of the random catalogue is fundamental to avoid biases in the measurement of the 3D and angular correlation functions. The random catalogue must take into account all the possible selection effects associated with the construction of the data catalogue, for example irregularities in the mask due to satellite tracks, stars, or the Milky Way plane. In this case, such effects are limited to the angular and redshift distribution of dark matter haloes within the footprint mask of the simulation. So we build our random catalogue by randomly extracting cluster positions and redshifts from the full set of 1000 halo mocks, a number large enough to avoid the presence of a spurious signal of clustering induced by the coordinate repetition. This process is repeated for real, redshift and photo-z space. In this way the random light cone reproduces the angular positions, average redshift and mass distributions of PINOCCHIO, which are independent of the individual realisations. To limit the shot noise effects, the size of the random catalogue is ten times larger than the halo catalogue.

3. 3D two-point correlation function

3.1. Characterisation of the measurement and covariance matrix

The 3D two-point correlation function is defined as the excess probability of finding a pair of haloes in the comoving volume elements δV₁ and δV₂, at the comoving distance r, relative to that expected from a random distribution. It is expressed by

$$\begin{array}{c}\hfill \delta P_{12}(r)=n_V^2[1+\xi (r)]\delta V_1\delta V_2,\end{array}$$(3)

where n_V is the mean number density of objects per unit comoving volume. The two-point correlation function can be measured with the Landy & Szalay (1993) (LS) estimator,

$$\begin{array}{c}\hfill {\xi }_{\mathrm{LS}}(r)=\frac{\mathrm{DD}(r)+\mathrm{RR}(r)-2\mathrm{DR}(r)}{\mathrm{RR}(r)},\end{array}$$(4)

where DD(r), RR(r), and DR(r) are the number of data–data, random–random and data–random pairs with separation r ± Δr/2, respectively, normalised for the total number of data–data, random–random and data–random pairs. Equation (4) provides an unbiased measurement of the two-point correlation function for a number of random objects N_R → ∞.

In this work, the 3D two-point correlation is measured in ten linearly equispaced radial bins, in the range 15 − 150 Mpc h⁻¹. At these scales the BAO peak is probed and the halo bias can be approximated as a constant (Manera & Gaztañaga 2011; Mandelbaum et al. 2013; Euclid Collaboration 2024a); in other words, it does not depend on r, but only on redshift and mass. The first three even multipole moments of the 3D correlation function at the redshift-space comoving distance s, ξ₀(s), ξ₂(s) and ξ₄(s), that is, the monopole, quadrupole and hexadecapole, respectively, are shown in Fig. 1 for all the redshift bins of our analysis. Focusing on the monopole, we verified that, as expected, the RSDs cause an enhancement of the clustering signal at all scales, with respect to real space, quantifiable by a factor of 1.1 − 1.2, due to the high mass and thus high bias of the haloes. Conversely, the effect of photo-z errors is dominant, leading to a complete deletion of the BAO feature and to a scale-dependent suppression of the clustering signal, which is particularly important at low scales. This leads to a change in the slope of the correlation function with respect to the unperturbed one. The impact of photo-z errors is even more important on the higher-order multipoles of the correlation function. The quadrupole changes sign, providing a positive contribution, while the hexadecapole, which in redshift space is always compatible with zero for s > 20 Mpc h⁻¹, becomes dominant. From the full set of mocks we also measured the numerical covariance matrix:

Fig. 1. First three even multipoles of the 3D two-point correlation function measured from one PINOCCHIO PLC. The errors are computed as the square root of the diagonal elements of the covariance matrix. We show results for monopoles (white circles), quadrupoles (orange squares), and hexadecapoles (blue triangles). The different columns refer to redshift space with σz = 0, and to photo-z space with σz = 0.02(1 + z) and σz = 0.05(1 + z). The solid lines show the theoretical two-point correlation function multipoles, computed from the power spectrum estimated with CAMB, assuming a linear model with ΣNL = 0 (i.e. without non-linear damping effects at the BAO scale). They include systematic effects, such as RSDs and photo-z damping. The colour of the lines corresponds to the multipoles indicated by the symbols.

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{abij}}=\frac{1}{N_{\mathrm{mocks}}-1}{\displaystyle \sum _{m=1}^{N_{\mathrm{mocks}}}}({\xi }_{\mathit{ia}}^m-{\overline{\xi }}_{\mathit{ia}})({\xi }_{\mathit{jb}}^m-{\overline{\xi }}_{\mathit{jb}}).\end{array}$$(5)

Here a and b refer to the spatial bins involved, while i and j mark the considered redshift bins. The square root of the diagonal elements gives us the error bars shown in Fig. 1. By normalising the covariance matrix with its diagonal elements, we can obtain the correlation matrix,

$$\begin{array}{c}\hfill R_{\mathit{abij}}=\frac{{\widehat{C}}_{\mathit{abij}}}{\sqrt{{\widehat{C}}_{\mathit{aaij}}{\widehat{C}}_{\mathit{bbij}}}},\end{array}$$(6)

which is shown in Fig. 2 for real space, redshift space and photo-z space measurements. The cross-correlation between different redshift bins is compatible with zero even in the presence of photo-z errors, which cause an overlap between the redshift distributions of the haloes, thus it is not reported and we only investigate the auto-correlation within the same redshift bin. In particular, if we focus on the block diagonal submatrices, which represent the auto-correlation of the multipoles between different radial bins, we find non-negligible off-diagonal terms that are large even when the bins are widely separated. This is a known feature, reflecting the fact that the 3D and angular correlation functions, in configuration space, are more correlated than their power spectrum counterparts in Fourier and spherical harmonic space (Crocce et al. 2011; Chan et al. 2022). However, this correlation decreases with both increasing multipoles and redshift. Moreover, while we see no significant difference between real and redshift space, due to the limited impact of RSDs, we find a significant reduction in the correlation with larger photo-z errors, which is also common to the w(θ) covariance matrices (see Sect. 4.1).

Fig. 2. Numerical correlation matrices of the 3D two-point correlation function first three even multipoles, derived from 1000 PLCs. The corresponding colour-bar is shown on the right. In the top panels the red diagonal separates the results for real space (upper triangles) and redshift space, with σz = 0 (lower triangles). In the bottom panels the red diagonal separates the results in photo-z space, where RSDs are included, for photo-z errors given by σz = 0.02(1 + z) (upper triangles) and σz = 0.05(1 + z) (lower triangles). In each panel the horizontal and vertical white solid lines separate the correlation sub-matrices of the monopole ξ0(r), the quadrupole ξ2(r), and the hexadecapole ξ4(r), and their corresponding cross-correlations.

On the other hand, from the inspection of the off-block diagonal terms, we can assess the cross-correlation between different multipoles. As expected, in real space, it is consistent with zero for every considered redshift bin, since there is no contribution from quadrupole and hexadecapole. In redshift and photo-z space, the correlation progressively reduces at large z and it is generally higher between adiacent multipoles. In particular, in presence of photo-z errors, we find that positive deviations of ξ₀(s) and ξ₂(s) from their mean value are associated with positive deviations in ξ₂(s) and ξ₄(s), if s₀ < s₂ or s₂ < s₄, respectively, while for s₀ > s₂ and s₂ > s₄ we find an anti-correlation.

3.2. Modelling the 3D two-point correlation function

Given the density contrast field δ(x)≡δρ(x)/ρ_b, where ρ_b is the mean background density, it is possible to define the 3D two-point correlation function,

$$\begin{array}{c}\hfill \xi ({\mathit{x}}_{\mathbf{1}}\mathbf{,}{\mathit{x}}_{\mathbf{2}})=⟨\delta ({\mathit{x}}_{\mathbf{1}})\delta ({\mathit{x}}_{\mathbf{2}})⟩,\end{array}$$(7)

as the expectation value, ⟨…⟩, of the product of the overdensities positioned at points x₁ and x₂. For the cosmological principle we can assume isotropic density fluctuations, so that ξ(x₁, x₂) = ξ(r), which means that it depends only on the distance between the sources, r ≡ |x₂ − x₁|, and not on their positions. Dark matter haloes can be treated as biased tracers of the underlying matter distribution. Therefore, in the range of masses and linear scales that we are analysing this translates into a linear relationship between the halo and the dark matter density contrast, δ_h(x) and δ(x):

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h}}(\mathit{x})=b_{\mathrm{h}}(z)\delta (\mathit{x}),\end{array}$$(8)

through a scale-independent bias, b_h(z), which implies a halo–halo correlation function, ξ_hh(r) = b_h²(z)ξ(r). Since the power spectrum is the Fourier transform of the correlation function, the linear bias can be expressed in the same way, as the ratio between the power spectra: P_hh(k) = b_h²P(k).

For P(k) we adopt the linear model, but taking into account the non-linear damping effects which produce a smearing of the two-point correlation function at the BAO peak, summarised in the parameter Σ_NL, which is left free to vary. Therefore, following Veropalumbo et al. (2016), Avila et al. (2018) and Chan et al. (2018), among others, we modelled P(k) as

$$\begin{array}{c}\hfill P(k)=[P_{\mathrm{lin}}(k)-P_{\mathrm{nw}}(k)]e^{-k^2{\mathrm{\Sigma }}_{\mathrm{NL}}^2/2}+P_{\mathrm{nw}}(k),\end{array}$$(9)

where P_lin(k, z) = P_lin(k, 0)D²(z) is provided by CAMB (Lewis et al. 2000), and is strictly valid only in linear theory (Blake et al. 2007), D(z) is the linear growth factor, normalised so that D(0) = 1, while P_nw(k) is the power spectrum without the BAO peak, as obtained by the parametrisation presented in Eisenstein & Hu (1998). The real-space correlation function, ξ(r), is then obtained by Fourier transforming the corresponding power spectrum.

To compare ξ_hh(r) to the 3D correlation measured from the simulations, we consider the effective bias (i.e. the linear bias integrated above our mass threshold). In practice, b_eff is the halo bias weighted with the halo mass function,

$$\begin{array}{c}\hfill b_{\mathrm{eff}}(\ge M,z)=\frac{{\displaystyle {\int }_{\mathrm{m}}^{\infty }}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M,z)b(M,z)}{{\displaystyle {\int }_{\mathrm{m}}^{\infty }}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M,z)},\end{array}$$(10)

where $\frac{\mathrm{d}n(M,z)}{\mathrm{d}M}$ is computed with the parametrisation provided by Despali et al. (2016), which was also used for the calibration of the halo mass function of the PINOCCHIO catalogues, as mentioned in Sect. 2.1. For b(M, z) we adopt the theoretical halo-bias model presented in Tinker et al. (2010).

3.2.1. Redshift-space and geometric distortions

In real galaxy and cluster photometric surveys, the observed redshifts of the tracers, z_obs, are linked to the cosmological redshifts, z_c, as (neglecting all relativistic corrections and redshift measurement uncertainties)

$$\begin{array}{c}\hfill z_{\mathrm{obs}}=z_{\mathrm{c}}+\frac{v_{\Vert }}{c}(1+z_{\mathrm{c}}),\end{array}$$(11)

where z_c is related to the cosmological Hubble flow, while the second term is due to the line-of-sight peculiar velocities, v_∥, which introduce distortions in the clustering signal if ignored. We neglect the so-called fingers-of-God (FoG) effect, generated by the non-linear dynamics at small scales (Jackson 1972; Peacock et al. 2001), while we consider the large-scale squashing determined by the coherent infall towards collapsing structures, namely the Kaiser effect (Kaiser 1984). In the latter case, the linear regime velocity field can be directly determined from the density field and the anisotropic redshift space power spectrum can be parametrised as

$$\begin{array}{c}\hfill P_{\mathrm{hh}}(k,\mu )=b_{\mathrm{eff}}^2{(1+\frac{f}{b_{\mathrm{eff}}}{\mu }^2)}^{2}P(k),\end{array}$$(12)

where P(k) is computed with Eq. (9), f ≡ dlnD/dlna is the linear growth rate and μ is the cosine of the angle between k and the line of sight. Here, the fμ² term reproduces at all scales the coherent motions of large-scale structure, which acts as a bulk flows from underdense to overdense regions. RSDs in the Kaiser limit are proportional to the parameter β = f/b_eff, thus their effect is less important in more clustered (i.e. in large mass haloes).

The anisotropic redshift–space 3D correlation function, ξ_hh(s, μ), can be expressed in terms of multipoles ξ_ℓ(s) and Legendre polynomials L_ℓ(μ) (Hamilton 1992):

$$\begin{array}{c}\hfill \xi (s,\mu )={\xi }_0(s)+{\xi }_2(s)L_2(\mu )+{\xi }_4(s)L_4(\mu )+\mathcal{O}(s^4).\end{array}$$(13)

To model the multipoles, we can start from the expansion of the anisotropic power spectrum in Eq. (12),

$$\begin{array}{c}\hfill P_{\ell }(k)=\frac{2\ell +1}{2{\alpha }_{\perp }^2{\alpha }_{\Vert }}{\displaystyle {\int }_{-1}^{+1}}\mathrm{d}\mu P(k\prime ,\mu \prime )L_{\ell }\left(\mu \right),\end{array}$$(14)

where μ is the line-of-sight cosine, while k′ and μ′ are rescaled according to the prescription of Beutler et al. (2014)

$$\begin{array}{c}\hfill k^{\prime }=\frac{k}{{\alpha }_{\perp }}{[1+{\mu }^2(\frac{{\alpha }_{\perp }^2}{{\alpha }_{\Vert }^2}-1)]}^{1/2},\end{array}$$(15)

and

$$\begin{array}{c}\hfill {\mu }^{\prime }=\mu \frac{{\alpha }_{\perp }}{{\alpha }_{\Vert }}{[1+{\mu }^2(\frac{{\alpha }_{\perp }^2}{{\alpha }_{\Vert }^2}-1)]}^{-1/2}.\end{array}$$(16)

The geometric distortions are corrected through

$$\begin{array}{cc}\hfill {\alpha }_{\Vert }& =\frac{H^{\mathrm{fid}}(z)}{H(z)}\frac{r_{\mathrm{s}}^{\mathrm{fid}}\left(z_{\mathrm{d}}\right)}{r_{\mathrm{s}}(z_{\mathrm{d}})},\hfill \\ \hfill {\alpha }_{\perp }& =\frac{D_{\mathrm{A}}(z)}{D_{\mathrm{A}}^{\mathrm{fid}}(z)}\frac{r_{\mathrm{s}}^{\mathrm{fid}}\left(z_{\mathrm{d}}\right)}{r_{\mathrm{s}}\left(z_{\mathrm{d}}\right)},\hfill \end{array}$$(17)

where H^fid(z) and $D_{\mathrm{A}}^{\mathrm{fid}}(z)$ represent the fiducial Hubble parameter and angular diameter distance, respectively, while r_s(z_d) is the sound horizon at the drag redshift, z_d, and $r_{\mathrm{s}}^{\mathrm{fid}}(z_{\mathrm{d}})$ is its fiducial value. Then, the multipoles ξ_ℓ are expressed as

$$\begin{array}{c}\hfill {\xi }_{\ell }(s)=\frac{i^{\ell }}{2{\pi }^2}{\displaystyle \int \mathrm{d}}k{k}^2{P}_{\ell }(k)j_{\ell }(ks),\end{array}$$(18)

where i is the imaginary unit and j_ℓ(x) is the ℓth order spherical Bessel function. The monopole can be simply written as a function of the real-space correlation function, ξ(r):

$$\begin{array}{c}\hfill {\xi }_{\mathrm{hh},0}(s)=[b_{\mathrm{eff}}^2+\frac{2}{3}{b}_{\mathrm{eff}}f+\frac{1}{5}{f}^2]\xi (\alpha r).\end{array}$$(19)

Here α is the parameter that corrects for the geometric distortions,

$$\begin{array}{c}\hfill \alpha =\frac{D_V(z)}{D_V^{\mathrm{fid}}(z)}\frac{r_{\mathrm{s}}^{\mathrm{fid}}(z_{\mathrm{d}})}{r_{\mathrm{s}}(z_{\mathrm{d}})},\end{array}$$(20)

where D_V is the average distance at redshift z, and D_V^fid is the same quantity, but estimated at the fiducial cosmology. Their ratio rescales the distances at which the correlation function model is evaluated.

3.2.2. Photo-z errors

Redshift errors introduce a σ_z term in Eq. (11), perturbing the distance measurements along the line of sight:

$$\begin{array}{c}\hfill z_{\mathrm{obs}}=z_{\mathrm{c}}+\frac{v_{\Vert }}{c}(1+z_{\mathrm{c}})\pm {\sigma }_z.\end{array}$$(21)

Following the approach presented in Blake & Bridle (2005), Asorey et al. (2012), Pezzotta et al. (2017), we can modify Eq. (12) as

$$\begin{array}{c}\hfill P_{\mathrm{hh}}(k,\mu )=P(k)b_{\mathrm{eff}}^2{(1+\frac{f}{b_{\mathrm{eff}}}{\mu }^2)}^2\exp (-k^2{\mu }^2{\sigma }^2).\end{array}$$(22)

The exponential factor in Eq. (22) affects the clustering signal in a way similar to that of random peculiar motions (Marulli et al. 2012). It translates in a Gaussian damping term, analogous to the FoG effect (García-Farieta et al. 2020), which causes a scale-dependent suppression of the signal for k > 1/σ, where σ quantifies the impact of photometric random perturbations on cosmological redshifts. In particular,

$$\begin{array}{c}\hfill \sigma =\frac{c{\sigma }_z}{H(z)},\end{array}$$(23)

where c is the speed of light and H(z) is the Hubble parameter. This is equivalent to assuming that the redshift errors follow a Gaussian distribution with zero mean and with dispersion given by σ_z (Hütsi 2010).

The monopole of the 3D correlation function can be computed after integrating over μ, following the redshift-space model presented in Sereno et al. (2015), and used by Lesci et al. (2022b):

$$\begin{array}{c}\hfill {\xi }_{\mathrm{hh},0}(s)=b_{\mathrm{eff}}^2\xi ″(s)+b_{\mathrm{eff}}\xi ″(s)+{\xi }^{″\prime }(s).\end{array}$$(24)

Here different terms of the correlation function are obtained from the inverse Fourier transforms of the corresponding power spectra. In particular,

$$\begin{array}{c}\hfill P″(k)=P(k)\frac{\sqrt{\pi }}{2k\sigma }\mathrm{erf}(k\sigma ),\end{array}$$(25)

$$\begin{array}{c}\hfill P″(k)=\frac{f}{{(k\sigma )}^3}P(k)[\frac{\sqrt{\pi }}{2}\mathrm{erf}(k\sigma )-k\sigma \exp (-k^2{\sigma }^2)],\end{array}$$(26)

and

$$\begin{array}{cc}\hfill P^{″\prime }(k)& =\frac{f^2}{{(k\sigma )}^5}P(k)\{\frac{3\sqrt{\pi }}{8}\mathrm{erf}(k\sigma )\hfill \\ \hfill & \pm \frac{k\sigma }{4}[2{(k\sigma )}^2+3]\exp (-k^2{\sigma }^2)\},\hfill \end{array}$$(27)

respectively. We note that in real space P″(k) and P^″′(k) are both zero. Finally, the multipoles are simply computed by using the power spectrum of Eq. (22) into the Eqs. (14) and (18).

4. Angular two-point correlation function

4.1. Characterisation of the measurement and covariance matrix

The angular correlation function measures the average correlation between any two objects placed in the solid angle elements δΩ₁ and δΩ₂, separated by an angle θ. It can be defined in a completely analogous way with respect to its full 3D counterpart,

$$\begin{array}{c}\hfill \delta P_{12}(\theta )=n_{\mathrm{\Omega }}^2[1+w(\theta )]\delta {\mathrm{\Omega }}_1\delta {\mathrm{\Omega }}_2,\end{array}$$(28)

where n_Ω is the mean number of objects per unit solid angle. We measure w(θ) in ten linearly spaced bins with the LS estimator (Landy & Szalay 1993),

$$\begin{array}{c}\hfill w_{\mathrm{LS}}(\theta )=\frac{DD(\theta )+RR(\theta )-2DR(\theta )}{RR(\theta )},\end{array}$$(29)

where DD(θ), RR(θ), and DR(θ) are the normalised number of data–data, random–random, and data–random pairs in the angular bin θ ± Δθ/2, respectively. The advantage of the angular correlation measurement is that it requires only the random catalogue in RA and Dec, without any information about the redshift distribution. This allows us to greatly simplify the measurement process.

We choose angular ranges that vary as a function of redshift. In particular, we explore scales above the maximum angular size obtained from the virial radii of the halo sample, while for z > 0.4 we adopt a more conservative limit of 20 arcmin, in order to exclude comoving separations below 10 Mpc h⁻¹, where there is a significant deviation of the clustering signal from the linear model prediction, both in ξ(r) and in w(θ). In Fig. 3 we show our results for each redshift bin. For graphical requirements, we plot θ w(θ) and we set independent x-axes. The same angle covers greater physical distances as the redshift increases. This results in a shift of w(θ) towards smaller scales. As for the 3D case, the RSDs have only a slight impact. Again, we find the dominant effect in the photo-z space, where photo-z errors cause a drop in the angular clustering.

Fig. 3. Angular two-point correlation function, measured from one PINOCCHIO PLC. The errors are computed as the square root of the diagonal elements of the covariance matrix. We show results for redshift space with σz = 0 (blue circles), and for photo-z space with RSDs, for photo-z errors equal to σz = 0.02(1 + z) (red squares) and σz = 0.05(1 + z) (green triangles). The solid lines show the theoretical two-point correlation function, computed from the power spectrum estimated with CAMB, assuming a linear model with ΣNL = 0 (i.e. without non-linear damping effects at the BAO scale). They include RSDs, and their colours correspond to the photo-z damping indicated by the symbols.

It is noteworthy that, in contrast to the 3D case, the BAO feature is preserved, suffering only a minor degradation when passing to σ_z, 0 = 0.02 − 0.05. This is because redshift uncertainties affect only the radial distances and not the transverse angular ones. Thus, their effect is greater on the reconstruction of the 3D correlation function or power spectrum, which depend on r or k, rather than on their angular counterparts.

In Fig. 4 we show the numerical covariance matrices, derived from the full set of 1000 PLCs. Their features are similar to the 3D case: no significant level of correlation between different redshift bins, even after the introduction of photo-z errors; positive auto-correlation between different radial bins, which decreases progressively with redshift and with larger photo-z errors.

Fig. 4. Numerical correlation matrices of the angular two-point correlation function derived from 1000 PLCs (see colour-bar on the right). Left: results for real space (upper triangle) and redshift space, with σz = 0 (lower triangle). Right: results in photo-z space, where RSDs are included, for photo-z errors given by σz = 0.02(1 + z) (upper triangle) and σz = 0.05(1 + z) (lower triangle). In each panel the horizontal and vertical white solid lines separate the cross-correlation sub-matrices measured between different redshift bins.

4.2. Modelling the angular two-point correlation function

The study of angular clustering considers the projection of the spatial halo density field onto the celestial sphere, in a given direction of the sky $\widehat{\mathit{n}}$:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h}}^i(\widehat{\mathit{n}})={\displaystyle \int \mathrm{d}z}{\varphi }^i(z){\delta }_{\mathrm{h}}^i(\mathit{x}).\end{array}$$(30)

Here ϕⁱ(z) is the normalised selection function in the ith redshift bin. Given the projected density field, the angular correlation function can be defined as (Peebles 1973)

$$\begin{array}{c}\hfill w_{\mathrm{hh}}(\theta )=⟨{\delta }_{\mathrm{h}}(\widehat{\mathit{n}}){\delta }_{\mathrm{h}}^{\ast }(\widehat{\mathit{n}}+\widehat{\mathit{\theta }})⟩,\end{array}$$(31)

and can be computed as the projection of the 3D spatial correlation function, ξ(s),

$$\begin{array}{c}\hfill w_{\mathrm{hh}}^{\mathit{ij}}(\theta )={\displaystyle \int \int \mathrm{d}}{z}_1\mathrm{d}{z}_2{\varphi }^i(z_1){\varphi }^j(z_2){\xi }_{\mathrm{hh}}(s,\mu ),\end{array}$$(32)

where $s=\sqrt{r^2(z_1)+r^2(z_2)-2r(z_1)r(z_2)cos\theta }$, μ is the same of Eq. (12) and takes the form $\mu =\frac{r(z_2)-r(z_1)}{s}cos(\theta /2)$, and r(z) is the comoving distance at redshift z. The case i = j refers to the auto correlation, while i ≠ j to the cross correlation between different redshift bins. RSDs can be naturally incorporated, by including a multiplicative factor at the level of the power spectrum.

The anisotropic redshift-space 3D correlation function, ξ_hh(s, μ) can be expressed with an infinite series of even multipoles, ξ_ℓ(s), and Legendre polynomials, L_ℓ(μ), such that ξ(s, μ) = ∑_ℓξ_ℓ(s)L_ℓ(μ), as in Eq. (13). Most of the information is contained in the monopole (García-Farieta et al. 2020), but following Salazar-Albornoz et al. (2014), to avoid systematics, we include also the correction given by the quadrupole. We do not account for successive terms, because we have verified that the hexadecapole contribution is negligible, for the adopted P(k) model.

4.3. Redshift selection function model

The radial selection function, ϕ(z), encodes the probability of including a halo in a given redshift bin. In the absence of redshift errors, ϕ(z) reduces to the true underlying redshift distribution,

$$\begin{array}{c}\hfill \varphi (z)=\frac{1}{N}\frac{\mathrm{d}N}{\mathrm{d}z},\end{array}$$(33)

where

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}z}={\mathrm{\Omega }}_{\mathrm{sky}}\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}n(M,z)}{\mathrm{d}M}\mathrm{d}M,\end{array}$$(34)

which is related to the survey area, Ω_sky, in our case 10 313 deg². Equation (34) has a cosmological dependence given by the derivative of the comoving volume, $\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}$, and by the halo mass function, dn(M, z)/dM, for which we use the functional form provided by Despali et al. (2016).

The selection function in a given redshift bin is determined by the window function W(z),

$$\begin{array}{c}\hfill {\varphi }^i(z)=\varphi (z|W)=\varphi (z)W(z),\end{array}$$(35)

which in our redshift-space mock can simply be considered as a top-hat function in the chosen redshift shell. Complications arise when we consider photometric redshifts, which have non-negligible errors that can significantly alter the clustering signal. When dealing with angular clustering, redshift uncertainties enter into the redshift selection function as a modification of the halo redshift distribution, rather than as a damping term at the level of the power spectrum, as we have seen in Sect. 3.2.2 (Asorey et al. 2012). So we need to consider the conditional probability P(z|z_phot) of having a halo at the true redshift, z, given the photometric redshift, z_phot. The selection of an object in a photometric redshift bin is determined by a top-hat window function in z_phot, reflecting the adopted binning strategy (Asorey et al. 2012):

$$\begin{array}{c}\hfill W(z_{\mathrm{phot}})=\{\begin{array}{c}0z_{\mathrm{phot}}\le z_{\min }^i\mathrm{or}{z}_{\mathrm{phot}}>z_{\max }^i\hfill \\ 1z_{\min }^i<z_{\mathrm{phot}}\le z_{\max }^i\hfill \end{array}.\end{array}$$(36)

Here z_min and z_max represent the limits of our photometric redshift interval, respectively.

Thus, the normalised redshift distribution in the i-th photometric bin is obtained with the convolution (Budavári et al. 2003; Crocce et al. 2011; Hütsi et al. 2014)

$$\begin{array}{c}\hfill {\varphi }^i(z)=\varphi (z|W)=\varphi (z){\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}_{\mathrm{phot}}W(z_{\mathrm{phot}})P(z_{\mathrm{phot}}|z),\end{array}$$(37)

where P(z_phot|z) is given by Eq. (1). This quantity can assume a complicated form, depending on several parameters, including the fraction of catastrophic outliers, f_out, namely systems with severely misestimated photometric redshifts (e.g. Hütsi et al. 2014; Euclid Collaboration 2020). For example, in optically selected clusters, if the catastrophic outlier fraction of the member galaxies is large, the cluster properties could be altered, as an increased number of false detections or a large uncertainty in the redshift estimate (Euclid Collaboration 2019). However, in our case, for construction, a simple Gaussian distribution is assumed, whose mean is z, while the standard deviation is equal to σ_z = σ_z, 0(1 + z), with σ_z, 0 = 0 − 0.02 − 0.05, depending on the photo-z errors of our mocks. This choice was also adopted in Romanello et al. (2024) and is motivated by the fact that we do not expect a significant shift in the estimated cluster redshift of real catalogues, assuming the outlier rate of recent photometric surveys (e.g. Hildebrandt et al. 2017). Finally, the redshift distribution in the ith redshift bin is expressed by

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{d{z}^i}={\mathrm{\Omega }}_{\mathrm{sky}}\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}{\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n(M,z)}{\mathrm{d}M}{\displaystyle {\int }_{\mathrm{\Delta }{z}_i}}\mathrm{d}{z}_{\mathrm{phot}}P(z_{\mathrm{phot}}|z).\end{array}$$(38)

In Fig. 5 we show the redshift distribution of haloes, in 0.2 < z < 1.0. For the sake of simplicity, we report only two extreme cases. The coloured, hatched histogram, represented by a solid thick line, shows the real-space redshift distribution, while the dashed histogram shows the photo-z space redshift distribution, in which the redshifts are perturbed with σ_z, 0 = 0.05. Due to this perturbation, some haloes leave their redshift bin, while others enter it. The net effect would be zero in the case of a uniform initial distribution in z. However, given the initial shape of the dN/dz (solid line in Fig. 5), that can be modelled by integrating the Despali et al. (2016) mass function with Eq. (34), between the minimum and maximum masses of the mock catalogue, the introduction of photo-z errors determines a flattening of the overall distribution, leading to a migration of haloes from the most populated to the less populated bins. The result is shown with the partially overlapping dashed curves in Fig. 5.

Fig. 5. Redshift distribution of dark matter haloes, dN/dz, from a single realisation of PINOCCHIO PLC, with an angular radius of approximately 60 deg, in the redshift range 0.2 < z < 1.0. The coloured hatched histogram represented with a solid line shows the real-space redshift distribution, while the dashed histogram shows the photometric distribution, after the introduction of a Gaussian perturbation with σz = 0.05(1 + z). The light shaded regions give the limits of the photometric redshift bins of our tomographic analysis, with Δz = 0.2. The solid line gives the true dN/dz in the window function W, derived by integrating the halo mass function of Despali et al. (2016), assuming the cosmological parameters of the simulation. The dashed lines are derived through Eq. (38) and show the expected halo distribution in the ith redshift bin, selected by W.

5. Angular power spectrum

5.1. HEALPIX pixelated map

Given the angular position and redshift of the haloes, we generate tomographic density maps, using the HEALPIX pixelisation scheme (Górski et al. 2005), with a resolution of N_side = 512. The smoothing of the density field within equal area pixels results in a loss of information for ℓ ≳ 1500, namely the corresponding pixel size. However, this value is well above the upper limits in which the measurements are made (i.e. ℓ = 150 − 200), depending on the redshift bin. In each pixel, the density contrast, δ_h, p, is measured as

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h},p}=\frac{n_{\mathrm{h},p}}{{\overline{n}}_{\mathrm{h}}}-1,\end{array}$$(39)

where n_h, p is the halo count in the pth pixel, while ${\overline{n}}_{\mathrm{h}}$ is the average halo count computed within the circle of 60 degrees angular radius that delimits our PLCs.

5.2. Angular power spectrum estimator

The angular power spectrum can be obtained from the projected density field expressed in Eq. (30). The density contrast, being defined on a 2D sphere, can be expanded into a series of spherical harmonics, Y_ℓm, computed at the direction on the sky $\widehat{\mathit{n}}\equiv (\theta ,\phi )$:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h}}^i(\widehat{\mathit{n}})={\displaystyle \sum _{\ell =0}^{\infty }}{\displaystyle \sum _{m=-\ell }^{m=+\ell }}{a}_{\ell m}^i{Y}_{\ell m}(\widehat{\mathit{n}}).\end{array}$$(40)

Here a_ℓm are the harmonic coefficients. The orthonormality of the Y_ℓm implies that the last equation can be reversed and approximated for a pixellised map:

$$\begin{array}{c}\hfill a_{\ell m}^i={\displaystyle \int \mathrm{d}}\widehat{\mathit{n}}{\delta }_{\mathrm{h}}^i(\widehat{\mathit{n}})Y_{\ell m}^{\ast }(\widehat{\mathit{n}})\simeq {\displaystyle \sum _p^{N_{\mathrm{pix}}}}{\delta }_{\mathrm{h},p}^i{Y}_{\ell m}^{\ast }({\theta }_p,{\phi }_p)\mathrm{\Delta }{\mathrm{\Omega }}_p.\end{array}$$(41)

The angular power spectrum can then be easily evaluated with an estimator already introduced by Hauser & Peebles (1973), Peebles (1973) and widely used in the literature (e.g. Balaguera-Antolínez et al. 2018),

$$\begin{array}{c}\hfill K_{\ell }^{\mathit{ij}}=\frac{1}{w_{\ell }^2}[\frac{1}{f_{\mathrm{sky}}(2\ell +1)}{\displaystyle \sum _{m=-\ell }^{+\ell }}|a_{\ell m}^i{a}_{\ell m}^{\ast j}|-\frac{{\mathrm{\Omega }}_{\mathrm{sky}}}{N_{\mathrm{h}}}{\delta }_{\mathrm{K}}^{\mathit{ij}}],\end{array}$$(42)

where w_ℓ is the HEALPIX pixel window function, which corrects for the loss of power due to the pixelisation, f_sky is the sky fraction covered by our PLC, N_h is the number of haloes in the considered redshift bin and ${\delta }_{\mathrm{K}}^{\mathit{ij}}$ is the Kronecker delta that cancels the shot noise term in the cross-correlation case.

The results are shown in Fig. 6. Coloured symbols refer to the measurements in redshift and photo-z space, while error bars are estimated as the square root of the diagonal of the covariance matrix. RSDs have an impact that increases with redshift, and affect the measured angular power spectrum at large angular scales, corresponding to ℓ < 10 − 40, depending on the redshift bin. However, since we select haloes with large virial masses, their effect is generally negligible in the angular range of our analysis. On the other hand, photo-z errors cause a general suppression of the clustering signal and need to be modelled in the same way as for the angular correlation function. We averaged the C_ℓs in bands of Δℓ = 20, which allowed us to make the covariance matrix more diagonal, since it is blurred due to the partial sky coverage, and to reduce its size (Blake et al. 2007; Balaguera-Antolínez et al. 2018; Loureiro et al. 2019). In Fig. 7 we show the correlation matrices, estimated by substituting C_ℓs in Eq. (5). As expected, assuming that the a_ℓms follow a Gaussian distribution, different multipoles of the power spectrum are independent, and the correlation in spherical harmonics space is lower than in configuration space (see Appendices B.1 and B.2).

Fig. 6. Angular power spectrum, measured from one PINOCCHIO PLC. The errors are computed as the square root of the diagonal elements of the covariance matrix. We show results for redshift space with σz = 0 (blue circles), and photo-z space with RSDs for photo-z errors equal to σz = 0.02(1 + z) (red squares) and σz = 0.05(1 + z) (green triangles). The solid lines show the theoretical angular power spectrum, computed from the power spectrum estimated with CAMB, assuming a linear model with ΣNL = 0 (i.e. without non-linear damping effects at the BAO scale). They include RSDs, and their colours correspond to the photo-z damping indicated by the symbols.

Fig. 7. Numerical correlation matrices of the angular power spectrum, derived from 1000 PLCs (see colour bar on the right). Left: results for real space (upper triangle) and redshift space, with σz = 0 (lower triangle). Right: results in photo-z space where RSDs are included, for photo-z errors given by σz = 0.02(1 + z) (upper triangle) and σz = 0.05(1 + z) (lower triangle). In each panel the horizontal and vertical white solid lines separate the cross-correlation sub-matrices measured between different redshift bins.

5.3. Shot noise correction

The second term in the right-hand side of Eq. (42) quantifies the part of the measured signal related to the discreteness distribution of tracers, which does not contain clustering information. This is called shot noise and must be removed in order to correctly estimate the angular power spectrum. We have seen that the shot noise can be modelled as the ratio between the survey area Ω_sky and the number of haloes N_h in a given redshift bin, while it is equal to zero in the cross-correlation case.

Depending on the number of objects, the shot noise is sensitive to the adopted binning strategy. It plays a minor role in large galaxy catalogues, while it generally becomes important when working with smaller catalogues of clusters. In our case, due to the M > 10¹⁴ M_⊙ h⁻¹ mass cut, we have approximately 10⁴ haloes in each redshift bin, between 0.2 < z < 1.0. This fact enhances the relative importance of the shot noise, which becomes completely dominant over the signal even at relatively large angular scales. For this reason, we set the upper limit of our measurements at the scale where it is about twice as large as the signal, which corresponds to ℓ = 150 − 200, depending on the redshift bin. These values are equivalent to the lower limits θ = 1.2 deg and 0.9 deg, respectively. Thus, the final angular range of the C_ℓ analysis is smaller than that considered for w(θ).

The Poissonian approximation assumes point-like objects and so may not hold for finite-size tracers, such as haloes. We expect mass-dependent deviations from Ω_sky/N_h, due to the halo exclusion, namely the fact that it is not possible for the distance between two haloes to be smaller than the sum of their physical sizes, and non-linear effects (Giocoli et al. 2010; Baldauf et al. 2013; Paech et al. 2017). In Romanello et al. (2024) we verified that this approximation holds for galaxy clusters in the angular range ℓ ∈ [10 − 175], which roughly corresponds to the angular scales we are focusing on. In particular, we use a method already introduced in Ando et al. (2018), Makiya et al. (2018) and Ibitoye et al. (2022), for the analysis of the galaxy-galaxy angular power spectrum. In practice, we randomly split our PLC into two sub-catalogues each containing, by construction, roughly the same number of dark matter haloes. Then we build the corresponding density maps, δ_1, h and δ_2, h, which are slightly different in each realisation, being based on a random extraction process. By combining them, we can build the half-sum, $\mathrm{HS}=\frac{1}{2}({\delta }_{1,\mathrm{h}}+{\delta }_{2,\mathrm{h}})$, and the half-difference density fluctuation maps, $\mathrm{HD}=\frac{1}{2}({\delta }_{1,\mathrm{h}}-{\delta }_{2,\mathrm{h}})$. The half-sum map contains both signal and noise, while only the shot noise contributes to the half-difference map, where the signal, common to both δ_1, h and δ_2, h, cancels out. We estimate C_ℓ^HS and C_ℓ^HD, the half-sum and the half-difference power spectra, respectively, which are then averaged over 100 different realisations. In Fig. 8 we summarise the result. The purple squares show the averaged C_ℓ^HD and the purple shaded regions indicate the stardard deviation of the estimated shot noise. Interestingly, we notice that the theoretical Poissonian approximation, which we indicate with the cyan, horizontal lines, is valid in every redshift bin of our analysis and in a wide range of multipoles, excluding the largest angular scales, namely ℓ ≲ 10, where the measured shot noise is lower than the theoretical expectation. This method shows us that we can simply subtract the Poissonian noise from the measured C_ℓ, as in Eq. (42), and that the uncertainties related to this process can be quantified, and eventually included as a free parameter of the power spectrum model. Following Loureiro et al. (2019) and Romanello et al. (2024), we include in the theoretical model a nuisance parameter, such that C_ℓⁱⁱ → C_ℓⁱⁱ + 𝒮ⁱ. Then, 𝒮ⁱ is forward modelled with a uniform prior centred in zero, and is allowed to vary within a multiple of the standard deviation estimated in the ith redshift bin.

Fig. 8. Shot noise estimation over 100 realisations, using one realisation of the PINOCCHIO PLCs, in four redshift bins in the range 0.2 < z < 1.0 (panels from left to right). The red circles represent the angular power spectrum of the half-sum (HS) maps, which include the contribution of both signal and noise. The purple squares instead show the angular power spectrum of the half-difference (HD) maps, which provide a direct estimate of the shot noise, with their average (purple dashed line) and standard deviation (purple shaded band) in agreement with the theoretical Poissonian value (cyan solid line).

5.4. Modelling the angular power spectrum

The angular power spectrum is modelled from the radial projection of the spatial power spectrum, using the relationship (Padmanabhan et al. 2007; Thomas et al. 2011; Asorey et al. 2012; Camacho et al. 2019)

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{ij}}=\frac{2}{\pi }{\displaystyle \int \mathrm{d}}k{k}^{2}P(k)[{\mathrm{\Psi }}_{\ell }^i(k)+{\mathrm{\Psi }}_{\ell }^{i,r}(k)][{\mathrm{\Psi }}_{\ell }^j(k)+{\mathrm{\Psi }}_{\ell }^{j,r}(k)],\end{array}$$(43)

where the projection kernel, Ψ_ℓⁱ(k), describes the mapping of k to ℓ in real space. As discussed in Sect. 4.3, it includes the radial selection, which is summarised in ϕⁱ(z) (Asorey et al. 2012), see Eq. (37), and the redshift evolution

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\ell }^i(k)={\displaystyle \int \mathrm{d}}zb(z){\varphi }^i(z)D(z)j_{\ell }(kr(z)),\end{array}$$(44)

where j_ℓ(x) is the spherical Bessel function of order ℓ. The second term of Eq. (43) accounts for the enhancement in the 2D power spectrum due to the linear large scale infall motion towards overdense regions (i.e. the Kaiser effect) (Padmanabhan et al. 2007; Asorey et al. 2012) and can be written as

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}^{i,r}(k)& ={\displaystyle \int \mathrm{d}z}{\varphi }^i(z)f(z)D(z)[\frac{2{\ell }^2+2\ell -1}{(2\ell +3)(2\ell -1)}{j}_{\ell }(kr(z))\hfill \\ \hfill & -\frac{\ell (\ell -1)}{(2\ell -1)(2\ell +1)}{j}_{\ell -2}(kr(z))-\frac{(\ell +1)(\ell +2)}{(2\ell +1)(2\ell +3)}{j}_{\ell +2}(kr(z))].\hfill \end{array}$$(45)

Equation (45) derives from the assumption that the magnitudes of the peculiar velocities are small and their impact on z_obs is minor with respect to the radial binwidth (Padmanabhan et al. 2007; Thomas et al. 2011). This is even more accurate considering the peculiar velocities of galaxy clusters, whose effect on the angular power spectrum is limited only to the lower ℓs (Romanello et al. 2024). For ℓ ≫ 0, Ψ^i, r tends to zero (Padmanabhan et al. 2007; Thomas et al. 2011), reflecting the fact that RSDs are erased during the radial integration, and therefore it is not possible to resolve radial perturbations on scales smaller than the thickness of the redshift slices (Padmanabhan et al. 2007).

For the angular scales of our interest, we make use of the Limber approximation (Limber 1953), in order to reduce the computational cost related to the evaluation of spherical Bessel functions. Under this condition, Eq. (43) reduces to

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{ij}}=b_{\mathrm{eff}}^i{b}_{\mathrm{eff}}^j{\displaystyle {\int }_0^{\infty }}\mathrm{d}z{\varphi }^i(z){\varphi }^j(z)P(\frac{\ell +\frac{1}{2}}{r(z)},z)\frac{H(z)}{r^2(z)c},\end{array}$$(46)

where b_eff is computed for each redshift bin with Eq. (10) and P(k) is evaluated with Eq. (9). Strictly speaking this equation is valid for small angles, namely ℓ ≫ 1, but in practice for such massive tracers it coincides with the exact computation expression already at ℓ ≈ 25 (Romanello et al. 2024). As a consequence, the Limber approximation is accurate enough to reproduce the observed power spectra, given the uncertainties on the measurements of the lowest multipoles. Specifically, coloured lines in Fig. 6 are obtained with this simplified model, after the convolution with the mixing matrix (see Eq. (A.4)).

6. Bayesian analysis

We now focus on comparing the posteriors of ξ(s), w(θ) and C_ℓ. The aim is to test whether it is possible to extract consistent amount of cosmological information from these three different probes and to understand to what extent they can be complementary. In particular, we want to investigate if a strategy based on tomographic angular clustering can provide competitive results compared to the full 3D clustering, but with the advantage of not requiring cosmological assumptions to convert from redshifts to distances, during the measurement.

6.1. Priors and likelihood functions

We estimate the set of cosmological parameters the clustering is more sensitive to, namely Ω_m and σ₈. For them we choose large, flat priors, between [0.1, 0.8] and [0.4, 1.6], respectively. The baryon density, Ω_b, the primordial spectral index, n_s, and the normalised Hubble constant, h, are kept fixed to the values of the PINOCCHIO simulation (i.e. from Planck Collaboration XVI 2014). We also choose a flat prior, between [0, 100], for the non-linear damping at the BAO scale, Σ_NL. In the angular power spectrum analysis we have introduced some extra shot noise parameters, 𝒮_k, one for each independent redshift bin, which quantify the variation of the shot noise around the Poissonian value. In order not to lose constraining power, we have chosen narrow prior intervals, namely [ − 3α_k, +3α_k], where α_k is the standard deviation of the shot noise distribution, derived from the study in Sect. 5.3. Its typical value ranges from 2 × 10⁻⁶ to 6 × −10⁻⁶, depending on the binwidth, Δz, and on the considered redshift bin. This choice also takes into account the fact that the Poissonian approximation holds for all the angular scales of our interest, so that the marginalised posterior distributions on 𝒮_k cannot deviate significantly from this value. The full list of parameters is summarised in Table 1.

The clustering analysis is performed through Bayesian statistics by adopting a Gaussian likelihood in the kth redshift bin:

$$\begin{array}{c}\hfill {\mathcal{L}}_k\propto \exp (-{\chi }_k^2/2).\end{array}$$(47)

Here

$$\begin{array}{c}\hfill {\chi }_k^2={\displaystyle \sum _{a=1}^{N_k}}{\displaystyle \sum _{b=1}^{N_k}}{({\mu }_a^d-{\mu }_a^m)}_{(k)}{\widehat{C}}_{\mathit{abkk}}^{-1}{({\mu }_b^d-{\mu }_b^m)}_{(k)},\end{array}$$(48)

where N_k is the number of bins in the kth redshift range, a and b refer to the bin involved, μ is the clustering statistic (i.e. the correlation function or the power spectrum) with superscripts d and m if the quantity is measured from the data or computed with the model, respectively. ${\widehat{C}}_{\mathit{abkk}}^{-1}$ is the inverse of the numerical covariance matrix in the kth redshift bin, estimated directly from 1000 mocks, using ξ(r), w(θ) or C_ℓ in Eq. (5). As already discussed, from an inspection of the blocks outside the diagonal we conclude that the cross-correlation between different redshift bins is consistent with zero; therefore, the total likelihood simply reduces to the product of the individual likelihoods, ℒ_k.

Fig. 9. Marginalised posterior distributions in the Ωm − σ8 plane, with 68% and 95% confidence intervals, for redshift and photo-z space, using a binwidth of Δz = 0.2 and Δz = 0.4, in the upper and lower panels, respectively. The contours are for w(θ) and Cℓ. The dashed lines indicate the Ωm and σ8 values of the PINOCCHIO simulation, while the dot-dashed line shows the predicted S8 degeneracy curve.

6.2. Results

In this section we present the results of the cosmological analyses of ξ(s), w(θ) and C_ℓ. We sampled the posterior distribution with a Monte carlo Markov chain (MCMC), assuming different redshift bins as statistically independent. For all the probes considered, the cosmological parameters of the simulations, corresponding to Planck Collaboration XVI (2014), are correctly recovered, falling within the range at 95% of the corner plot in the Ω_m − σ₈ plane. Systematic effects causing a slight shift of the marginalised posterior distribution can rise due to non-linearities or to the inaccuracy of the halo bias model. The theoretical prediction provided by Tinker et al. (2010) overestimates the numerical bias by about 5%, for low masses and redshifts, see Fig. 2 in Euclid Collaboration (2021). On the other hand, at high masses and redshifts, there is an underestimation of about 10%. However, these scales are excluded from the redshift cut of our sample.

In Fig. 9 we summarise the result of the angular tomographic analysis. First, we note that w(θ) and C_ℓ produce largely overlapping results, giving us the same cosmological information. This is due to the fact that our catalogue has a large coverage of the sky, allowing us to investigate a wide range of different wavelengths. The relative shift of their posteriors, for example, for Δz = 0.2 and σ_z = 0.02(1 + z) or Δz = 0.4 and σ_z = 0.05(1 + z), can be interpreted as an effect of the variance of the specific PLC, and it is not found when instead of measuring on a single mock we consider the average measure from the full set of 1000 PINOCCHIO simulations.

Second, by widening the binwidth from Δz = 0.2 to Δz = 0.4, the marginalised posterior distributions of w(θ) and C_ℓ suffer a broadening on both Ω_m and σ₈. For example, the constraints from the angular correlation function slightly change from $0.{31}_{-0.03}^{+0.03}$ to $0.{32}_{-0.04}^{+0.04}$, and from $0.{80}_{-0.03}^{+0.03}$ to $0.{81}_{-0.03}^{+0.04}$, respectively. With a binwidth of Δz = 0.4 we have both more statistics and less shot noise, as we include a larger number of dark matter haloes. However, the effective bias is considered as a constant, within a given redshift bin, to be taken out of the integrals in Eqs. (32) and (46). This may have systematic consequences for the parameters estimation. Furthermore, the 2D clustering is modified by projection effects that aggregate structures that are not physically bound. In addition, the deeper radial projection erases the BAO peak, because its physical size ends up in a wider range of angular scales. This highlights the importance of defining an appropriate tomographic strategy in clustering studies.

Fig. 10. Marginalised posterior distributions in the Ωm − σ8 plane, with 68% and 95% confidence intervals, for redshift and photo-z space, using a binwidth of Δz = 0.2 and Δz = 0.4, in the upper and lower panels, respectively. The contours are for the different multipoles of ξ(s). The dashed lines indicate the Ωm and σ8 values of the PINOCCHIO simulation, while the dot-dashed line shows the predicted S8 degeneracy curve.

The introduction of photo-z errors, on the other hand, determines a general reduction of the clustering signal, resulting in a loss of constraining power. The posterior distributions given by w(θ) broaden from $0.{32}_{-0.03}^{+0.03}$ to $0.{33}_{-0.04}^{+0.04}$, and from $0.{80}_{-0.03}^{+0.03}$ to $0.{81}_{-0.03}^{+0.03}$, if Δz = 0.2, while they widen from $0.{32}_{-0.04}^{+0.04}$ to $0.{36}_{-0.06}^{+0.05}$, and from $0.{81}_{-0.03}^{+0.04}$ to $0.{79}_{-0.03}^{+0.05}$, if Δz = 0.4, respectively. These results show us that the constraints on σ₈ are generally more stable than the ones on Ω_m.

As we can see in Table 1, we have also included some extra shot noise terms for the angular power spectrum. While in general we have verified that these shot noise parameters are characterised by a flat zero mean posterior, we notice that the Ω_m − σ₈ corner plots, obtained by including these additional free parameters, are slightly wider and better centred.

In Fig. 10 we can see the posterior distributions derived from the 3D correlation function. As in many studies based on real observations, we start from the analysis of the monopole, since it is more stable and less sensitive to the purity of the cluster sample. Then, we progressively take into account the contribution of the quadrupole and the hexadecapole. Focusing on ξ₀(s), we found that enlarging the redshift shell from Δz = 0.2 to Δz = 0.4 does not change significantly the redshift-space marginalised distribution, with a maximum shift of 0.01 in the standard deviations of both Ω_m and σ₈. This result should be interpreted with the fact that in the 3D correlation function the redshift information is not integrated out, thus it contributes to the clustering signal, with only a minor degradation due to the assumption of a constant effective bias in each tomographic bin.

On the other hand, we note that ξ₀(s), in contrast to w(θ) and C_ℓ, exhibits a significant widening of the confidence intervals as the photo-z errors increase. While, from the point of view of the model, the damping of the power spectrum in Eq. (24) is fully equivalent to the introduction of Gaussian errors in the dN/dz, as shown in Hütsi (2010), it is the measurement that is mainly altered, since the photo-z errors significantly modify the radial component of the clustering, leaving the angular part unaffected. For example, while the 3D correlation function suffers from the complete erosion of the BAO peak, this is still visible in w(θ), which only experiences a reduction in the signal due to differences in the radial distribution within a given photometric redshift interval (see Fig. 5). A second consequence is that, in photo-z space, the parameter Σ_NL is not constrained, due to the smearing that affects the BAO peak. Therefore, it behaves as a merely computational expensive complication of the model.

Furthermore, the constraints are better centred on the true cosmology of the simulation as the photo-z errors increase. This is because we start from the assumption of a linear model for P(k), with an infrared (IR) resummation at the BAO peak. As we can see in Fig. 1, this approximation does not hold at low redshift and small scales, because of deviations due to the FoG effect and the non-linear coupling between velocity and density fields (see e.g. Taruya et al. 2010). This could result in a systematic on the estimate of cosmological parameters. Conversely in photo-z space, the effect of non-linearities is diluted by the redshift perturbation, in a similar way to the redshift projection of 2D clustering. This makes it possible to still consider a simple linear power spectrum, avoiding complications in the modelling of the clustering signal.

As evident from Fig. 10, the inclusion of multipoles leads, in general, to a shrinking of the contours along the $S_8\equiv {\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ degeneration line. We should notice, in particular, that in redshift space the clustering signal comes from the monopole and the quadrupole only, while the contribution of the hexadecapole is negligible, as we can see in Fig. 1, being consistent with zero for s > 20 Mpc h⁻¹ in every bin of our analysis. On the other hand, in photo-z space, the inclusion of ξ₂(s) and ξ₄(s) is relevant. By way of example, for Δz = 0.4, the constraints on Ω_m and σ₈ change from $0.{35}_{-0.06}^{+0.08}$ and $0.{83}_{-0.07}^{+0.10}$, to $0.{34}_{-0.04}^{+0.04}$ and $0.{80}_{-0.05}^{+0.06}$, and finally to $0.{35}_{-0.03}^{+0.03}$ and $0.{78}_{-0.03}^{+0.04}$, as we add the even multipoles.

In Fig. 11, in order to summarise the cosmological results and to understand what contribution cluster clustering can make, we focus on the structure growth parameter, S₈. It is particularly relevant in cosmology, since various measurements in the local Universe (e.g. cosmic shear, galaxy clustering, galaxy and galaxy cluster number counts) have shown a 2 − 3σ tension with the high-z value inferred by studying the CMB anisotropies from Planck (see e.g. Bocquet et al. 2019; Asgari et al. 2021; Heymans et al. 2021; Abbott et al. 2022; Amon et al. 2022), which represents a potential stress for the current ΛCDM model.

Fig. 11. Summary plot for the parameter S8, in redshift and photo-z space with σz = 0.02(1 + z) and σz = 0.05(1 + z), panels from left to right, for w(θ) and ξ0(s). The solid lines refer to Δz = 0.2, the dashed lines to Δz = 0.4. The dot-dashed line indicates the S8 value of the PINOCCHIO simulation.

For the sake of simplicity, we compare only the two-point correlation function in both its variants (i.e. w(θ) and ξ(s)). These two probes lie in configuration space. Moreover, they are built from the same random catalogue, and therefore are subject to common systematics. While in redshift space their competitiveness essentially depends on the width of the analysed redshift bins, with a slight decentring of ξ₀(s) due to non-linearities at small scales, we note that in the presence of a large photo-z errors, such as 0.05(1 + z), w(θ) is always favoured over ξ₀(s), with $0.{85}_{-0.03}^{+0.03}$ versus $0.{87}_{-0.04}^{+0.06}$ and $0.{86}_{-0.03}^{+0.04}$, versus $0.{86}_{-0.05}^{+0.07}$, for Δz = 0.2 and Δz = 0.4, respectively. The latter requires the inclusion of the multipoles to produce results comparable or superior to the 2D clustering, depending on Δz. In this case, more care should be taken to avoid systematics due to non-linearities not included in the model, at scales s < 20 Mpc h⁻¹. The intermediate case, with σ_z = 0.02(1 + z), shows that the 2D correlation function is clearly preferred for Δz = 0.2, while it is slightly disfavoured for Δz = 0.4, due to the lack of the redshift information. This fact confirms the possible advantages and the competitiveness of the angular tomographic clustering within photometric redshift surveys. In particular, the potential and limitations of this approach have already been tested in the angular clustering analysis of the KiDS-DR3 cluster catalogue, discussed in Romanello et al. (2024), where the constraint on S₈ are mildly more stringent than the corresponding 3D study (Lesci et al. 2022b).

7. Conclusions

In this paper we analysed the clustering properties of dark matter haloes from the PINOCCHIO simulation in preparation for subsequent work to be carried out on forthcoming galaxy cluster catalogues from Stage III and Stage IV photometric surveys. We used a PLC with an angular size of 60 deg, thus covering a circular area of 10 313 deg², namely one-quarter of the sky. We selected dark matter haloes with a virial mass greater than 10¹⁴ M_⊙ h⁻¹ and redshift in the range 0.2 < z < 1.0. The mass selection reproduces to a good approximation what is expected for future cluster catalogues identified in the Stage IV photometric surveys, while the selection in z excludes the low-redshift regime, where it will be difficult to perform a robust lensing analysis for the mass calibration of real galaxy clusters, and the high-redshift regime, where the clustering signal is completely dominated by shot noise.

Our study used a tomographic approach adopting the binwidths Δz = 0.2 − 0.4 (i.e. four or two redshift bins) to compare 3D and 2D clustering. The former was analysed by means of the 3D two-point correlation function, while for the latter we made use of the angular correlation function and its spherical harmonic counterpart, the angular power spectrum. For each cosmological probe, we estimated the numerical covariance matrix from a set of 1000 mock catalogues, in both redshift space and photo-z space. For this purpose we introduced Gaussian errors on the halo redshifts, with zero mean and standard deviation given by σ_z = σ_z, 0(1 + z), with the typical σ_z, 0 of photometric redshift surveys, namely 0.02 and 0.05.

We used a linear model for the power spectrum, including the IR-resummation, which characterises non-linear damping phenomena at the BAO scale. The 3D correlation function was obtained via Fourier transform of the corresponding P(k), computing the effective bias of dark matter haloes by assuming the model presented in Tinker et al. (2010) for the halo bias, and the parametrisation of Despali et al. (2016) for the halo mass function. The angular correlation function and power spectrum are modelled from the radial projection of the corresponding 3D correlation function and power spectrum, respectively, with a selection kernel that takes into account the normalised redshift distribution of haloes. RSDs are modelled in the Kaiser limit for the correlation functions, while they are not considered in the C_ℓ as they only slightly affect the largest angular scales, being fully subdominant relative to the measurement uncertainties. We found that they have only a minor effect on ξ₀(s) and on w(θ), quantified in a boost factor of 1.1 - 1.2, due to the high bias of the tracers. On the other hand, photo-z errors play the major role in suppressing the clustering signal and the numerical errors. In particular, they modify the small-scale slope of the 3D correlation function monopole, and completely erase the BAO features, which conversely remain visible in w(θ). Moreover, photo-z errors determine a substantial increase in the contribution of the quadrupole and the hexadecapole. In particular, ξ₂(s) changes its sign, moving from redshift space to photo-z space, while ξ₄(s), which is consistent with zero in redshift space, becomes dominant for σ_z ≳ 0.02(1 + z). For the angular power spectrum, we verified that the Limber approximation and the Poissonian assumption for the shot noise hold over the full redshift interval and angular scales considered in our analysis.

Finally, we performed a Bayesian MCMC analysis to constrain fundamental cosmological parameters, such as Ω_m and σ₈. Even though all cosmological probes are able to accurately estimate the cosmological parameters of the simulations, we find that the posterior distributions obtained from the 3D correlation function degrade more in the presence of photo-z errors than from its angular counterparts. In the presence of large photo-z errors, ξ₀(s) produces wider constraints on the structure growth parameter, S₈, with respect to w(θ). Conversely, by also exploiting the redshift infomation, it is less affected by the specific tomographic redshift binning strategy adopted, but is more sensitive to non-linearities that characterise the small-scale clustering. The inclusion of the two-point correlation function multipoles determines a collapse of the Ω_m − σ₈ contours along the S₈ degeneracy curve, which is more relevant in photo-z space, while the contribution of ξ₄(s) is negligible in redshift space.

On the other hand, w(θ) and C_ℓ provide the same cosmological information, and unbiased cosmological constraints, which in turn complement that obtained from the 3D analysis, with the advantage of being able to rely on a linear model for P(k) since non-linearities are diluted by projection effects. Moreover, in the presence of photo-z errors, for σ_z, 0 ≳ 0.02, the angular clustering shows its full potential as a cosmological probe, with respect to the traditional 3D study, with the advantage of not requiring a fiducial cosmology to convert redshifts into distances in the measurements. This comparison was already performed in Romanello et al. (2024), where we successfully applied the models developed in this paper to the AMICO KiDS-DR3 galaxy cluster catalogue, finding consistent results with respect to the ξ₀(s) analysis, with slightly more stringent constraints on S₈. This highlights the significance of extending the analysis of 2D clustering to other cosmological parameters, for example studying the degeneracy in the w₀ − w_a dynamical dark energy plane and the galaxy cluster catalogues that will be available from ongoing Stage III and Stage IV photometric surveys, such as the forthcoming data releases of KiDS and Euclid.

Acknowledgments

We acknowledge the financial contribution from the grant PRIN-MUR 2022 20227RNLY3 “The concordance cosmological model: stress-tests with galaxy clusters” supported by Next Generation EU and from the grants ASI n.2018-23-HH.0 and n. 2024-10-HH.0 “Attivitá scientifiche per la missione Euclid – fase E”, and the use of computational resources from the parallel computing cluster of the Open Physics Hub (https://site.unibo.it/openphysicshub/en) at the Physics and Astronomy Department in Bologna. We would like to thank P. Monaco for the production of the PINOCCHIO mocks, R. Zangarelli and the anonymous referee for the constructive suggestions that improved this work.

References

Appendix A: Reducing the sky coverage

Fig. A.1. Mixing matrices for different sky coverages. Left panels: Colour plots of the mixing matrix, Rℓℓ′, computed up to ℓ = 500 with the public code NAMASTER (Alonso et al. 2019). Right panels: Plots of Rℓℓ*/Rℓ*ℓ*, centred in three different multipoles, at ℓ = 100 (blue), 150 (orange) and 200 (green). They quantify the mode-mode coupling related to the partial sky coverage, reflecting the blurring of the mixing matrix along the diagonal. The grey shaded regions indicate the multipole range of Δℓ in which we bin our measurements. The top panels are computed from the 10313 deg2 region covered by the full PINOCCHIO footprint mask, while at the bottom we focus on the 1004 deg2 subregion.

In partial-sky surveys, spherical harmonics no longer provide a complete orthonormal basis for the expansion of the angular overdensities (Camacho et al. 2019). The masked density field is related to its full-sky counterpart through a binary mask function, ${\stackrel{\sim }{\delta }}_{\mathrm{h}}(\widehat{\mathit{n}})=M(\widehat{\mathit{n}}){\delta }_{\mathrm{h}}(\widehat{\mathit{n}})$, which can be considered as an angular window function. Its effect propagates into the power spectrum evaluation, in the form of a coupling between different multipoles, which becomes progressively more important for lower sky coverages and survey fragmentation. In particular, the measured power spectrum at multipole ℓ depends on an underlying range of multipoles ℓ′ (Blake et al. 2007). The power transfer between different multipoles is described by the mixing matrix, R_ℓℓ′, which depends only on the geometry of the angular mask. It reduces to the identity matrix in full-sky surveys and can be expressed in terms of the Wigner 3j symbols

$$\begin{array}{c}\hfill R_{\ell {\ell }^{\prime }}=\frac{2{\ell }^{\prime }+1}{4\pi }{\displaystyle \sum _{{\ell }^{″}}}(2{\ell }^{″}+1)W_{{\ell }^{″}}{\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& {\ell }^{″}\\ 0& 0& 0\end{array}\right)}^2,\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill W_{\ell }={\displaystyle \sum _{m=-\ell }^{+\ell }}\frac{|I_{\ell m}{|}^2}{(2\ell +1)}.\end{array}$$(A.2)

Here I_ℓm represents the spherical harmonic coefficient of the mask, which can be approximated in a pixelated sky as

$$\begin{array}{c}\hfill I_{\ell m}\simeq {\displaystyle \sum _p^{N_{\mathrm{pix}}}}{Y}_{\ell m}^{\ast }({\theta }_p,{\phi }_p)\mathrm{\Delta }{\mathrm{\Omega }}_p.\end{array}$$(A.3)

With the mixing matrix we can link the ensemble average of the measured power spectrum with the theoretical power spectrum through (Balaguera-Antolínez et al. 2018)

$$\begin{array}{c}\hfill ⟨K_{\ell }^{\mathit{ij}}⟩=\frac{1}{f_{\mathrm{sky}}}{\displaystyle \sum _{{\ell }^{\prime }}}{R}_{\ell {\ell }^{\prime }}{C}_{\ell }.\end{array}$$(A.4)

We estimate the mixing matrix of the footprint masks using the publicly available code NAMASTER⁶ (Alonso et al. 2019), which provides a general framework for the pseudo-C_ℓ analysis.

To better reproduce the observational conditions of photometric surveys, which often scan separated parts of the sky, we produced an angular mask composed of three different patches, two near the equator (P1 and P2), and another further south, (P3). Their angular limits are given respectively by

• P1: RA ∈[−60, −30], Dec ∈[−10, 0];

• P2: RA ∈[0, 30], Dec ∈[4, 10];

• P3: RA ∈[−30, 30], Dec ∈[−40, −30].

The final effective area is equal to 1004 deg², which simulates the expected sky coverage of the cluster catalogue of KiDS-DR4, for which we plan to carry out a subsequent in-depth study.

Reducing the sky coverage has important consequences on our analysis. First, it increases the angular auto-correlation, with a blurring of the correlation matrix both in configuration and in spherical harmonic space. At the level of the angular power spectrum, we find that the mixing matrix calculated over an area of 1004 deg² is more blurred than the one computed over the full PLC, as we can see in Fig. A.1, with the convolution function R_ℓ^*ℓ′. The grey shaded regions highlight the Δℓ = 8 bandwidth, which allows us to include almost all of the clustering signal, reducing the power transfer between different ℓ bands due to the off-diagonal terms, while with the reduced sky coverage, the same band of Δℓ = 8 identifies only 50% of the transferred power.

Reducing the survey area, on the other hand, does not change the shot noise, since both Ω_sky and N_h decrease in the same way. However, a smaller range of wavelengths can be investigated within the reduced volume. This fact suggests that it is possible to exploit C_ℓ and w(θ) as complementary sources of cosmological information and specifically, to use their comparison to assess possible model systematics.

Appendix B: Analytical modelling of the covariance matrix

A fundamental aspect of cluster clustering cosmology is the correct evaluation of measurement errors. Biases in the determination of cosmological parameters can be present, if all significant observational systematic uncertainties are not properly taken into account.

The estimation of covariance matrices, which summarise the set of uncertainties associated with the observables, is therefore a fundamental challenge for the future. Analytical covariance matrix calculations are probably the best solution, as they do not require large computational resources. In addition, they are noise-free and naturally include the dependence on cosmological parameters (Euclid Collaboration 2024a). However, the accuracy of analytical models might be limited due to the complex physics of structure formation, which introduces non-linearities and non-Gaussianities in clustering, and by the necessity of appropriate models for shot noise and RSDs (Manera et al. 2013; Euclid Collaboration 2024a). These and other observational effects, such as angular masks and redshift uncertainties, can lead to systematic theoretical errors, and therefore their modelling needs to be validated against simulations in any case (Avila et al. 2018; Fumagalli et al. 2022; Euclid Collaboration 2024a).

In this section we test some simple theoretical models for the covariance matrix of the angular two-point correlation function and power spectrum. The validation of semi-analytic models for the covariance of the 2D clustering is beyond the scope of this paper, while for the 3D two-point correlation function we rely on the results presented in the works by Fumagalli et al. (2022), Euclid Collaboration (2024a).

B.1. Covariance matrix of C_ℓ

Fig. B.1. Errors in the angular power spectrum, derived as the square root of diagonal elements of Covℓℓ′, in real (blue squares), redshift (with σz = 0, black circles) and in photo-z space, with σz = 0.02, (red triangles) and σz = 0.05 (green triangles), from the numerical covariance matrix measured from 1000 mocks. The lines show the results obtained with the linear power spectrum in the Limber approximation (solid black line), with the exact integration given by Eq. (43), which also includes the modelling of the RSDs (dashed red line) and with the damped linear power spectrum in the Limber approximation (solid red and green lines).

Concerning angular clustering, several works in the literature considered internal subsampling error estimates, such as bootstrap or jackknife (e.g. Paech et al. 2017; Balaguera-Antolínez et al. 2018). These methods are quite accurate for 2D clustering, because the radial projection alleviates the tensions that emerge in the 3D correlation (Norberg et al. 2009), depending on the scales involved, ensuring a good agreement with theoretical models (Crocce et al. 2011).

In the linear regime we can assume a Gaussian distribution for the a_ℓms. As a consequence of the Gaussian theory, covariance matrices are limited to diagonal terms (Dodelson 2003). This analytical estimation has been used in several works (Blake et al. 2007; Thomas et al. 2011; Balaguera-Antolínez et al. 2018) and can be written as

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\ell {\ell }^{\prime }}=\frac{2}{(2\ell +1)\mathrm{\Delta }\ell f_{\mathrm{sky}}}{(C_{\ell }+\frac{{\mathrm{\Omega }}_{\mathrm{sky}}}{N_{\mathrm{h}}})}^2{\delta }_{\ell {\ell }^{\prime }}^K,\end{array}$$(B.1)

where the Kronecker delta, δ_ℓℓ′^K, ensures the diagonality of the covariance matrix. It consists of an average over 2ℓ+1 different realisations, one for each m in a given multipole, with the cosmic variance term, C_ℓ, and the shot noise term, Ω_sky/N_h, that depends on the number of haloes in the selected redshift bin. The partial sky coverage enters in Eq. (B.1) as a boost factor 1/f_sky, though the covariance matrix remains diagonal. This is accurate enough because, as we have already discussed, the off-diagonal elements of the correlation matrix become negligible when averaging in the ℓ bands and consequently dividing for their size Δℓ. Equation (B.1) depends on the theoretical power spectrum C_ℓ, which includes the modelling of observational effects, such as photo-z errors and RSDs, according to Eqs. (37) and (43), respectively. In Fig. B.1 we show the measured errors, which are obtained as the square root of the diagonal elements of the numerical covariance matrix, ${\widehat{C}}_{\mathit{abij}}$. These are compared with the theory, finding a generally good agreement, within 10%. By looking at Fig. 7, the fact that the numerical errors slightly underestimate the model predictions can be explained as follows. The partial sky coverage of the PLCs determines the emergence of a non-diagonal correlation, spread over an interval of multipoles given approximately by ±1/f_sky, as already found by Crocce et al. (2011). Since for our mocks f_sky = 0.25, the expected correlation width is Δℓ = 8 and the measured non-diagonal terms rapidly drop to a noise consistent with zero. As discussed in Crocce et al. (2011), the integral of the elements of the covariance matrix arranged along a row gives us

$$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}}{\ell }^{\prime }{\mathrm{Cov}}_{\ell {\ell }^{\prime }}\approx \frac{2}{(2\ell +1)\mathrm{\Delta }\ell f_{\mathrm{sky}}}.\end{array}$$(B.2)

In full-sky surveys, all the errors are summarised in the diagonal of the covariance matrix. However, the presence of boundaries causes a leakage of the diagonal errors to other ℓ modes, so that their values become lower than the simple f_sky boost given by Eq. (B.1). As we have seen in Fig. A.1, thanks to the large sky coverage offered by the PINOCCHIO PLCs, almost all the correlation is concentrated within the band Δℓ. Therefore, the diagonal model in Eq. (B.1) is accurate enough to describe the numerical covariance. For smaller and highly fragmented regions of the sky, the approximation in Eq. (B.1) does not hold. This is clearly visible in Fig. B.2, where we show the result, after repeating the analysis but considering the 1004 deg² area subsample. Here, the diagonal approximation (black lines) overestimates the numerical errors in all redshift bins, with a more dramatic effect in the analysis of the reduced area. This issue can be resolved by relaxing the assumption of diagonality of the covariance matrix. We can think of the Kronecker delta as the full-sky limit of the mixing matrix, R_ℓℓ′. Thus, we can make the mixing matrix block-diagonal and compute the covariance as (Loureiro et al. 2019)

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\ell {\ell }^{\prime }}=\frac{2R_{\ell {\ell }^{\prime }}}{(2\ell +1)f_{\mathrm{sky}}^2}{(C_{\ell }+\frac{{\mathrm{\Omega }}_{\mathrm{sky}}}{N_{\mathrm{h}}})}^2,\end{array}$$(B.3)

with the second f_sky in the denominator coming up as the normalisation factor of the mixing matrix.

Fig. B.2. Comparison between diagonal and non-diagonal errors on the real-space angular power spectrum, Cℓ, in different redshift bins. Grey squares and magenta pentagons are measured from the set of mocks, with the full PLCs sky coverage and with its 1000 deg2 subregions, respectively. The black lines are computed through the Limber approximation, with the diagonal $1/\sqrt{f_{\mathrm{sky}}}$ boost factor. The purple lines account for the mode-mode coupling induced by the mixing matrix, Rℓ, ℓ′. Solid lines refer to the full PLCs sky coverage, while dashed lines refer to 1000 deg2 subregions.

Furthermore, as we have already seen in Sect. 5.1, the role of the RSDs is only relevant at the largest angular scales and increases progressively with z. However, for 0.6 < z < 0.8 and 0.8 < z < 1.0, we find a percentage difference in errors between real and redshift space greater than 10% only at ℓ < 10. To account also for non-linearities at the BAO scale, we use the best-fit Σ_NL values of Eq. (9), derived by fixing the cosmological parameters at Planck Collaboration XVI (2014), the ones of the PINOCCHIO simulations. However, their role in the C_ℓ error estimate is completely negligible, so they are not considered in the plot. As noted by Crocce et al. (2011), the damping of the baryon acoustic features has no effect on the error predictions on the large angular scales.

Analogously to what we found in Sect. 5.1, photo-z uncertainties also determine a suppression in the error estimation, which is modelled by convolving the angular power spectrum model, C_ℓ, with the Gaussian distribution in Eq. (1).

B.2. Covariance matrix of w(θ)

The angular correlation function can be derived analytically from C_ℓ, by summation over an infinite set of multipoles, sampling progressively smaller angular scales,

$$\begin{array}{c}\hfill w(\theta )={\displaystyle \sum _{\ell }}\left(\frac{2\ell +1}{4\pi }\right)L_{\ell }(cos\theta )C_{\ell },\end{array}$$(B.4)

without the shot noise term, since it only contributes to the zero separation limit (Chan et al. 2018), and also propagating systematic effects, such as RSDs and partial sky coverage. The covariance matrix of the angular correlation function is related to the C_ℓ covariance through (Cabré et al. 2007; Crocce et al. 2011; Ross et al. 2011)

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{{\theta }_i{\theta }_j}={\displaystyle \sum _{\ell ,{\ell }^{\prime }}}{\left(\frac{2\ell +1}{4\pi }\right)}^2{L}_{\ell }({\mu }_i)L_{{\ell }^{\prime }}({\mu }_j){\mathrm{Cov}}_{\ell {\ell }^{\prime }},\end{array}$$(B.5)

where θ_i and θ_j indicate the angular bins, while μ = cos θ. This yields

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{{\theta }_i{\theta }_j}=\frac{2}{f_{\mathrm{sky}}}{\displaystyle \sum _{\ell }}\frac{2\ell +1}{{(4\pi )}^2}{L}_{\ell }({\mu }_i)L_{\ell }({\mu }_j){(C_{\ell }+\frac{1}{\overline{n}})}^2,\end{array}$$(B.6)

where we set the shot noise to ${\mathrm{\Omega }}_{\mathrm{sky}}/N_{\mathrm{h}}=1/\overline{n}$ to simplify the notation.

As for the power spectrum case, the w(θ) covariance matrix is modified by the fact that measurements are computed over finite angular width bins, which reduces the non-diagonal correlation. A correct estimate requires the bin-averaged Legendre polynomials, defined as (Salazar-Albornoz et al. 2014, 2017; Chan et al. 2018)

$$\begin{array}{cc}& {\widehat{L}}_{\ell }\equiv \frac{{\displaystyle {\int }_{{\theta }_{-}}^{{\theta }_{+}}}{L}_{\ell }(cos\theta )sin\theta \mathrm{d}\theta }{{\displaystyle {\int }_{{\theta }_{-}}^{{\theta }_{+}}}sin\theta \mathrm{d}\theta }=\hfill \\ \hfill & =\frac{L_{\ell +1}\left({\mu }_{+}\right)-L_{\ell +1}\left({\mu }_{-}\right)-L_{\ell -1}\left({\mu }_{+}\right)+L_{\ell -1}\left({\mu }_{-}\right)}{(2\ell +1)({\mu }_{-}-{\mu }_{+})},\hfill \end{array}$$(B.7)

where θ₊ = θ + Δθ/2 and θ₋ = θ − Δθ/2 represent the upper and the lower limits of the bin, while μ₊ and μ₋ are the corresponding cosines. This leads to (Salazar-Albornoz et al. 2014)

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{{\theta }_i{\theta }_j}=\frac{2}{f_{\mathrm{sky}}}{\displaystyle \sum _{\ell }}\frac{2\ell +1}{{(4\pi )}^2}{\widehat{L}}_{\ell }({\mu }_i){\widehat{L}}_{\ell }({\mu }_j){(C_{\ell }+\frac{1}{\overline{n}})}^2.\end{array}$$(B.8)

Notably, even the covariance matrix for partial sky coverage follows the inverse of f_sky. However, this approximation might break when measuring w(θ) in highly fragmented survey areas, where the number of pairs as a function of separation does not simply scale with the effective area, Ω_sky, and thus it cannot be reproduced by the shot noise term, Ω_sky/N_h (Chan et al. 2018).

Moreover, since the shot noise does not depend on the multipole ℓ, we can make the binomial square explicit, and take $1/{\overline{n}}^2$ out from the summation. In particular, we can exploit the relation

$$\begin{array}{c}\hfill {\displaystyle \sum _{\ell }}(2\ell +1)\widehat{L}(\mu )\widehat{L}({\mu }^{\prime })=\frac{2}{{\mu }_{-}-{\mu }_{+}}{\delta }_K^{\theta {\theta }^{\prime }},\end{array}$$(B.9)

obtaining (Chan et al. 2018)

$$\begin{array}{cc}\hfill {\mathrm{Cov}}_{{\theta }_i{\theta }_j}& =\frac{{\delta }_K^{{\theta }_i{\theta }_j}}{4{\pi }^2{f}_{\mathrm{sky}}{\overline{n}}^2({\mu }_{-}-{\mu }_{+})}+\hfill \\ \hfill & +{\displaystyle \sum _{\ell }}\frac{2\ell +1}{8{\pi }^2{f}_{\mathrm{sky}}}{\widehat{L}}_{\ell }({\mu }_i){\widehat{L}}_{\ell }\left({\mu }_j\right)(C_{\ell }^2+\frac{2C_{\ell }}{\overline{n}}),\hfill \end{array}$$(B.10)

where the Kronecker δ_K^θ_iθ_j implies that the first term only acts on the diagonal of the covariance matrix.

Actually, despite its mathematical equivalence, Eq. (B.10) is more convenient than Eq. (B.8). The convergence to the final result is ensured by taking into account the contribution of infinite multipoles. Dealing with high-mass tracers, the shot noise term becomes dominant over the signal for ℓ ≳ 150, depending on the redshift bin. From this scale, its contribution is dominant, while scale-dependent terms containing C_ℓ behave as secondary corrections. As a consequence, the analytical formulation of the first term of Eq. (B.10) allows the summation to be truncated to a reasonable ℓ_max, whereas with Eq. (B.8) the absence of higher multipoles leads to an underestimation of the error at small θ.

Fig. B.3. Errors in the angular correlation function, derived as the square root of diagonal elements of Covθiθj, in real (blue squares), redshift (with σz = 0, black circles) and photo-z space, with σz = 0.02, (red triangles) and σz = 0.05 (green triangles), from the numerical covariance matrix measured from 1000 mocks. The lines show the results obtained with the linear power spectrum in the Limber approximation (solid black line), with the exact integration given by Eq. (43), which also includes the modelling of the RSDs (dashed red line) and with the damped linear power spectrum in the Limber approximation (solid red and green lines).

In Fig. B.3 we can see the errors in the angular correlation function, derived as the square root of the diagonal elements of Cov_θiθj, in real, in redshift and in photo-z space. As for Cov_ℓℓ′, we notice that there is a good agreement between the numerical estimate and the theoretical model, with the former overestimating the latter by a maximum of 10%, in the interval 0.8 < z < 1.0. On the other hand, due to the range of angular scales involved, we see no differences between the Limber approximation and the C_ℓ model including the RSDs. This is valid also for the measurements, since the ratio between the real-space and redshift-space errors is always greater than 0.95, meaning that there is no substantial difference between the two covariance matrices. In addition, we find that the theoretical red and green lines of Fig. B.3, for σ_z, 0 = 0.02 and σ_z, 0 = 0.05 respectively, get progressively closer as the redshift increases. This is due to the fact that, at z > 0.6, according to Eqs. (B.1) and hence (B.5), the errors becomes progressively more dominated by the shot noise and variations due to the C_ℓ term, given by photo-z damping and RSDs, are no longer appreciable.

Fig. B.4. Comparison between measurements and model in the full set of the w(θ) correlation matrices. The rows correspond to the different redshift bins, while different columns refer to real, redshift, with σz, 0 = 0, and photo−z space, with σz, 0 = 0.02 and σz, 0 = 0.05, from left to right. The upper triangles show the numerical estimate made on 1000 mocks. The lower triangles show the normalised covariance model.

In Fig. B.4 we can see the full set of correlation matrices. The rows correspond to the different redshift bins and indicate that the correlation decreases with z, while the columns refer to the analysis in real, redshift or photo−z space, with angular bins being less correlated in the presence of photo-z errors. The upper triangles show the numerical estimate made on 1000 PLCs. The lower triangles show the bin-averaged model presented in Eq. (B.10), where C_ℓs are computed with the Limber approximation.

Appendix C: Investigating further observational systematics

Fig. C.1. Mass proxy in cosmological analysis. Left panel: Mass-observable distribution of galaxy clusters in our mock catalogue. The grey shaded region represents the mass selection of our sample, excluding clusters with Mvir < 1014 M⊙ h−1. Right panel: Marginalised posterior distributions in the Ωm − σ8 plane, with 68% and 95% confidence intervals, using a binwidth of Δz = 0.4, in redshift space. The red contour shows the result for the two-point correlation function monopole, while the blue contour accounts also for the effect of the mass-observable scaling relation. The dashed lines indicate the Ωm and σ8 values of the PINOCCHIO simulation, while the dot-dashed line shows the predicted S8 degeneracy curve.

Using observations at different wavelengths, it is possible to calibrate parametric scaling relations between the mass of galaxy clusters and a generic observable, O^obs, which can be used as mass proxy. The so-called mass-observable scaling relation represents a key link between the theoretical mass function and bias – with their cosmological dependence – and the distribution of clusters in the space of survey observables. However, this introduces a degeneracy with ΛCDM parameters, which might significantly degrade the constraining power of cluster clustering and other cluster cosmological probes.

Therefore, accurate constraints require us to properly take into account the dependence of the mass on these additional parameters. To quantify possible systematic uncertainties due to a non-accurate mass estimation, we model the scaling relation as

$$\begin{array}{c}\hfill log\frac{O^{\mathrm{obs}}}{O_{\mathrm{piv}}^{\mathrm{obs}}}=A+Blog\frac{M_{\mathrm{vir}}}{M_{\mathrm{piv}}}+Clog\frac{1+z}{1+z_{\mathrm{piv}}},\end{array}$$(C.1)

where we fix $O_{\mathrm{piv}}^{\mathrm{obs}}=10$, M_piv = 10¹⁴ M_⊙ h⁻¹ and z_piv = 0.60. Here, A, B and C represent the normalisation, the slope and the redshift evolution of the scaling relation, respectively. In practice, starting from the PINOCCHIO catalogue, we extract a value for O^obs according to a log-normal distribution, with a mean derived from Eq. (C.1) and a standard deviation equal to the intrinsic scatter, σ_intr:

$$\begin{array}{c}\hfill P(log{O}^{\mathrm{obs}}|M_{\mathrm{vir}},z)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{intr}}}\exp [-\frac{x^2(O^{\mathrm{obs}},M_{\mathrm{vir}},z)}{2{\sigma }_{\mathrm{intr}}^2}].\end{array}$$(C.2)