© The Authors 2026

1. Introduction

Galaxy clusters, as tracers of the large-scale structure geometry and evolution, are well-known cosmological probes and are sensitive to the properties of the initial density field of the Universe, to the nature of dark matter and dark energy, and to the laws of gravity on large scales (e.g. Allen et al. 2011; Kravtsov & Borgani 2012). So far, the primary method for extracting cosmological information from cluster catalogues has been through analysis of their mass- and redshift-dependent number counts (e.g. Borgani et al. 2001; Henry 2004; Vikhlinin et al. 2009; Mantz et al. 2015; Planck Collaboration XXIV 2016; Costanzi et al. 2021; Lesci et al. 2022; Bocquet et al. 2024; Ghirardini et al. 2024; Lesci et al. 2025), which allows constraints on the matter density parameter (Ω_m) and the amplitude of density fluctuations (σ₈) in the Universe to be derived.

Beyond traditional methods such as cluster counts, cluster clustering stands out as a powerful technique for cosmology (Borgani et al. 1999; Estrada et al. 2009; Marulli et al. 2018). By describing the spatial distribution of galaxy clusters, clustering captures features that are missed by simple cluster abundance, allowing us to deepen our understanding of the Universe’s large-scale structure. Despite the limitations of current statistics, which prevent cluster clustering from being a competitive stand-alone probe, this approach has proven to be a valuable tool in breaking degeneracies in cosmological and mass-observable parameters when combined with complementary observables such as cluster abundance (Schuecker et al. 2003; Majumdar & Mohr 2004; Mana et al. 2013; Lesci et al. 2022; Fumagalli et al. 2024).

The constraining power of number counts and clustering is expected to increase significantly with the upcoming datasets from Stage IV wide-surveys, such as eROSITA¹ (Predehl 2014), Euclid² (Laureijs et al. 2011), and Vera C. Rubin Observatory Legacy Survey of Space and Time³ (LSST; Abate et al. 2012). In particular, the successful launch of the ESA space mission Euclid (Euclid Collaboration: Scaramella et al. 2022; Euclid Collaboration: Mellier et al. 2025) in July 2023 marks the beginning of its ambitious goal of mapping 14 000 deg² of the extragalactic sky and collecting an unprecedented amount of cosmological data. For galaxy clusters, Euclid is expected to yield a sample of about 10⁵ objects out to redshift z = 2 (Sartoris et al. 2016), using photometric and spectroscopic data and through gravitational lensing. Handling such a large amount of data requires a thorough understanding of the systematic uncertainties, which will be the primary limitation in future cosmological analyses. Even with substantial progress made over the past decade in characterising these systematic effects, attaining the sub-percent accuracy required by Euclid remains a major challenge.

Accurate cosmological parameter constraints rely on effectively characterising uncertainties through covariance matrices. Direct calculation of these matrices requires computationally expensive simulations, and while approximate methods can help reduce the computational costs (Pope & Szapudi 2008; Monaco 2016), they can also yield noisy results unless a large number of mock realisations (on the order of 10³– 10⁴) are generated. Alternatively, covariance can be estimated using data-driven methods such as bootstrap or jackknife, but the results are noisy and can overestimate the covariance, particularly for two-point statistics (Norberg et al. 2009). Analytical methods offer a fast way to compute covariance matrices without heavy computational demands, producing noise-free, cosmology-dependent results (e.g. Scoccimarro et al. 1999; Meiksin & White 1999; Hu & Kravtsov 2003; Takada & Hu 2013). However, accurately describing all contributions remains challenging. This complexity necessitates validation against simulations to determine which contributions are significant at the required statistical precision, sometimes leading to the need for calibration of parameters that cannot be derived analytically (Xu et al. 2012; O’Connell et al. 2016; Fumagalli et al. 2022).

Fumagalli et al. (2021, hereafter F21) and Euclid Collaboration: Fumagalli et al. (2024, hereafter EC:F24) carried out the first steps for the characterisation of theoretical systematics in the cosmological analysis of cluster number counts (see also Payerne et al. 2023, 2024) and cluster clustering by validating (semi-)analytical covariance matrix models for the analysis of the upcoming Euclid catalogues. This work extends that process by assessing the impact of observational and modelling systematic effects, such as redshift-space distortions in the linear regime (RSDs; Kaiser 1987) and photometric redshift (photo-z) errors. We also investigate the possible cross-correlations between the two probes to develop a robust combined likelihood analysis. Another important contribution to the systematics of cluster cosmology arises from projection effects, which induce a correlated scatter between cluster observables and the surrounding large-scale environment. Such effects can bias cluster measurements, modifying the inferred cluster abundance and enhancing the clustering signal (e.g. Costanzi et al. 2019; Sunayama et al. 2020). Since the modelling of projection effects requires a dedicated and detailed treatment, we leave their investigation to future work and neglect them in the present analysis.

This paper is structured as follows. In Sect. 2 we present the theoretical framework to model the two statistical probes, that is, cluster counts and the two-point correlation function (2PCF), along with the corresponding covariance models as well as the likelihood models adopted to infer cosmological and mass-observable relation parameters. In Sect. 3 we describe the simulations analysed in this work and the measurement of the statistics and their numerical covariance matrices. In Sect. 4 we investigate the presence of cross-correlations (Sect. 4.1) and examine the probe combination (Sect. 4.2) as well as the impact of cosmology-dependent covariances (Sect. 4.3). Additionally, we provide a validation of cluster clustering (Sect. 4.4) and number counts (Sect. 4.5) modelling. Lastly, in Sect. 5 we discuss our results and summarise our main conclusions.

2. Theoretical background

In this section we introduce the theoretical formalism to describe number counts and cluster clustering, along with their covariance matrix models, originally presented in F21 and EC:F24, respectively. As demonstrated in Fumagalli et al. (2024), number counts and clustering alone do not provide sufficient constraints on scaling relations. As Euclid has been designed to deliver high quality space-based weak lensing measures, we also incorporate the weak lensing mean mass to further refine the constraints on the scaling relation parameters. We also describe the likelihood function adopted for the cosmological inference and the summary statistics used to compare the results. Throughout the paper, quantities labelled with ‘tr’ or ‘ob’ denote true intrinsic or observed quantities, respectively, while P(A|B) denotes the conditional probability of A given B. Note that we use the Dark Energy survey (DES; Abbott et al. 2005) as a reference to model certain quantities. While the values adopted for the Euclid analysis may differ, DES offers realistic benchmarks for this purpose.

2.1. Number counts

We describe the expected number counts of galaxy clusters in the a-th redshift bin and i-th richness bin as follows

$$\begin{array}{c}\hfill {⟨N⟩}_{\mathit{ai}}={\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}^{\mathrm{tr}}{\mathrm{\Omega }}_{\mathrm{sky}}\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}(z^{\mathrm{tr}}|\mathrm{\Delta }{z}_a^{\mathrm{ob}}){⟨n(z^{\mathrm{tr}})⟩}_i.\end{array}$$(1)

Here, Ω_sky is the survey area in steradians, assuming a homogeneous data quality over the sky area, and the observed comoving volume element per unit redshift and solid angle is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}(z^{\mathrm{tr}}|\mathrm{\Delta }{z}_a^{\mathrm{ob}})={\displaystyle {\int }_{\mathrm{\Delta }{z}_a^{\mathrm{ob}}}\mathrm{d}{z}^{\mathrm{ob}}}\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}(z^{\mathrm{tr}})P(z^{\mathrm{ob}}|z^{\mathrm{tr}},\mathrm{\Delta }{\lambda }_i^{\mathrm{ob}}),\end{array}$$(2)

where the term P(z^ob | z^tr, Δλ_i^ob) is the probability distribution that assigns an observed redshift z^ob to a cluster, given the true redshift and the observed richness bin Δλ_i^ob. This term incorporates the effect of photometric uncertainty that we assume to be modelled as a Gaussian distribution with a mean equal to z^tr and a scatter

$$\begin{array}{c}\hfill {\sigma }_z(z^{\mathrm{tr}})={\sigma }_{z0}(1+z^{\mathrm{tr}}),\end{array}$$(3)

where σ_z0 is the typical photometric redshift uncertainty of the survey associated with the cluster detection error.

The expected number density of haloes in the i-th richness bin is given by

$$\begin{array}{c}\hfill {⟨n(z^{\mathrm{tr}})⟩}_i={\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M,z^{\mathrm{tr}}){\displaystyle {\int }_{\mathrm{\Delta }{\lambda }_i^{\mathrm{ob}}}}\mathrm{d}{\lambda }^{\mathrm{ob}}P({\lambda }^{\mathrm{ob}}|M,z^{\mathrm{tr}}),\end{array}$$(4)

where dn/dM is the halo mass function and P(λ^ob | M, z^tr) denotes the observed richness-mass relation

$$\begin{array}{c}\hfill P({\lambda }^{\mathrm{ob}}|M,z^{\mathrm{tr}})={\displaystyle {\int }_0^{\infty }}\mathrm{d}{\lambda }^{\mathrm{tr}}P({\lambda }^{\mathrm{ob}}|{\lambda }^{\mathrm{tr}},z^{\mathrm{tr}})P({\lambda }^{\mathrm{tr}}|M,z^{\mathrm{tr}}).\end{array}$$(5)

The term P(λ^tr | M, z^tr) is the intrinsic richness-mass relation. In this work, we assumed a log-normal distribution with mean and scatter given by

$$\begin{array}{cc}\hfill ⟨ln{\lambda }^{\mathrm{tr}}|M_{\mathrm{vir}}⟩& =ln{A}_{\lambda }+B_{\lambda }ln\left(\frac{M_{\mathrm{vir}}}{M_{\mathrm{p}}}\right)+C_{\lambda }ln\left(\frac{1+z}{1+z_{\mathrm{p}}}\right),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\sigma }_{ln{\lambda }^{\mathrm{tr}}|M_{\mathrm{vir}}}^2& =D_{\lambda }^2+\frac{⟨{\lambda }^{\mathrm{tr}}|M_{\mathrm{vir}}⟩-1}{{⟨{\lambda }^{\mathrm{tr}}|M_{\mathrm{vir}}⟩}^2},\hfill \end{array}$$(7)

where M_p = 3 × 10¹⁴ h⁻¹ M_⊙ and z_p = 0.45 (Costanzi et al. 2021) are the pivot mass and redshift parameters, while A_λ, B_λ, C_λ, and D_λ are the mass-observable relation parameters. The variance consists of intrinsic scatter, D_λ, plus a Poisson contribution.

The second term in Eq. (5), P(λ^ob|λ^tr, z^tr), represents the observational scatter in the richness-mass relation, caused by photometric noise, background subtraction uncertainties, and projection or percolation effects. We assumed a normal distribution with scatter given by (Costanzi et al. 2021)

$$\begin{array}{c}\hfill {\sigma }_{{\lambda }^{\mathrm{ob}}}(\lambda ,z)=(0.9+0.1z){\lambda }^{0.4}.\end{array}$$(8)

For the number count covariance, we followed the model described in F21, adapted to account for richness-mass relation and selection functions. We describe the number count covariance as the sum of shot-noise and sample variance (Hu & Kravtsov 2003),

$$\begin{array}{c}\hfill C_{\mathit{abij}}={\langle N\rangle }_{\mathit{ai}}{\delta }_{\mathit{ab}}{\delta }_{\mathit{ij}}+{\langle Nb\rangle }_{\mathit{ai}}{\langle Nb\rangle }_{\mathit{bj}}{S}_{\mathit{ab}},\end{array}$$(9)

where ⟨N⟩_ai is the prediction for the number counts in the a-th redshift bin and i-th richness bin (Eq. 1), and ⟨Nb⟩_ai is the prediction for the number counts times the halo bias, computed in an analogous way. The term S_ab is the covariance of the linear density field between two redshift bins given by

$$\begin{array}{c}\hfill S_{\mathit{ab}}={\displaystyle \int \frac{{\mathrm{d}}^{3}k}{{(2\pi )}^3}}\sqrt{P_{\mathrm{m}}(k,{\overline{z}}_a)P_{\mathrm{m}}(k,{\overline{z}}_b)}{W}_a(\mathbf{k})W_b(\mathbf{k}),\end{array}$$(10)

where $P_{\mathrm{m}}(k,{\overline{z}}_a)$ is the linear matter power spectrum evaluated at the centre of the a-th redshift bin, and W_a(k) is the window function of the redshift bin. The window function for a redshift slice of a light cone is given in Costanzi et al. (2019) and takes the form

$$\begin{array}{c}\hfill W_a(\mathbf{k})=\frac{4\pi }{V_a}{\displaystyle {\int }_{\mathrm{\Delta }{z}_a}}\mathrm{d}z{\mathrm{\Omega }}_{\mathrm{sky}}\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}{\displaystyle \sum _{\ell =0}^{\infty }}{\displaystyle \sum _{m=-\ell }^{\ell }}{\mathrm{i}}^{\ell }{j}_{\ell }[kr(z)]Y_{\ell m}(\widehat{\mathbf{k}})K_{\ell },\end{array}$$(11)

where V_a is the volume of the redshift slice, j_ℓ[k r(z)] are the spherical Bessel functions, $Y_{\ell m}(\widehat{\mathbf{k}})$ are the spherical harmonics, $\widehat{\mathbf{k}}$ is the angular part of the wave-vector, and K_ℓ are the coefficients of the harmonic expansion.

2.2. Two-point correlation function

In this work we describe cluster clustering through the 3D 2PCF monopole. We considered the auto- and cross-correlation function between different richness samples, each one computed in redshift and radial separation bins. We modelled the 2PCF in the a-th redshift bin, s-th radial bin, and i-th and j-th richness bins as

$$\begin{array}{c}\hfill {⟨\xi ⟩}_{\mathit{asij}}={\displaystyle \int \frac{\mathrm{d}k{k}^2}{2{\pi }^2}}{⟨P_{\mathrm{h}}(k)⟩}_{\mathit{aij}}{W}_s(k),\end{array}$$(12)

where W_s(k) is the spherical shell window function given by

$$\begin{array}{c}\hfill W_s(k)={\displaystyle \int \frac{{\mathrm{d}}^{3}r}{V_s}}{j}_0(kr)=\frac{r_{s,+}^3{W}_{\mathrm{th}}(k{r}_{s,+})-r_{s,-}^3{W}_{\mathrm{th}}(k{r}_{s,-})}{r_{s,+}^3-r_{s,-}^3},\end{array}$$(13)

where W_th(kr) is the top-hat window function and V_s is the volume of the s-th spherical shell, delimited by r_s, − and r_s, +.

The term ⟨P_h(k)⟩_aij is the halo power spectrum integrated in the a-th redshift bin and in the i-th and j-th richness bins. It can be expressed as the combination of the auto-power spectra in each richness bin as

$$\begin{array}{c}\hfill {⟨P_{\mathrm{h}}(k)⟩}_{\mathit{aij}}=\sqrt{{⟨P_{\mathrm{h}}(k)⟩}_{\mathit{ai}}{⟨P_{\mathrm{h}}(k)⟩}_{\mathit{aj}}},\end{array}$$(14)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {⟨P_{\mathrm{h}}(k)⟩}_{\mathit{ai}}=\frac{1}{{⟨N⟩}_{\mathit{ai}}}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}^{\mathrm{tr}}& {\mathrm{\Omega }}_{\mathrm{sky}}\frac{\mathrm{d}V}{\mathrm{d}\mathrm{\Omega }\mathrm{d}z}(z^{\mathrm{tr}}|\mathrm{\Delta }{z}_a^{\mathrm{ob}})\hfill \\ \hfill & \times {⟨n(z^{\mathrm{tr}})⟩}_i{⟨b_{\mathrm{eff}}(z^{\mathrm{tr}})⟩}_i{P}_{\mathrm{m}}(k,z^{\mathrm{tr}}).\hfill \end{array}\end{array}$$(15)

In the above equation, ⟨b_eff(z^tr)⟩_i is the effective bias for all clusters having an observed richness in the i-th richness bin, and it is computed by weighting the halo mass function in Eq. (4) by the halo bias.

To take into account the RSD effect and the uncertainty in the photometric redshift measurements, the matter power spectrum in Eq. (15) had to be modified as (Marulli et al. 2012; Sereno et al. 2015)

$$\begin{array}{c}\hfill P_{\mathrm{m}}^{\mathrm{obs}}(k)=P\prime (k)+\beta P″(k)+{\beta }^{2}P″\prime (k),\end{array}$$(16)

with

$$\begin{array}{cc}\hfill P\prime (k)& =P_{\mathrm{m}}(k)\frac{\sqrt{\pi }}{2k\sigma }\mathrm{erf}(k\sigma ),\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill P″(k)& =P_{\mathrm{m}}(k)\frac{\beta }{{(k\sigma )}^3}[\frac{\sqrt{\pi }}{2}\mathrm{erf}(k\sigma )-k\sigma exp(-k^2{\sigma }^2)],\hfill \end{array}$$(18)

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

where β = f_g/b_eff⁴ with f_g ≈ [Ω_m(z)]^0.55 and σ depends on the photo-z uncertainty σ_z as

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

Equations (17)–(19) exclusively incorporate the linear part of the RSDs, commonly referred to as the Kaiser effect (Kaiser 1987). The so-called finger-of-god non-linear distortion introduces perturbations to the power spectrum similar in shape to those caused by photometric redshift uncertainties. However, its impact is significantly smaller compared to that of photo-z uncertainties and can therefore be safely neglected (Veropalumbo et al. 2014; Sereno et al. 2015; Romanello et al. 2025). The Kaiser effect produces a scale-independent boost of the 2PCF, while the effect of the photo-z uncertainties is to decrease the correlation on small scales and increase it on larger scales, smearing the baryon acoustic oscillation (BAO) peak.

The BAO scales are also subject to broadening and a shift of the peak due to non-linear damping, which can be corrected through the infrared resummation (IR; Senatore & Zaldarriaga 2015; Baldauf et al. 2015). At lowest order, the matter power spectrum is corrected as

$$\begin{array}{c}\hfill P_{\mathrm{m}}(k,z)\simeq P_{\mathrm{nw}}(k,z)+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2(z)}{P}_{\mathrm{w}}(k,z),\end{array}$$(21)

where P_w and P_nw are the wiggled and smooth parts of the linear power spectrum, respectively, and

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}^2(z)={\displaystyle {\int }_0^{k_s}}\frac{\mathrm{d}q}{6{\pi }^2}{P}_{\mathrm{nw}}(q,z)[1-j_0\left(q{r}_{\mathrm{BAO}}\right)+2j_2\left(q{r}_{\mathrm{BAO}}\right)].\end{array}$$(22)

Lastly, when measuring the 2PCF from data, it is necessary to make assumptions about the underlying cosmology to convert redshifts into distances. To take this assumption into account, a parameter was introduced to model the geometric distortions (GDs) in the 2PCF shape (Marulli et al. 2012, 2016):

$$\begin{array}{c}\hfill \xi (r)⟶\xi (\alpha r),\end{array}$$(23)

where r is the radial separation and

$$\begin{array}{c}\hfill \alpha =\frac{D_{\mathrm{V}}}{r_{\mathrm{s}}}\frac{r_{\mathrm{s}}^{\mathrm{fid}}}{D_{\mathrm{V}}^{\mathrm{fid}}}.\end{array}$$(24)

Here, r_s is the position of the sound horizon at decoupling, D_V is the isotropic volume distance, and the label ‘fid’ indicates the quantities evaluated at the fiducial cosmology assumed in the measurement process.

We describe the clustering covariance with the semi-analytical model presented in EC:F24. The cross-covariance between different richness samples i, j, k, and l, in the a-th redshift bin and between the s-th and r-th separation bins is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill C_{\mathit{asr}}^{\mathit{ijkl}}& =\frac{1}{V_a}{\displaystyle \int \frac{\mathrm{d}k{k}^2}{2{\pi }^2}}[{⟨P_{\mathrm{h}}(k)⟩}_{\mathit{aik}}+\frac{{\delta }_{\mathit{ik}}}{{⟨n⟩}_{\mathit{ai}}}]\hfill \\ \hfill & \times [{⟨P_{\mathrm{h}}(k)⟩}_{\mathit{ajl}}+\frac{{\delta }_{\mathit{jl}}}{{⟨n⟩}_{\mathit{aj}}}]W_s(k)W_r(k)\hfill \\ \hfill & +\frac{1}{V_a}{\displaystyle \int \frac{\mathrm{d}k{k}^2}{2{\pi }^2}}{⟨P_{\mathrm{h}}(k)⟩}_{\mathit{aij}}\frac{{\delta }_{\mathit{ik}}}{{⟨n⟩}_{\mathit{ai}}}\frac{{\delta }_{\mathit{jl}}}{{⟨n⟩}_{\mathit{aj}}}\frac{W_r(k)}{V_s}{\delta }_{\mathit{sr}}\hfill \\ \hfill & +(k⟷l),\hfill \end{array}\end{array}$$(25)

where the term (k ↔ l) denotes that the preceding expression should be repeated with the indices k and l swapped. The first integral represents the standard Gaussian covariance, while the second one is the low-order non-Gaussian contribution.

According to EC:F24, the bias and shot-noise terms in the clustering covariance must be modified by the addition of three parameters {α_c,  β_c,  γ_c} such that

$$\begin{array}{cc}\hfill b& ⟶{\beta }_{\mathrm{c}}b,\hfill \\ \hfill \frac{1}{n}& ⟶\frac{1+{\alpha }_{\mathrm{c}}}{n}\text{in the Gaussian term},\hfill \\ \hfill \frac{1}{n}& ⟶\frac{1+{\gamma }_{\mathrm{c}}}{n}\text{in the non-Gaussian term}.\hfill \end{array}$$

These parameters can be fitted from a few simulations (of the order of 10²; Fumagalli et al. 2022) in each redshift and richness bin; the reference values are α_c = 0, β_c = 1, and γ_c = 0. As shown in Fumagalli et al. (2025, their figure 8), the cosmology dependence of such parameters is absent for β_c and γ_c, and it is almost negligible for α_c, thus allowing us to perform a single calibration at the fiducial cosmology of the simulations (Sect. 3.1).

2.3. Weak lensing mean mass

Similar to the number counts approach, the expected value of the mean cluster mass within the a-th redshift bin and i-th richness bin is given by

$$\begin{array}{c}\hfill {\overline{M}}_{\mathit{ai}}=\frac{⟨M^{\mathrm{tot}}(\mathrm{\Delta }{\lambda }_i,\mathrm{\Delta }{z}_a)⟩}{⟨N(\mathrm{\Delta }{\lambda }_i,\mathrm{\Delta }{z}_a)⟩},\end{array}$$(26)

where ⟨M^tot⟩ represents the total mass associated with clusters identified within a given redshift and richness range, expressed as

$$\begin{array}{c}\hfill ⟨M^{\mathrm{tot}}(\mathrm{\Delta }{z}_i,\mathrm{\Delta }{\lambda }_j)⟩={\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}^{\mathrm{tr}}{\mathrm{\Omega }}_{\mathrm{sky}}\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}(z^{\mathrm{tr}}|\mathrm{\Delta }{z}_a^{\mathrm{ob}}){⟨Mn(z^{\mathrm{tr}})⟩}_i,\end{array}$$(27)

where the term ⟨M n(z^tr)⟩_i is computed by weighting the halo mass function in Eq. (4) by the halo mass and accounts for the expected mass contribution from clusters at a given redshift and richness interval.

2.4. Likelihood function

We evaluated the impact of systematic effects on cosmological constraints by performing a Bayesian inference on Euclid-like cluster catalogues derived from 1000 simulated light cones sharing the same underlying cosmological model (see Sect. 3.1). We explored the posterior distribution by using the python wrapper for the nested sampling PolyChord (Handley et al. 2015).

For number counts and clustering, as well as for weak lensing log-masses, we adopted Gaussian likelihoods:

$$\begin{array}{c}\hfill \mathcal{L}(\mathbf{x}|\mathit{\mu },\mathsf{C})=\frac{exp\{-\frac{1}{2}{(\mathbf{x}-\mathit{\mu })}^T{C}^{-1}(\mathbf{x}-\mathit{\mu })\}}{\sqrt{2\pi det[\mathsf{C}]}},\end{array}$$(28)

where x and μ are, respectively, the measurement and the model prediction, and C is the covariance matrix. We assumed weak lensing masses to be uncorrelated with both counts and clustering based on the assumptions made in Costanzi et al. (2019) and Fumagalli et al. (2024). We also considered number counts and cluster clustering as independent probes, following Mana et al. (2013) and Fumagalli et al. (2024), and we further test the absence of cross-correlation through the use of simulations (see Sect. 4.1). Unless differently specified, for cluster counts and clustering, we used cosmology-dependent covariances computed through the analytical models of Eqs. (9) and (25), respectively.

Both the cosmological parameters and the parameters describing the mass-observable relation are determined by maximising the average log-likelihood computed across all the N_s simulated catalogues; this approach helps mitigate the cosmic variance effect characterising each Universe realisation by a factor of $\sqrt{N_{\mathrm{s}}}$. As a result, we are able to detect potential systematic effects in the analysis, identified by any residual shift of the posteriors with respect to the input cosmology.

To assess the difference in the final cosmological posteriors, we use two summary statistics: the figure of merit (FoM; Albrecht et al. 2006), which measures the relative amplitude of the posteriors compared to a reference case; and the index of inconsistency (IoI; Lin & Ishak 2017), which quantifies their shift from the fiducial cosmology. By running multiple realisations of the same chain, we find that the statistical noise associated with the convergence of likelihood maximisation produces an error of less than 4% on the FoM and 1% on the IoI.

3. Simulated data

In this section we describe the simulated light cones adopted in this work. We also outline the procedure for the measuring of numerical covariances.

3.1. Simulations

Accurate covariance matrix estimation necessitates extensive datasets, often exceeding 10³ catalogues (Taylor et al. 2013; Dodelson & Schneider 2013). Instead of computationally expensive N-body simulations, we employed the PINOCCHIO code (Monaco et al. 2002; Munari et al. 2017), which utilises third-order Lagrangian perturbation theory and ellipsoidal collapse to quickly generate dark matter halo catalogues. PINOCCHIO generates initial density fields on a regular grid, and uses Lagrangian perturbation theory to first compute particle collapse times, and then to identify haloes and trace their evolution to final positions. Large simulation boxes were used to construct past light cones, capturing haloes causally connected to an observer at present time. As shown in Fumagalli et al. (2025, their Sect. 4.2), on linear scales PINOCCHIO produces numerical covariances that accurately match those obtained from N-body simulations.

Our dataset consists of 1000 light cones covering 10 313 deg² and redshifts z = 0 to z = 2.5. These catalogues, generated from boxes of side 3870 h⁻¹ Mpc and 2160³ particles, contain haloes with M_vir ≥ 10^13.4 h⁻¹ M_⊙. The simulations are based on a flat ΛCDM cosmology consistent with Planck Collaboration XVI (2014): Ω_m = 0.30711, Ω_b = 0.048254, h = 0.6777, n_s = 0.96, ∑m_ν = 0.06 eV, A_s = 2.21 × 10⁻⁹, and σ₈ = 0.8288. Similarly to F21, we rescaled the halo masses to match, on average, a reference halo mass function model, while preserving shot-noise and sample variance fluctuations in each mock. In this work we assume the Euclid Collaboration: Castro et al. (2023) halo mass function. We describe the halo bias with the Euclid Collaboration: Castro et al. (2024) model, computed and calibrated for the assumed halo mass function. Both models have been calibrated to a percent level of accuracy, ensuring that any uncertainty in the formalism used for modelling these two quantities is negligible. The resulting catalogues contain around 10⁵ haloes with M_vir ≥ 10^13.5h⁻¹ M_⊙ each. Each object is labelled with a rescaled virial mass, a true undistorted redshift, and a redshift with RSDs, along with angular coordinates.

To mimic the observational properties of cluster surveys, we attached a true richness to each mass by assuming a scaling relation with mean and variance from Eqs. (6) and (7), with parameters A_λ = 50, B_λ = 1.0, C_λ = 0.01, and D_λ = 0.2, obtained from Costanzi et al. (2021)⁵. Such a richness is scattered to produce an observed richness by assuming a log-normal distribution with scatter from Eq. (8).

In addition, two observed redshifts are assigned assuming a normal distribution with scatter given by Eq. (3): for the typical photo-z uncertainty, we consider an optimistic (σ₀ = 0.005) and a more realistic (σ₀ = 0.01) estimate of the expected Euclid photo-z errors (Euclid Collaboration: Bhargava et al. 2026). We define the values based on the photo-z uncertainties of Stage-III surveys (see, for instance, Maturi et al. 2019; Abbott et al. 2020; DES Collaboration 2025), which serve as a benchmark that Stage-IV surveys aim to improve. Unless otherwise specified, we assume the conservative scenario as the default in the analysis.

Apart from introducing scatter around the true value, photo-z measurements can also bias cluster redshifts. However, since the Euclid requirement of Δz = 0.002 (Laureijs et al. 2011) results in a negligible impact on our analysis (see Appendix A), we did not apply any shift to our redshifts.

Lastly, it should be noted that while the catalogues effectively replicate observational conditions, this work does not account for the effects of any survey mask or inhomogeneous sky coverage, and we assume a simplified model of P(λ^ob|λ^tr, z^tr) that does not take into account the correlation with the LSS, since this function has not yet been calibrated for the Euclid data.

3.2. Measurements of the summary statistics

We measure the number counts assuming twenty redshift bins of width Δz = 0.1 in the range z = 0–2, and five log-spaced richness bins in the range λ = 20–400.

The same binning scheme is adopted for weak lensing masses. However, since our simulations do not include lensing information, we generate synthetic measurements of the mean weak lensing mass. This is done by predicting a true mean mass in each richness and redshift bin and perturbing it with a Gaussian distribution, assuming a diagonal constant error of 1%, which aligns with the precision target for Stage IV surveys⁶. To match the number of simulations, we generate 1000 sets of synthetic measurements by varying the seed of the Gaussian scatter. According to Costanzi et al. (2019), we model the cosmological dependence on Ω_m of the weak lensing mass estimates as

$$\begin{array}{c}\hfill {log}_{10}{\widehat{M}}_{\mathrm{WL}}({\mathrm{\Omega }}_{\mathrm{m}})={log}_{10}{\widehat{M}}_{\mathrm{WL}}{|}_{{\mathrm{\Omega }}_{\mathrm{m}}=0.3}+\frac{\mathrm{d}{log}_{10}{M}_{\mathrm{WL}}}{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}}({\mathrm{\Omega }}_{\mathrm{m}}-0.3),\end{array}$$(29)

where the slopes derived in each bin from fitting such an equation to the data are listed in Table 1 of Costanzi et al. (2019), and are used in our cosmological analysis to re-scale the ‘observed’ mean mass at each step of the MCMC.

We measured the 2PCF by using the Landy & Szalay (1993) estimator, implemented via the CosmoBolognaLib package (Marulli et al. 2016). To construct the random catalogue, we randomly extracted a subset of objects from each mock, shuffled their coordinates, and combined them into a single catalogue containing N_R = 20 N_D objects randomly distributed within the light cone volume. We assumed five redshift bins of width Δz = 0.4 and two richness bins defined as λ = [20, 40, 400]. For the standard analysis, we used 25 separation bins in the range r = 20–130 h⁻¹ Mpc.

4. Results

In this section we present the results of the covariance and joint analysis validation. Our analysis begins with an assessment of potential cross-correlations between the two probes, that is, number counts and 2PCF monopole, using numerical simulations. We then compare the cosmological posteriors from different probe combinations to evaluate their constraining power for a Euclid-like survey. Next, we examine the impact of photo-z uncertainties and redshift distortions on the clustering and number counts modelling, and we assess the impact on the cosmological constraints obtained through the likelihood analysis of clustering alone, and in combination with number counts. In both cases, we add the weak lensing mass likelihood to break the degeneracies between cosmological and scaling-relation parameters. For the likelihood analysis, we assume flat priors on both the cosmological parameters Ω_m and log₁₀ A_s, and on the four scaling-relation parameters of Eqs. (6) and (7). Specifically, we use Ω_m ∈ [0.2, 0.4], log₁₀ A_s ∈ [ − 9.0, −8.5], A_λ ∈ [20, 80], B_λ ∈ [0.5, 1.5], C_λ ∈ [ − 0.5, 0.5], and D_λ ∈ [0.01, 0.5]. All the other parameters are fixed to the fiducial values, while σ₈ is a derived parameter.

4.1. Cross-covariance

To properly combine cluster counts and cluster clustering in the likelihood analysis, we should account for the cross-correlation between the two probes. We compute the auto- and cross-covariance matrices for the two probes from simulations, as shown in the left panel of Fig. 1. The cross-correlation term is fully dominated by noise and consistent with zero. The result holds independently of redshift. As a further proof, we compute the number counts and clustering likelihoods for each light cone, and in the right panel of Fig. 1 we show their difference with respect to the log-likelihood averaged over the entire sample of simulations. We quantify the correlation by computing the Pearson correlation coefficient between the log-likelihood residuals. For this purpose, we assumed that the residuals are Gaussian distributed around zero, with a covariance given by

$$\begin{array}{c}\hfill C=\left(\begin{array}{cc}{\sigma }_{\mathrm{NC}}^2& \rho {\sigma }_{\mathrm{NC}}{\sigma }_{\mathrm{CL}}\\ \rho {\sigma }_{\mathrm{NC}}{\sigma }_{\mathrm{CL}}& {\sigma }_{\mathrm{CL}}^2\end{array}\right),\end{array}$$(30)

Fig. 1. Cross-covariance between number counts and clustering. For better visualisation, we used redshift bins of width Δz = 0.2 for counts and Δz = 0.5 and a single richness bin for the 2PCF. Left: Auto- and cross-correlation matrix of number counts and clustering, computed from 1000 mocks. Right: Log-likelihood residuals for number counts and clustering for each one of the 1000 light cones, with respect to the mean value assuming the fiducial model parameters.

where σ_NC and σ_CL are, respectively, the standard deviations computed from the number counts and clustering residuals, and ρ is the correlation coefficient. We computed the Pearson coefficient from the data, obtaining ρ = −0.011 ± 0.031, where the error was estimated using bootstrap. The absence of correlation in the distribution of points confirms the absence of correlation between the two observables. According to the result, number counts and clustering can be considered as two independent probes. The independence of number counts and clustering is due to the fact that the former is a measure of the average comoving number density of tracers, while the latter measures the fluctuations around this mean value. Since fluctuations around the mean can be arbitrarily varied within a sufficiently large volume of the Universe without altering the mean density of tracers, the two probes turn out to be independent of each other.

4.2. Probe combination

We explore different combinations of the three probes – number counts, 2PCF monopole, and weak lensing mean mass – to evaluate their constraining power for a Euclid-like survey. As shown in Fig. 2, our results align with those of Fumagalli et al. (2024), confirming the benefits of incorporating cluster clustering as a cosmological probe.

Fig. 2. Parameter posterior distributions with 68% and 95% confidence intervals from different combinations of probes: number counts and weak lensing mass (in blue), clustering and weak lensing mass (in orange), number counts and clustering (in green), and all the three probes (in red). Dotted grey lines are the input cosmology of the catalogues.

Examining the number counts and weak lensing mass posteriors (blue contours), we note a slight bias, though at a low significance (IoI = 0.24), relative to the fiducial cosmology. This could stem from prior-volume effects. Smaller biases are also visible in other two-probe combinations. However, the full-probe combination is correctly centred, demonstrating that integrating the three probes produces not only more precise, but also more accurate results. This combination yields the tightest constraints, with clustering and weak lensing mass providing the second strongest constraints (see Table 1). Therefore, in the subsequent analysis, we primarily focus on these two cases.

4.3. Cosmology dependence of the covariances

Here, we assess the impact of cosmology-dependent covariances on cosmological posteriors. As demonstrated in F21 and EC:F24, fixing the covariance matrix at an incorrect cosmology can alter the amplitude of the posteriors, whereas using a cosmology-dependent covariance matrix avoids this issue. Additionally, allowing the clustering covariance to vary freely in the likelihood analysis can provide extra cosmological information. As shown in the lower-left panel of Fig. 3, we confirm that despite the broader contours caused by uncertainties in the scaling relation parameters, the clustering covariance still improves the cosmological constraints from cluster clustering, leading to a 40% increase in the FoM. Also, the use of a fixed covariance computed with a wrong cosmology can significantly affect the posteriors’ amplitude. As an example, a covariance evaluated at Ω_m = 0.32 and σ₈ = 0.775 produces a 50% difference in the FoM with respect to the input-cosmology case, despite the small parameter shifts: ΔΩ_m = 0.013 and Δσ₈ = 0.053.

Fig. 3. Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, when using a fixed-cosmology (filled contours), a cosmology-dependent (empty contours), or the wrong (Ωm = 0.320 and σ8 = 0.775, dashed empty contours) covariance matrix. Lower left: Clustering and weak lensing masses alone. Upper right: Full-probe combination.

On the other hand, the impact of the cosmological dependence of the covariance decreases when combining clustering and counts (upper right panel). The improvement observed by using the cosmology-dependent clustering covariance disappears when incorporating the number counts likelihood since the latter already extracts all the additional information encoded in the clustering covariance. This confirms that there is no double counting of the halo mass function information. If such an effect were present, we would expect an additional tightening of the contours when combining number counts with a cosmology-dependent clustering covariance, compared to the fixed-covariance case. However, this is not observed. Further discussion on this can be found in EC:F24 and Fumagalli et al. (2024). Lastly, fixing the covariance matrix to the same wrong cosmology only produces a 6% difference in the FoM. In the following, we always use cosmology-dependent covariances.

4.4. Cluster clustering validation

In this section we examine the impact of photo-z distortions (Sect. 4.4.1), and secondary distortions (Sect. 4.4.2) on the 2PCF monopole, as well as on its covariance matrix. In Sect. 4.4.3, we assess how the choice of separation range in the 2PCF measurement influences the constraining power in the cosmological analysis.

4.4.1. Impact of photo-z

In the case of a photometric survey, the primary dilution in the 2PCF arises from photo-z uncertainties. As shown in the left panel of Fig. 4, the 2PCF shape is significantly altered as photo-z uncertainty increases, nearly erasing the BAO peak, which could potentially lead to a substantial loss of information. The right panel of Fig. 4 shows that the constraints on Ω_m and σ₈ from clustering and weak lensing masses are only modestly affected in the case of σ_z0 = 0.005, with a decrease in the FoM below 10% compared to the case without photo-z uncertainties, while a difference larger than 20% is obtained when considering σ_z0 = 0.01. Such differences slightly increase when the full probe combination is considered, as a consequence of the tightening of the posteriors (see Table 2). More considerations about the impact of photo-z uncertainties can be found in Appendix B.

Fig. 4. Impact of photo-z uncertainties on the 2PCF. Left: Prediction of the 2PCF for different values of photo-z uncertainty. The values are indicated as follows: σz0 = 0 in red, σz0 = 0.005 in orange, and σz0 = 0.01 in blue. This is for the redshift bin z = 0.4–0.8. Right: Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, from the combination of clustering and weak lensing masses for the three cases in the left panel.

4.4.2. Impact of secondary redshift distortions on clustering

As described in Sect. 2.2, the 2PCF is also affected by additional redshift distortions, namely RSDs, IR resummation, and geometric distortions. Although their effects are secondary to photo-z uncertainties, neglecting them can still lead to significant biases in the cosmological posteriors. The left panel of Fig. 5 illustrates the distortions that occur when each of these effects is omitted from the 2PCF modelling but included in the measurement. In the case of RSDs, such a difference translates into a significant bias in the Ω_m–σ₈ posteriors (right panel), with an IoI = 1.16, and a decrease of 5% in the FoM with respect to the full-modelling case (Table 3). Instead, ignoring the IR resummation or the geometrical distortion results in almost no impact on the cosmological FoM nor in the position of the contours.

Fig. 5. Systematic effects introduced by neglecting redshift distortions in 2PCF modelling. Left: Residuals from the 2PCF with respect to the full model (i.e. photo-z, RSDs, IR resummation, and geometric distortion correction) for the three cases. The line colours are as follows: the model without RSDs, blue line; model without IR resummation, red line; and model without geometric distortion correction, green line). Shaded areas are given by the square root of the diagonal covariance divided by the number of mocks. Right: Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, from the combination of clustering and weak lensing masses for the three cases represented in the left panel (same colour code) compared to the full redshift effects case (black empty contours). Dotted grey lines are the input cosmology of the catalogues.

When considering the full-probe analysis, the lack of RSDs modelling produces a shift in the posterior that is even larger (IoI = 2.18), due to the tightening of the contours. In all the other cases, the contours are well centred on the input cosmology. The difference in the FoM increases in all the cases, up to around 10%, indicating that neglecting each of the redshift effects may affect the final posteriors’ amplitude.

The 2PCF covariance is largely dominated by shot-noise (see the bottom panel of Fig. F.1 in EC:F24) and is therefore less affected by redshift distortions. As shown in the top panel of Fig. 6, RSDs slightly increase the covariance between radial bins, while photo-z uncertainties reduce it. Both effects are more pronounced at low redshift and diminish at higher redshift, where the cluster density decreases and shot-noise becomes more dominant.

Fig. 6. Residuals of the clustering covariance diagonal terms for the lowest (blue) and highest (red) redshift bins and λ ∈ [20, 30]. The grey area marks the 10% region. Top: Residuals of the distorted numerical covariance relative to the undistorted redshift case. Solid lines indicate the RSD case, while dashed lines include both RSDs and photo-z effects (σz0 = 0.01). Bottom: Residuals between analytical and numerical covariances for different cases. Solid lines are for true redshift, dashed for RSDs, and dotted for both RSDs and photo-z effects.

The lower panel of Fig. 6 compares the residuals between model predictions (Eq. 25) and numerical results for different redshift distortions. The similarity of residuals across cases suggests that RSDs and photo-z uncertainties do not alter the relationship between the model and numerical covariance. Thus, any differences between the model and the numerical matrix can be resolved by fitting the covariance parameters (defined in Sect. 2.2) as described in EC:F24 and Fumagalli et al. (2022), maintaining their validity independent of redshift effects.

4.4.3. Impact of radial scales in the clustering analysis

In this section we investigate the impact of radial scales on the cosmological analysis of cluster clustering. As shown in Fumagalli et al. (2024, their Appendix C), the majority of the cosmological information comes from the region around the BAO peak (approximately r = 60–130 h⁻¹ Mpc), with small and large scales offering negligible contributions. We repeat this analysis with a more detailed study, accounting for the properties of a Euclid-like survey, such as extended redshift range and larger sample. By taking as reference the separation range r = 20–130 h⁻¹ Mpc, Fig. 7 and Table 4 show that increasing the lower limit to larger scales, such as r = 40 h⁻¹ Mpc or r = 60 h⁻¹ Mpc, leads to a progressive broadening of the parameter constraints. This suggests that, for a survey such as Euclid, scales between 20 and 60 h⁻¹ Mpc carry sizeable cosmological information. On the other hand, we confirm that not much more information is extracted when extending the 2PCF analysis to larger scales, with a difference in the FoM of 3% when increasing the upper limit to 200 h⁻¹ Mpc. Similar results are obtained from the full-probe analysis.

Fig. 7. Marginalised posterior distributions with 68% and 95% confidence intervals from combined analysis of clustering and weak lensing masses for different clustering scales: r = 20–200 h−1 Mpc (in black), r = 20–130 h−1 Mpc (in blue), r = 40–130 h−1 Mpc (in red), and r = 60–130 h−1 Mpc (in green). Dotted grey lines are the input cosmology of the catalogues.

4.5. Number count validation

This section examines the effects of RSDs and photo-z redshift errors on cluster number counts and its covariance. The presence of redshift effects, especially photo-z uncertainties, shifts redshifts from the centre of the distribution to the edges, potentially altering the number of objects in each redshift bin. However, this does not affect the overall shape of the cluster abundance distribution, leaving the shot-noise term unchanged. In contrast, the sample variance signal, which is determined by the clustering signal integrated over the redshift bins (Eq. 10), reflects the distortions in the power spectrum caused by RSDs and photo-z uncertainties. As shown in Fig. 8, RSDs tend to amplify sample variance, whereas photo-z uncertainties reduce it. The increase due to RSDs arises because large-scale velocity flows, which drive RSDs on linear scales, are predominantly sourced by long-wavelength (small-k) density fluctuations – the same modes that contribute to sample variance. As a result, RSDs enhance the impact of these fluctuations, leading to a higher sample variance. Conversely, photo-z uncertainties blur the positions of clusters along the line of sight, effectively smoothing the density field in redshift space. This reduces the contrast of density fluctuations, leading to a decrease in sample variance.

Fig. 8. Number count sample covariance terms from numerical simulations for different redshift settings: undistorted (red lines), RSDs (purple lines), RSDs and small photo-z uncertainty (σz0 = 0.005, blue lines), RSDs and large photo-z uncertainty (σz0 = 0.01, green lines), and redshift with only large photo-z uncertainty (σz0 = 0.01, yellow lines). Solid lines represent the diagonal sample variance, dashed lines are the sample covariance between two adjacent mass bins (first and second) and same redshift bin, and dotted lines are the sample covariance between adjacent redshift bins and same mass bin (first one).

Incorporating in the model the RSDs and photo-z effects solely through the monopole of the power spectrum proves insufficient to recover the numerical results, as shown in Fig. 9. On the other hand, including the correct modelling of the anisotropic power spectrum in Eq. (10) would break many symmetries and lead to a computationally expensive calculation. Dashed lines in Fig. 9 show that an effective power spectrum given simply by the sum of the power spectrum multipoles, i.e.

$$\begin{array}{c}\hfill P_{\text{eff}}(k)={\displaystyle \sum _{\ell =0,2,4}}{P}_{\ell }(k),\end{array}$$(31)

Fig. 9. Number count covariance terms: numerical matrix (shaded areas, representing the 1σ region), analytical model with power spectrum monopole (solid lines), and model with effective power spectrum (see Eq. (31), dashed lines). The colour-coded terms represent different components: Grey is for diagonal shot-noise; blue is for diagonal sample variance; red and yellow are for the first and second off-diagonal sample covariance between mass bins, respectively; and green is for the first off-diagonal sample covariance in redshift bins.

effectively captures the sample covariance in a computationally efficient manner, balancing accuracy with practical feasibility. Furthermore, since the hexadecapole (ℓ = 4) is found to be negligible, the sum can be restricted to the monopole and quadrupole terms.

We find that neglecting photo-z and RSD effects in the covariance significantly impacts cosmological constraints, as shown in Fig. 10. In the photo-z case with σ_z0 = 0.01, this leads to an underestimation of the posterior amplitude by about 20%. The effect is expected to be even larger for smaller photo-z uncertainties, as the difference in sample variance increases. A similar trend is observed in the full-probe analysis, where neglecting RSDs and photo-z in the number count covariance degrades the FoM by 10% compared to proper covariance modelling.

Fig. 10. Marginalised posterior distributions with 68% and 95% confidence intervals from number counts (blue and cyan contours) and number counts with clustering (red and orange contours), with (empty contours) and without (filled contours) the proper modelling of redshift distortions in the number count covariance. Dotted grey lines are the input cosmology of the catalogues.

5. Discussion and conclusions

In this paper, we have shown a framework to extract cosmological information from the combination of number counts and clustering of galaxy clusters in a survey with a mass selection, sky coverage, and depth similar to those expected from Euclid. Building on the models from F21 and EC:F24, we extended the methodology to explicitly include the effect of RSDs as well as observational uncertainties, such as errors in photo-z and in the mass-richness relation. We assessed the impact of including these effects on number counts, the monopole of the 2PCF, and their respective covariance matrices.

Our validation was conducted using 1000 Euclid-like light cones generated with the PINOCCHIO code, where we introduced observational effects such as photometric and richness estimate errors. To quantify the influence of these distortions, we performed a likelihood analysis across different configurations and evaluated their impact on cosmological posteriors. The analysis focused on constraining Ω_m, σ₈, and the four scaling relation parameters (Eqs. 6 and 7).

The key findings of our analysis can be summarised as follows.

• We demonstrated in Sect. 4.1 that numerical covariances show no significant cross-correlation between number counts and 2PCF, with a Pearson correlation factor consistent with zero. The two statistics are therefore independent. Such a result is in agreement with Mana et al. (2013) and Fumagalli et al. (2024).

• We find that the full-probe analysis is the most powerful combination for extracting cosmological information and the less affected by prior-volume effects.

• We confirmed in Sect. 4.3 that the cosmology-dependence of the clustering covariance brings useful information, increasing the constraining power of clustering by about 40%. Such an improvement disappears when number counts are combined, which is in line with the result that number counts already encapsulate additional information from the clustering covariance. Similarly, evaluating the covariance at a wrong cosmology has a significant impact on the clustering-only case, an effect that decreases when adding number counts.

• The primary effect of redshift effects on the 2PCF arises from photo-z uncertainties. We showed in Sect. 4.4.1 that an optimistic estimate for photo-z uncertainties of σ_z0 = 0.005 results in a small impact on parameter constraints, reducing the cosmological FoM by less than 10%. However, for a more conservative value of σ_z0 = 0.01, the posteriors broaden significantly, with a decrease in the FoM of 20–30%.

• We showed in Sect. 4.4.2 that secondary redshift distortions – namely, RSDs, IR resummation, and geometric distortions – have an almost negligible impact (≲5%) on the FoM when clustering is considered alone. This impact slightly increases to around 10% when combining all probes in the analysis. In addition, the absence of proper RSD modelling leads to a non-negligible shift in the contours (about 2σ).

• Although most of the cosmological information in the linear regime is contained at large scales (r = 60–130 h⁻¹ Mpc), we showed in Sect. 4.4.3 that smaller scales also contribute to increases in the constraining power of clustering. In contrast, expanding the analysis to scales larger than 130 h⁻¹ Mpc, does not bring any significant improvement.

• While RSDs and photo-z uncertainties have no direct impact on the number counts themselves, we found that they do influence the number count covariance (Sect. 4.5). Specifically, RSDs increase the sample variance, whereas photo-z uncertainties decrease it. Neglecting these effects in the covariance modelling can lead to a 10–20% underestimation of the posterior amplitude. We propose an effective model that accurately incorporates these effects into the covariance without adding unnecessary computational burden to the model.

Our work underscores the benefits of combining cluster clustering with traditional number counts. First, we note that the use of a wide redshift range allows the amount of extracted information to be greatly increased. This is observed, in particular, through the combination of number counts and clustering. In Fumagalli et al. (2024), where a single low-redshift bin was analysed, the combination of these two probes failed to constrain σ₈. In contrast, with multiple bins spanning a wider redshift range, the constraints become much tighter and almost competitive with the combination of number counts and weak lensing mass, highlighting how the redshift dependence of the halo mass function and halo bias provides substantial additional information.

Second, in addition to confirming that the combination of the three probes provides the most stringent constraints – with an FoM improvement of more than 300% with respect to the standard combination of number counts and weak lensing masses – we also find that these results are better centred around the input cosmology. This indicates that the addition of clustering not only enhances the extraction of cosmological information but also mitigates the impact of systematics, leading to more precise and more accurate constraints. The absence of significant correlations between number counts and clustering simplifies the analysis, further optimising the information gain. Moreover, this combination is less sensitive to the cosmological dependence of the covariances, as the cosmological dependence of the observables saturates the information that can be extracted from the sample.

However, our findings also highlight the need for caution when performing full-probe analyses, as narrower parameter constraints make them more sensitive to systematic effects. Some uncertainties that are negligible when considering clustering or number counts independently may become significant when both probes are combined, stressing the importance of carefully modelling systematic uncertainties.

It is important to emphasise that the results of this work do not fully reflect the complexity of the future Euclid cluster-cosmology analysis. The presented analysis is rather conservative, assuming mass calibration at the 1% level out to redshift z = 2 and neglecting potential model systematics, such as uncertainties in the P(λ^ob|M) relation at high redshift. In particular, the impact of projection effects on cluster number counts, clustering, and their combination remains to be quantified. In addition, some key elements still need to be explored. Firstly, current analyses neglect the cross-correlation between weak lensing mass and both counts and clustering. While this assumption is supported by reasonable arguments, it needs to be tested. Additionally, the impacts of the survey mask and the presence of different homogeneous regions on the covariances of clustering and counts still require investigation. Finally, the clustering analysis is not fully explored due to the absence of higher-order multipoles. Given that photo-z uncertainties induce a non-negligible signal in these components, their inclusion could significantly reduce the final uncertainties on the parameters. Nevertheless, these results provide valuable insight into the influence of observational and modelling systematics on cosmological constraints from upcoming Euclid data, highlighting the necessity of accounting for effects that were neglected in previous studies.

Acknowledgments

AF acknowledges support by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311 and from the Ludwig-Maximilians-Universität in Munich. This paper is supported by the Agenzia Spaziale Italiana (ASI) under – Euclid-FASE D Attivita’ scientifica per la missione – Accordo attuativo ASI-INAF n. 2018-23-HH.0, by the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 1409 published on 14.9.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU– Project Title “Space-based cosmology with Euclid: the role of High-Performance Computing” – CUP J53D23019100001 – Grant Assignment Decree No. 962 adopted on 30/06/2023 by the Italian Ministry of Ministry of University and Research (MUR); by the Italian Research centre on High-Performance Computing Big Data and Quantum Computing (ICSC), a project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2, by the INFN INDARK PD51 grant, and by the PRIN 2022 project EMC2 – Euclid Mission Cluster Cosmology: unlock the full cosmological utility of the Euclid photometric cluster catalogue (code no. J53D23001620006). MR and FM acknowledge the financial contribution from the PRIN-MUR 2022 20227RNLY3 grant “The concordance cosmological model: stress-tests with galaxy clusters” supported by Next Generation EU and from the grant ASI n. 2024-10-HH.0 “Attività scientifiche per la missione Euclid – fase E”. ET acknowledges support by STFC through Imperial College Astrophysics Consolidated Grant ST/W000989/1. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org).

References

Appendix A: Impact of redshift photometric bias

We check here the impact of the photo-z bias on cosmological constraints. For this test, we modified the redshift selection function P(z^ob | z^tr, Δλ_i^ob) in Eq. (2) by adding a shift to the true redshift of Δz^tr = 0.002. Figure A.1 shows that, even in the optimistic photo-z case (σ₀ = 0.005), such a redshift bias does not impact the cosmological constraints significantly. Both the clustering and weak lensing masses and the full-probe combinations are characterised by no difference in the IoI between the bias and no-bias cases, as well as by a ΔFoM consistent with zero. We therefore do not add any bias to the observed redshift in out simulations (see Sect. 3.1), as the effect is negligible.

Fig. A.1. Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, with (filled contours), and without (empty contours), photo-z bias. The results correspond to the optimistic photo-z uncertainty amplitude. Lower left: Clustering and weak lensing masses alone. Upper right: Full-probe combination.

Appendix B: Impact of non-linearities

When including the effect of photo-z uncertainties, we notice a discrepancy between our model and the measured 2PCF. Although the model outlined in Sect. 2.2 accounts for photo-z-induced distortions in the 2PCF shape, Fig. B.1 highlights a scale-dependent mismatch between the predicted and observed 2PCF when photo-z uncertainties are included. The effect, which increases with redshift, becomes significant for z ≳ 1. A similar discrepancy characterises the undistorted 2PCF, but is confined to the non-linear regime (r ≲ 20 h⁻¹ Mpc). Photo-z uncertainties mix the power spectrum modes parallel to the line of sight, propagating non-linearities to larger scale and making the linear theory inadequate for describing clustering up to scales of approximately 60 h⁻¹ Mpc at high redshift. Non-linear distortions also explain the redshift dependence of this mismatch: Objects of fixed mass are characterised by a higher bias and, therefore, increasing non-linearity at higher redshifts.

Fig. B.1. Two-point correlation function without (red) and with (blue) photo-z uncertainty (σz0 = 0.01). Solid lines are the linear model, and circles are the measurements, averaged over the 1000 mocks. The standard error on the mean is included but too small to be noticeable. The grey area indicates scales not included in the cosmological analysis. Top: Low-redshift bin. Bottom: High-redshift bin.

Addressing this requires a non-linear model for the halo bias, which is beyond the scope of this paper. Instead, we introduce a correction factor, defined as

$$\begin{array}{c}\hfill B_2=\frac{{\widehat{\xi }}_{\mathrm{obs}}}{{\xi }_{\mathrm{fid}}},\end{array}$$(B.1)

where ${\widehat{\xi }}_{\mathrm{obs}}$ is the measured 2PCF, ξ_fid is the prediction at the fiducial cosmology, and B₂ is a scale-dependent correction factor for each redshift and richness bin. This factor is included in the likelihood analysis to enhance the model accuracy. However, it serves primarily as a ‘safety net’ and is not strictly essential: as shown in Fig. B.2, this discrepancy does not lead to significant differences in the cosmological posteriors.

Fig. B.2. Parameter posteriors with 68% and 95% confidence intervals, from clustering and weak lensing mass analysis, with (blue contours) and without (empty red contours) the bias correction factor B2.

With low-redshift and large scales dominating the clustering constraining power, the high-redshift, small-scale regime contributes only marginally to the cosmological constraints.

We also notice a small discrepancy between model and measurements around the BAO peak, when true redshifts are considered (top panel of Fig. B.3). Such a difference increases when RSDs are included, as shown in the central panel. Again, the cause can be attributed to the failure of linear theory in IR re-summation and RSDs modelling. However, as can be seen from the bottom panel, in the presence of photo-z uncertainties this effect is completely washed out.

Fig. B.3. Impact of non-linearities around the BAO peak. Model (solid lines) and measured (circles) 2PCF with true redshift (top), true redshift with RSDs (centre), and observed redshift with RSDs and photo-z (top).

All Tables

All Figures

Fig. 1. Cross-covariance between number counts and clustering. For better visualisation, we used redshift bins of width Δz = 0.2 for counts and Δz = 0.5 and a single richness bin for the 2PCF. Left: Auto- and cross-correlation matrix of number counts and clustering, computed from 1000 mocks. Right: Log-likelihood residuals for number counts and clustering for each one of the 1000 light cones, with respect to the mean value assuming the fiducial model parameters.

In the text

	Fig. 2. Parameter posterior distributions with 68% and 95% confidence intervals from different combinations of probes: number counts and weak lensing mass (in blue), clustering and weak lensing mass (in orange), number counts and clustering (in green), and all the three probes (in red). Dotted grey lines are the input cosmology of the catalogues.
In the text

	Fig. 3. Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, when using a fixed-cosmology (filled contours), a cosmology-dependent (empty contours), or the wrong (Ωm = 0.320 and σ8 = 0.775, dashed empty contours) covariance matrix. Lower left: Clustering and weak lensing masses alone. Upper right: Full-probe combination.
In the text

Fig. 4. Impact of photo-z uncertainties on the 2PCF. Left: Prediction of the 2PCF for different values of photo-z uncertainty. The values are indicated as follows: σz0 = 0 in red, σz0 = 0.005 in orange, and σz0 = 0.01 in blue. This is for the redshift bin z = 0.4–0.8. Right: Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, from the combination of clustering and weak lensing masses for the three cases in the left panel.

In the text

Fig. 5. Systematic effects introduced by neglecting redshift distortions in 2PCF modelling. Left: Residuals from the 2PCF with respect to the full model (i.e. photo-z, RSDs, IR resummation, and geometric distortion correction) for the three cases. The line colours are as follows: the model without RSDs, blue line; model without IR resummation, red line; and model without geometric distortion correction, green line). Shaded areas are given by the square root of the diagonal covariance divided by the number of mocks. Right: Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, from the combination of clustering and weak lensing masses for the three cases represented in the left panel (same colour code) compared to the full redshift effects case (black empty contours). Dotted grey lines are the input cosmology of the catalogues.

In the text

Fig. 6. Residuals of the clustering covariance diagonal terms for the lowest (blue) and highest (red) redshift bins and λ ∈ [20, 30]. The grey area marks the 10% region. Top: Residuals of the distorted numerical covariance relative to the undistorted redshift case. Solid lines indicate the RSD case, while dashed lines include both RSDs and photo-z effects (σz0 = 0.01). Bottom: Residuals between analytical and numerical covariances for different cases. Solid lines are for true redshift, dashed for RSDs, and dotted for both RSDs and photo-z effects.

In the text

	Fig. 7. Marginalised posterior distributions with 68% and 95% confidence intervals from combined analysis of clustering and weak lensing masses for different clustering scales: r = 20–200 h−1 Mpc (in black), r = 20–130 h−1 Mpc (in blue), r = 40–130 h−1 Mpc (in red), and r = 60–130 h−1 Mpc (in green). Dotted grey lines are the input cosmology of the catalogues.
In the text

Fig. 8. Number count sample covariance terms from numerical simulations for different redshift settings: undistorted (red lines), RSDs (purple lines), RSDs and small photo-z uncertainty (σz0 = 0.005, blue lines), RSDs and large photo-z uncertainty (σz0 = 0.01, green lines), and redshift with only large photo-z uncertainty (σz0 = 0.01, yellow lines). Solid lines represent the diagonal sample variance, dashed lines are the sample covariance between two adjacent mass bins (first and second) and same redshift bin, and dotted lines are the sample covariance between adjacent redshift bins and same mass bin (first one).

In the text

Fig. 9. Number count covariance terms: numerical matrix (shaded areas, representing the 1σ region), analytical model with power spectrum monopole (solid lines), and model with effective power spectrum (see Eq. (31), dashed lines). The colour-coded terms represent different components: Grey is for diagonal shot-noise; blue is for diagonal sample variance; red and yellow are for the first and second off-diagonal sample covariance between mass bins, respectively; and green is for the first off-diagonal sample covariance in redshift bins.

In the text

	Fig. 10. Marginalised posterior distributions with 68% and 95% confidence intervals from number counts (blue and cyan contours) and number counts with clustering (red and orange contours), with (empty contours) and without (filled contours) the proper modelling of redshift distortions in the number count covariance. Dotted grey lines are the input cosmology of the catalogues.
In the text

	Fig. A.1. Marginalised posterior distributions in the Ωm–σ8 plane, with 68% and 95% confidence intervals, with (filled contours), and without (empty contours), photo-z bias. The results correspond to the optimistic photo-z uncertainty amplitude. Lower left: Clustering and weak lensing masses alone. Upper right: Full-probe combination.
In the text

Fig. B.1. Two-point correlation function without (red) and with (blue) photo-z uncertainty (σz0 = 0.01). Solid lines are the linear model, and circles are the measurements, averaged over the 1000 mocks. The standard error on the mean is included but too small to be noticeable. The grey area indicates scales not included in the cosmological analysis. Top: Low-redshift bin. Bottom: High-redshift bin.

In the text

	Fig. B.2. Parameter posteriors with 68% and 95% confidence intervals, from clustering and weak lensing mass analysis, with (blue contours) and without (empty red contours) the bias correction factor B2.
In the text

	Fig. B.3. Impact of non-linearities around the BAO peak. Model (solid lines) and measured (circles) 2PCF with true redshift (top), true redshift with RSDs (centre), and observed redshift with RSDs and photo-z (top).
In the text