© The Authors 2024

1. Introduction

Clusters of galaxies are known to be powerful probes for studying the geometry and evolution of the large-scale structure of the Universe. Originating from high-density regions in the initial matter density field, they grow through a hierarchical process of accretion and merging of small objects into increasingly larger systems, and they currently represent the most massive virialized objects in the Universe (for reviews, see e.g., Allen et al. 2011; Kravtsov & Borgani 2012). Over the past two decades, the statistical tool that was most frequently adopted to extract cosmological information from cluster catalogs is the study of their abundance in mass, namely the number counts (e.g., Borgani et al. 2001; Vikhlinin et al. 2009; Planck Collaboration XX 2014; Bocquet et al. 2016; Costanzi 2021). This allows us to derive constraints on the average matter density (Ω_m) and the amplitude of density fluctuations (σ₈) in the Universe. Several examples can be provided about cosmological constraints from cluster-count analyses from optical cluster catalogs (e.g., Costanzi et al. 2019; Dark Energy Survey Collaboration 2022; Lesci et al. 2022b), X-ray catalogs (e.g., Mantz et al. 2010, 2015; Schellenberger & Reiprich 2017; Chiu et al. 2023) or millimeter catalogs obtained through the Sunyaev-Zeldovich signature (e.g., Planck Collaboration XXIV 2016; Bocquet et al. 2019), and also from the combination of multiwavelength data (e.g., Costanzi 2021) or different cosmological probes, such as CMB (Salvati et al. 2022).

The study of the clustering properties of galaxy clusters provides another means for their cosmological exploitation (Borgani et al. 1999; Moscardini et al. 2000; Estrada et al. 2009; Marulli et al. 2018). This probe offers several advantages: The fact that clusters are highly biased tracers of the matter density field, their detectable clustering signal on large scales where linear theory is applicable, the influence of cosmological parameters on the bias-mass relation (Mo & White 1996; Tinker et al. 2010), and the minimum impact of baryonic effects on clustering statistics (Castro et al. 2020) collectively contribute to the enhanced constraining power and reliability of cluster clustering as a valuable tool for cosmological investigations. Although the available statistics is still too low to use the clustering of clusters as a competitive stand-alone probe, it has been shown to be a valuable source of information for breaking parameter degeneracies when combined with other observables such as cluster counts (Schuecker et al. 2003; Majumdar & Mohr 2004; Sereno et al. 2015; Sartoris et al. 2016), and it improves the calibration of mass-observable relations (Mana et al. 2013; To et al. 2021b; Lesci et al. 2022a).

While cluster abundance and halo bias can be predicted with good accuracy as a function of mass from first principles or simulations (Mo & White 1996; Sheth & Tormen 1999; Tinker et al. 2010; Euclid Collaboration 2023), cluster masses are not directly observable. The mass inference process must in fact rely on observational proxies that exhibit correlations with mass, such as properties of galaxies (richness or velocity dispersion) or of the gas within the cluster (total gas mass, temperature, or pressure). The scaling relations linking these observables to cluster masses provide a statistical estimate of the latter, but require a careful calibration to avoid biases in the cosmological inference (Kravtsov & Borgani 2012; Pratt et al. 2019). Weak gravitational lensing measurements are based on the estimate of the tangential alignment of background galaxies around foreground clusters due to gravitational lensing. They provide the gold-standard technique for estimating cluster masses (Johnston et al. 2007; Hoekstra et al. 2013; von der Linden et al. 2014; Murata et al. 2018; Simet et al. 2017; Melchior et al. 2017). The advantages in using weak-lensing mass information is that weak-lensing measurements capture both dark and baryonic matter, and do not rely on assumptions about the dynamical state of the cluster. However, the calibration process is challenging because of various biases and observational uncertainties. Systematic errors, including shear and photometric redshift biases, halo triaxiality, and projection effects, hamper the interpretation of the weak-lensing measurements and dominate the overall error budget of the mass calibration (von der Linden et al. 2014; Hoekstra et al. 2015; Simet et al. 2017; Miyatake et al. 2019; Wu et al. 2022).

Cosmological0]Please provide the missing City and Country name for the Affiliation 6. constraints obtained from cluster counts and cluster clustering, as well as the accuracy of the mass calibration from weak lensing, are expected to improve strongly with the most recent and forthcoming wide-surveys such as the Dark Energy Survey¹ (DES, Dark Energy Survey Collaboration 2005), eROSITA² (Predehl 2014), Euclid³ (Laureijs et al. 2011) or the Vera C. Rubin Observatory LSST⁴ (LSST Dark Energy Science Collaboration 2012). Nevertheless, the cosmological gain provided by the clustering statistics has not been thoroughly investigated in currently available cluster catalogs. For this purpose, we derive cosmological constraints from the combination of abundances, weak lensing, and clustering of galaxy clusters from the redMaPPer cluster catalog in the Sloan Digital Sky Survey data release 8 (SDSS DR8, Aihara et al. 2011). The analysis of the cluster abundances and the weak lensing is based on Costanzi et al. (2019, hereafter C19). The aim is to show that the inclusion of clustering helps to narrow down the uncertainties on the inference of cosmological parameters from cluster counts and weak-lensing masses, and it also reduces biases and systematic uncertainties in the mass calibration.

This paper is structured as follows. In Sect. 2 we present the theoretical framework within which we modeled the three statistics considered here, that is, cluster counts, cluster clustering, and weak-lensing masses, and we also provide the likelihood models we adopted to infer the cosmological and mass-observable relation parameters. In Sect. 3 we describe the data we analyzed in this work. In Sect. 4 we present the results, and in Sect. 5 we discuss our conclusions.

2. Theory and methods

In this section, we introduce the theoretical formalism to describe number counts, weak-lensing masses, and cluster clustering, along with their covariance matrices. We also describe the likelihood function adopted for the parameter inference.

2.1. Cluster number counts

We model the cluster number counts in the ith redshift bin and jth richness bin as

$$\begin{array}{cc}& ⟨N(\mathrm{\Delta }{z}_i^{\mathrm{ob}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}})⟩={\displaystyle {\int }_0^{\infty }\mathrm{d}{z}^{\mathrm{tr}}{\mathrm{\Omega }}_{\mathrm{mask}}(z^{\mathrm{tr}})\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}(z^{\mathrm{tr}})}\hfill \\ \hfill & \times ⟨n(z^{\mathrm{tr}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}})⟩{\displaystyle {\int }_{\mathrm{\Delta }{z}_i^{\mathrm{ob}}}\mathrm{d}{z}^{\mathrm{ob}}P(z^{\mathrm{ob}}|z^{\mathrm{tr}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}}),}\hfill \end{array}$$(1)

where Ω_mask is the redshift-dependent survey area (see C19), and dV/dΩ dz is the comoving volume element per unit redshift and solid angle. $P(z^{\text{ob}}|z^{\text{tr}},\mathrm{\Delta }{\lambda }_j^{\text{ob}})$ is the probability distribution of assigning an observed redshift z^ob given the true redshift z^tr of a cluster and observed richness, to account for the observational scatter due to the photometric uncertainty. As described in C19, it is modeled as a Gaussian distribution with a mean z^tr and a redshift and richness-dependent variance. Last, the term $\langle n(z^{\text{tr}},\mathrm{\Delta }{\lambda }_j^{\text{ob}})\rangle$ describes the expected number density of halos in the ith richness bin, and it is given by

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

where dn/dM is the halo-mass function, and P(λ^ob | M, z^tr) is 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}$$(3)

The term P(λ^tr | M, z^tr) is the redshift-dependent intrinsic richness-mass relation. Following C19, we modeled it as a skew-normal distribution with a mean and variance given by

$$\begin{array}{cc}\hfill ⟨{\lambda }^{\mathrm{tr}}|M⟩& ={\left(\frac{M-M_{\min }}{M_1-M_{\min }}\right)}^{\alpha },\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\sigma }_{ln{\lambda }^{\mathrm{tr}}}^2& ={\sigma }_{\mathrm{intr}}^2+\frac{⟨{\lambda }^{\mathrm{tr}}|M⟩-1}{{⟨{\lambda }^{\mathrm{tr}}|M⟩}^2}.\hfill \end{array}$$(5)

Here, M_min is the minimum mass for a halo to form a central galaxy, while M₁ is the characteristic mass at which halos acquire one satellite galaxy. The variance is composed of an intrinsic scatter, σ_intr, plus a Poisson contribution.

The other term in Eq. (3), that is, P(λ^ob | λ^tr, z^tr), represents the observational scatter in the richness-mass relation due to photometric noise, uncertainties in the background subtraction, and projection or percolation effects. Projection effects occur when multiple foreground and background objects along the same line-of-sight are mistakenly associated with a galaxy cluster, which increases the apparent richness of the cluster. As a consequence, percolation is the reduction in the count of member galaxies for clusters that are masked by these projection effects. Following C19, the background subtraction and photometric noise terms are modeled as Gaussian components, while projection and percolation effects are modeled as an exponential and a uniform distribution, respectively. The corresponding model parameters were calibrated by means of a real and a simulated data analysis (see C19 for details).

The covariance matrix associated with cluster number counts was computed analytically following the model by Hu & Kravtsov (2003), validated in Euclid Collaboration (2021). This model accounts for the effects of shot-noise and sample variance, and in the case of a single redshift bin, it is given by

$$\begin{array}{c}\hfill C^{\mathrm{NC}}(\mathrm{\Delta }{z}_i^{\mathrm{ob}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}},\mathrm{\Delta }{\lambda }_k^{\mathrm{ob}})={\langle N\rangle }_{\mathit{ij}}{\delta }_{\mathit{jk}}+{\langle Nb\rangle }_{\mathit{ij}}{\langle Nb\rangle }_{\mathit{ik}}{\sigma }^2(\mathrm{\Delta }{z}_i^{\mathrm{ob}}),\end{array}$$(6)

where ⟨N⟩_ij is the prediction for the number counts in the ith redshift bin and jth richness bin (Eq. (1)), and ⟨Nb⟩_ij is the prediction for the number counts times the halo bias, computed in an analogous way. ${\sigma }^2(\mathrm{\Delta }{z}^{\text{ob}}{}_i)$ is the variance of the matter linear density field,

$$\begin{array}{c}\hfill {\sigma }^2(\mathrm{\Delta }{z}_i^{\mathrm{ob}})={\displaystyle \int }\frac{{\mathrm{d}}^{3}k}{{(2\pi )}^3}{P}_{\mathrm{m}}(k,\mathrm{\Delta }{z}_i^{\mathrm{ob}}){|W_{\mathrm{\Delta }{z}_j^{\mathrm{ob}}}(\mathbf{k})|}^2,\end{array}$$(7)

where $W_{\mathrm{\Delta }{z}^{\text{ob}}{}_j}$ is the window function of the light-cone redshift slice (C19).

In addition, we considered a contribution due to the uncertainty in the miscentering corrections, as described in C19 (see their Appendix A).

2.2. Weak-lensing masses

Following C19, we relied for the mass calibration on the mean cluster mass measurements from the stacked shear analysis of Simet et al. (2017) and did not model the cluster shear profiles directly. Similarly to the number counts, the expectation value of the mean cluster mass within the ith redshift bin and jth richness bin is given by

$$\begin{array}{c}\hfill \overline{M}(\mathrm{\Delta }{z}_i^{\mathrm{ob}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}})=\frac{⟨M^{\mathrm{tot}}(\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}},\mathrm{\Delta }{z}_i^{\mathrm{ob}})⟩}{⟨N(\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}},\mathrm{\Delta }{z}_i^{\mathrm{ob}})⟩},\end{array}$$(8)

where ⟨M^tot⟩ is the total mass associated with clusters identified in given redshift and richness intervals, given by

$$\begin{array}{cc}& ⟨M^{\mathrm{tot}}(\mathrm{\Delta }{z}_i^{\mathrm{ob}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}})⟩={\displaystyle {\int }_0^{\infty }\mathrm{d}{z}^{\mathrm{tr}}{\mathrm{\Omega }}_{\mathrm{mask}}(z^{\mathrm{tr}})\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}(z^{\mathrm{tr}})}\hfill \\ \hfill & \times ⟨Mn(z^{\mathrm{tr}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}})⟩{\displaystyle {\int }_{\mathrm{\Delta }{z}_i^{\mathrm{ob}}}\mathrm{d}{z}^{\mathrm{ob}}P(z^{\mathrm{ob}}|z^{\mathrm{tr}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}}),}\hfill \end{array}$$(9)

with

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

2.3. Cluster clustering

We modeled the two-point correlation function in the ath radial bin and ith redshift bin as (Matarrese et al. 1997)

$$\begin{array}{c}\hfill ⟨\xi (\mathrm{\Delta }{r}_a,\mathrm{\Delta }{z}_i^{\mathrm{ob}})⟩={\displaystyle \int }\frac{\mathrm{d}k{k}^2}{2{\pi }^2}{\langle \overline{b}\sqrt{P_{\mathrm{m}}(k)}\rangle }_i^2{W}_a(k),\end{array}$$(11)

where

$$\begin{array}{cc}& \langle \overline{b}\sqrt{P_{\mathrm{m}}(k)}{\rangle }_i=\frac{1}{⟨N(\mathrm{\Delta }{z}_i^{\mathrm{ob}})⟩}{\displaystyle {\int }_0^{\infty }\mathrm{d}{z}^{\mathrm{tr}}{\mathrm{\Omega }}_{\mathrm{mask}}(z^{\mathrm{tr}})\frac{\mathrm{d}V}{\mathrm{d}\mathrm{\Omega }\mathrm{d}z}(z^{\mathrm{tr}})}\hfill \\ \hfill & \times n(z^{\mathrm{tr}})\overline{b}(z^{\mathrm{tr}})\sqrt{P_{\mathrm{m}}(k,z^{\mathrm{tr}})}{\displaystyle {\int }_{\mathrm{\Delta }{z}_i^{\mathrm{ob}}}\mathrm{d}{z}^{\mathrm{ob}}P(z^{\mathrm{ob}}|z^{\mathrm{tr}},\mathrm{\Delta }{\lambda }_j^{\mathrm{ob}}).}\hfill \end{array}$$(12)

In the above equation, $\overline{b}(z^{\text{tr}})$ is the effective bias for all clusters with an observed richness above a given threshold value ${\lambda }_{\text{th}}^{\text{ob}}$,

$$\begin{array}{c}\hfill \overline{b}(z^{\mathrm{tr}})=\frac{1}{n(z^{\mathrm{tr}})}{\displaystyle {\int }_0^{\infty }\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M,z^{\mathrm{tr}})b(M,z^{\mathrm{tr}}){\int }_{{\lambda }_{\mathrm{th}}^{\mathrm{ob}}}^{\infty }\mathrm{d}{\lambda }^{\mathrm{ob}}P({\lambda }^{\mathrm{ob}}|M,z^{\mathrm{tr}}).}\end{array}$$(13)

Similarly, $N(\mathrm{\Delta }{z}_i^{\text{ob}})$ and n(z^tr) are given by Eqs. (1) and (2), with integrals over the richness bins substituted with integrals above the richness threshold⁵.

Finally, the term W_a(k) in Eq. (11) is the spherical shell window function, given by

$$\begin{array}{c}\hfill W_a(k)={\displaystyle \int }\frac{{\mathrm{d}}^{3}r}{V_a}{j}_0(kr)=\frac{r_{a,+}^3{W}_{\mathrm{th}}(k{r}_{a,+})-r_{a,-}^3{W}_{\mathrm{th}}(k{r}_{a,-})}{r_{a,+}^3-r_{a,-}^3},\end{array}$$(14)

where the W_th(kr) is the top-hat window function, V_a is the volume of the ath spherical shell, and r_a, −,  r_a, + are the extremes of the separation bin.

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

$$\begin{array}{c}\hfill P_{\mathrm{m}}^{\prime }(k)=P_{\mathrm{m}}(k)\frac{\sqrt{\pi }}{2k\sigma }\mathrm{erf}(k\sigma ),\end{array}$$(15)

where σ depends on the photo-z error σ_z as

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

Photo-z errors perturb the power spectrum in the same way as redshift-space distortions. At this level of statistics, the former represent the dominant error source, so that redshift-space distortions can safely be neglected (Veropalumbo et al. 2014; Sereno et al. 2015). We verified that the high-order correction terms (Sereno et al. 2015) have a negligible impact on our cosmological constraints. The effect of the photo-z uncertainty is to decrease the correlation on small scales and to increase it at large scales, at about the BAO peak. This allows us to neglect the infrared resummation as well (IR, Senatore & Zaldarriaga 2015; Baldauf et al. 2015), which is otherwise needed to correct for the broadening and shift of the peak due to nonlinear damping.

Last. when measuring the two-point correlation function from the data, it is necessary to make assumptions about the cosmology to convert redshifts into distances. However, because it is not feasible to assume the true cosmology in the measurements, a parameter is introduced in the model (specifically, into Eq. (14)) to account for geometric distortions in the model (Marulli et al. 2012, 2016),

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

with

$$\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}$$(18)

where 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 matrix by applying the semi-analytical model presented and validated in Euclid Collaboration (2024, see Eq. (12)). The model consists of a Gaussian covariance plus a low-order non-Gaussian term, and is given by

$$\begin{array}{cc}& C^{\mathrm{CL}}(\mathrm{\Delta }{z}_i^{\mathrm{ob}},\mathrm{\Delta }{r}_a,\mathrm{\Delta }{r}_b)=\hfill \\ \hfill & \frac{2}{V_i}{\displaystyle \int }\frac{\mathrm{d}k{k}^2}{2{\pi }^2}{\langle {\overline{b}}^2{P}_{\mathrm{m}}(k)+\frac{1}{\overline{n}}\rangle }_i^2{W}_a(k)W_b(k)\hfill \\ \hfill & +\frac{2}{V_a{V}_i}{\displaystyle \int }\frac{\mathrm{d}k{k}^2}{2{\pi }^2}{\langle {\overline{b}}^2{P}_{\mathrm{m}}(k){\left(\frac{1}{\overline{n}}\right)}^2\rangle }_i{W}_b(k){\delta }_{\mathit{ab}},\hfill \end{array}$$(19)

where V_i is the volume of the light-cone redshift slice. The model presented in Euclid Collaboration (2024) also contains nuisance parameters to correct inaccuracies due to modeling approximations. These corrections are only required at medium and high redshift (z ≳ 0.5), however, while at low redshift the model predicts the clustering covariance within an accuracy of 10%. Thus, in this work we do not need to fit any additional parameter for the clustering covariance.

2.4. Likelihood function

We considered the three statistics as independent. Following C19, we assumed that cluster counts and weak-lensing data are not correlated because the dominant systematic errors in the latter (shape noise, multiplicative shear, and photo-z bias) do not affect the number counts. In support of this, C19 also showed that the cross-correlation between the number counts and the mean weak-lensing cluster mass obtained from jackknife estimates is consistent with zero. For the same reason, we assumed that weak lensing is not correlated with clustering either. Moreover, the scales considered for the two observables in this work do not overlap. Finally, we show in Appendix A by means of 1000 mock catalogs that cluster counts and clustering present a negligible cross correlation for this type of survey.

We thus considered the total likelihood as the product of three independent Gaussian likelihood functions, each taking the form

$$\begin{array}{c}\hfill \mathcal{L}(\mathbf{d}|\mathbf{m}(\mathit{\theta }),C)=\frac{exp\{-\frac{1}{2}{[\mathbf{d}-\mathbf{m}(\mathit{\theta })]}^T{C}^{-1}[\mathbf{d}-\mathbf{m}(\mathit{\theta })]\}}{\sqrt{{(2\pi )}^N|C|}},\end{array}$$(20)

where d is the data vector, m is the predicted quantity depending on cosmological and nuisance parameters θ, and C is the covariance matrix. Following Euclid Collaboration (2021, 2024), we assumed for cluster counts and clustering that the covariance depends on cosmology, computed through the analytical models presented in the two works, respectively, as described in the previous sections. For weak-lensing log-masses, we assumed a Gaussian likelihood with fixed-cosmology covariance that accounts for shared multiplicative shear and photo-z biases, blended sources, and cluster triaxiality and projection effects. This is a common choice in cluster cosmology studies (Bocquet et al. 2019; Costanzi et al. 2019; Dark Energy Survey Collaboration 2020; Lesci et al. 2022a; Sunayama et al. 2023), as the dominant sources of uncertainties in the lensing-mass estimates, that is, shear and photo-z noise, are independent of cosmology, while other source of systematic error such as line-of-sight projections are weakly dependent on cosmology (e.g., Gruen et al. 2015; McClintock et al. 2019). Moreover, C19 explicitly tested by means of mock catalogs that the contribution of sample variance to the total error budget is negligible.

3. Data

The data set on which our analysis is based, as described in C19, is composed of 6964 photometrically selected galaxy clusters identified in the SDSS DR8, covering approximately 10 000 deg². The cluster identification process relies on the redMaPPer cluster-finding algorithm (Rykoff et al. 2014), which models the red sequence of galaxies and uses a probabilistic richness estimation to identify clusters. The algorithm iteratively refines its model and employs percolation to connect galaxies into clusters. It estimates the purity and completeness of the identified clusters and produces a catalog with cluster properties. The photometric redshift uncertainties for clusters are typically σ_z/(1 + z)≈0.01, and the analysis is confined to the redshift range z ∈ [0.1, 0.3] to ensure the accuracy of our photometric measurements and maintain a well-defined volume-limited catalog. Only clusters with a richness λ ≥ 20 were considered, ensuring that 99% of the redMaPPer galaxy clusters are unambiguously mapped to individual dark matter halos.

The observed number counts and weak-lensing mass measurements were described in C19. The number counts are measured in five richness bins λ^ob = {20, 27.9, 37.6, 50.3, 69.3, 140} and a single redshift bin z^ob ∈ [0.1, 0.3]. The weak-lensing masses were measured within the same intervals, from the shear catalog presented in Reyes et al. (2012), which includes about 39 000 000 galaxies over 9000 deg² of the SDSS footprint. The mass estimates, obtained by fitting an NFW profile (Navarro et al. 1997) to the stacked lensing profile of clusters in the radial range from 0.3 h⁻¹ Mpc to 3 h⁻¹ Mpc, are slightly updated from those presented in Simet et al. (2017). The measurements account for a broad range of systematic uncertainties, including shear calibration and photo-z biases, dilution by member galaxies, source obscuration, magnification bias, incorrect assumptions about the cluster mass profile, cluster miscentering, halo triaxiality, and projection effects. Because the recovered weak-lensing mass may not be identical to the mean mass of the clusters predicted by Eq. (8), the relation between the recovered weak-lensing mass and the mean mass in a bin was calibrated through the use of simulations due to the nonlinear relation between the stacked density profile ΔΣ and the mass M. The relation was directly applied to the observed data vector (see C19 for details).

For a flat ΛCDM cosmology, the cosmological dependence on Ω_m of the weak-lensing mass estimates can be approximated following C19 by the linear relation in log-space,

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

where the slopes derived in each bin from fitting this equation to the data are listed in Table 1 of C19. They were used in our cosmological analysis to rescale ${\widehat{M}}_{\mathrm{WL}}$ at each step of the MCMC.

To measure the two-point correlation function, we considered 30 log-spaced separation bins in the range r = 20 –130 h⁻¹ Mpc. This interval includes linear scales, where the bias is almost constant (Manera et al. 2010), plus the BAO peak, whose position and amplitude are sensitive to the density parameters and h. Because the available statistics are relatively low, we only considered a single richness threshold ${\lambda }_{\text{th}}^{\text{ob}}\ge 20$ instead of binning over richness. We performed the measurement by applying the Landy & Szalay (1993) estimator

$$\begin{array}{c}\hfill {\widehat{\xi }}_{\mathrm{h}}^{\mathit{aj}}=\frac{{\mathrm{DD}}_{\mathit{aj}}-2{\mathrm{DR}}_{\mathit{aj}}+{\mathrm{RR}}_{\mathit{aj}}}{{\mathrm{RR}}_{\mathit{aj}}},\end{array}$$(22)

where DD_aj,  DR_aj, and RR_aj are the number of pairs in the data-data, data-random, and random-random catalogs within the ath separation bin and jth redshift bin, normalized for the number of objects in the data and random catalogs, N_R and N_D (Kerscher et al. 2000). The measurement process was performed with the CosmoBolognaLib package (Marulli et al. 2016). Differently from number counts and weak-lensing masses, which were extracted from v5.10 of the SDSS redMaPPer cluster catalog, we measured the cluster clustering from the v6.3 catalog, for which the cluster random catalog is available. We verified that it has no significant impact on the results. The random catalog was generated by requiring that f_mask < 0.2 and λ/S > 20, where f_mask is the local mask fraction, and S is a scale factor, as defined in Rykoff et al. (2014). The random catalog is 120 times denser than the data catalog.

4. Results

Assuming a flat ΛCDM cosmological model with massive (degenerate) neutrinos, we constrained the six cosmological parameters Ω_m, ln(10¹⁰A_s), n_s, h, Ω_bh², and Ω_νh², plus the four mass-observable relation parameters of Eqs. (4) and (5), namely M_min, M₁, α, and σ_intr. In addition, following C19, we added the two parameters s and q to characterize the systematic uncertainty in the halo-mass function fit from simulations, such that

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}M}=\frac{\mathrm{d}n}{\mathrm{d}M}(s(M/M^{\ast })+q),\end{array}$$(23)

where (M^*) = 13.8 h⁻¹ M_⊙, and s = 1, q = 0 are the reference values. We assumed the Euclid Collaboration (2023) model for the halo-mass function (see Appendix B for details) and the Tinker et al. (2010) model for the halo bias. We modeled the effect of neutrinos following the cold dark matter prescription proposed by Costanzi et al. (2013), where the cold dark matter power spectrum was used instead of the total matter power spectrum.

We adopted the same cosmological prior as in C19, which was based on Dark Energy Survey Collaboration (2018). In brief, we set uninformative priors on parameters to which our data set was highly sensitive (Ω_m, σ₈, and scaling relation parameters), and informative priors on the others. C19 explicitly tested that the sizes of the adopted parameter priors do not affect the final posteriors. The priors we adopted for the parameter inference are listed in Table 1, and the sampler employed for the analysis is the Python wrapper for the nested sampler PolyChord (Handley et al. 2015).

4.1. Cosmological constraints

In the following, we present cosmological constraints derived by carrying out four analyses. We tested all possible combinations of the three cosmological probes introduced in Sect. 2:

• i)

  number counts and weak-lensing masses (NC+M_WL), which represent the standard analysis performed by C19;

• ii)

  number counts and clustering (NC+CL);

• iii)

  clustering and weak-lensing masses (CL+M_WL);

• iv)

  number counts, clustering, and weak-lensing masses (NC+CL+M_WL).

In Appendix B we compare our NC+M_WL analysis with the results of C19. The comparison demonstrates that the two analyses are fully consistent, and they provide some preliminary tests of the halo-mass function and halo bias models.

In Table 2 we report the best-fit values with 1σ uncertainty obtained by the four analyses, and in Fig. 1 we compare the related posterior distributions on the cosmological parameters Ω_m and σ₈ (which is a derived parameter), and of the mass-observable relation parameters M₁, α, σ_intr. The results were obtained by marginalizing over all the other cosmological and mass-richness parameters, and we only report results for parameters whose posterior is not dominated by a prior. The figure shows that the different degeneracies between parameters are constrained by different combinations of the three observables. In particular, cluster clustering is effective in constraining the matter density parameter Ω_m, breaking the Ω_m − σ₈ degeneracy proper of cluster counts, and shifting its constraints towards higher values. The main reason for this is the sensitivity of the power spectrum to Ω_m, which determines its shape. More specifically, this information is mainly extracted from the BAO scales: We verified, as shown in Appendix C, that restricting the clustering analysis to a separation range of r ∈ [20 − 60] h⁻¹ Mpc significantly broadens the final posteriors, while considering the range r ∈ [60 − 130] h⁻¹ Mpc allows us to recover almost all the information. No additional information is gained by extending the clustering analysis to r = 200 h⁻¹ Mpc.

Fig. 1. Contour plots at 68 and 95% of confidence level for different observable combinations. The number counts and weak-lensing masses are plotted in blue, the number counts and clustering are shown in green, the clustering and weak-lensing masses are shown in orange, and all the three probes are shown in red. The posteriors are obtained by marginalizing over the other parameters listed in Table 1; σ8 is a derived parameter. We only show the parameter posteriors that are not dominated by a prior.

Compared to other data combinations, the joint analysis NC+CL (green contours) has the least cosmological constraining power: The missing ability to constrain the slope of the mass-observable relation translates into almost uninformative posteriors on σ₈. The amplitude of the posterior on this parameter mainly depends on the amplitude of the prior on A_s. On the other hand, cluster clustering is extremely effective in constraining cosmology when it is combined with weak-lensing mass information (orange contours). In particular, the Ω_m and σ₈ contours shrink by 22% and 35% compared to the NC+M_WL analysis. Even though the posteriors in the scaling relation parameters are somewhat weaker than for the NC+M_WL combination, the break in the σ₈ − α degeneracy allows this data combination to provide the tightest constraints in the Ω_m − σ₈ plane (from the joint analysis of two observables). It should be noted that part of this constraining power comes from the clustering covariance matrix, which depends on shot-noise and thus contains information on the integrated halo-mass function (see Euclid Collaboration 2024, and discussion in Sect. 4.3).

As expected by the presence of different degeneracies and parameter dependences, the combination of the three observables (red contours) provides the tightest posteriors in whole parameter space. We obtain Ω_m = 0.28 ± 0.03 and σ₈ = 0.82 ± 0.05. Differently from the CL+M_WL case, which shows similar cosmological constraints, the three-probe combination also provides improved constraints on the scaling relation parameters. Because of the stacked approach used to estimate the mean cluster masses in the bins, the only scaling relation parameter that remains almost unconstrained is that describing the intrinsic scatter, as expected.

In addition to breaking degeneracies when they are combined with the other statistics, cluster clustering also helps us to constrain the Hubble parameter, as shown in Fig. 2. As indicated by the gray dashed contours, most of the information on the h parameter is provided by the clustering alone, whose dependence on the expansion rate of the Universe is captured by the shape of the two-point correlation function, in particular, around the BAO scales⁶. The best-fit value from the NC+M_WL+CL analysis is $h=0.{64}_{-0.04}^{+0.03}$, where the difference is only 1.90σ to the last Planck results (H₀ = 67.4 ± 0.5 km s⁻¹ Mpc⁻¹, Planck Collaboration VI 2020). Our constraints are opposite to those that were required to solve the tension with local probes, which prefer a higher value of H₀ (e.g., H₀ = (73.3 ± 1.1) km s⁻¹ Mpc⁻¹, Riess et al. 2022).

Fig. 2. Contour plots at the 68 and 95% confidence level in the Ωm − h plane for different observable combinations.

On the other hand, we verified that the cluster clustering does not appear to carry additional information on the total neutrino mass, and even for the full data combination, we retrieved posteriors on Ω_ν consistent with the prior.

4.2. Goodness of the fit and internal tensions

Figure 3 shows the measured number counts (left column), two-point correlation function (middle column), and mean cluster masses (right column) compared to the corresponding predictions computed by propagating the posteriors of Fig. 1 for the four different combinations of probes. For each data combination we considered, we obtained a good fit to data. In particular, the reduced χ² for our best-fit models are 1.91, 0.88, 0.94, and 1.04 for the NC+M_WL, NC+CL, CL+M_WL, and NC+M_WL+CL joint analyses, respectively. When only two data sets were combined, we derivde predictions for the third observable (first three rows of Fig. 1) to highlight possible internal tensions. For all these joint analyses, the predictions for the observable that is excluded from the fit are consistent within 2σ with the data. The worst agreement with the data is obtained for the mean cluster masses predicted by NC+CL, and it is driven by the broad posterior on σ₈.

Fig. 3. Observed number counts (left), two-point correlation function (middle), and weak-lensing masses (right), compared to the corresponding predicted quantities evaluated at the best-fit cosmology of the contours in Fig. 1 (same color code). The black dots are the observed points, and the colored lines and shaded areas represent the mean and standard deviation of each statistics, respectively, evaluated at 100 cosmologies that were randomly extracted from the posterior distribution, convoluted with the covariance matrix through a multivariate Gaussian distribution.

To quantitatively assess the level of tension among the data sets, we computed the posterior agreement⁷ (Bocquet et al. 2019), which quantifies whether the difference between two posterior distributions is consistent with zero. We computed the agreement for the set of parameters {Ω_m, σ₈, α, M₁}, and we compared the posteriors between pairs of probes. The results are reported in Table 3. The p-value indicates the probability of obtaining a difference between paired samples as extreme as or more extreme than the observed difference, under the assumption that there is no difference; a high (close to unity) p-value indicates a good agreement of the data sets. The worst level of agreement is given by the comparison of NC+M_WL and NC+CL, which can be interpreted as a small tension between clustering and weak-lensing masses. Similarly, a difference of about 1σ is obtained between number counts and weak-lensing masses (comparison of NC+CL and M_WL+CL), while a better agreement is reached by counts and clustering (comparison of NC+M_WL and M_WL+CL). However, none of the data combinations we considered exhibits a statistically significant tension (i.e., > 3σ). This is also consistent with the rightmost panel in the second line of Fig. 1, which shows the largest tension between weak-lensing masses and the prediction from the combination NC+CL. While not statistically significant, these results might indicate some unmodeled systematic effects in the derivation of the mean weak-lensing masses.

One systematic that was not accounted for in the analysis of Simet et al. (2017) is the optical selection bias induced by the correlation between richness and lensing signal (see, e.g., Dark Energy Survey Collaboration 2020; Sunayama et al. 2020; Wu et al. 2022; Salcedo et al. 2023). We explicitly tested this hypothesis and repeated the full joint analysis including a correction term in our mean mass model (Eqs. (8)–(10)) that accounted for this lensing-mass bias (see Appendix D for details). The analysis did not provide statistically significant evidence for the presence of biases in the mass estimates. On the other hand, this null result is in line with the fact that the full data combination does not show any internal tension.

4.3. Cosmology dependence of the covariance

In this section, we discuss the cosmology dependence of the covariance matrices. More specifically, the covariance models for the number counts and clustering were assumed to depend on the cosmology in the likelihood analysis, that is, they were recomputed at each step of the MCMC process. This should represent the proper way to perform the analysis because both the mean value and the covariance are needed to fully characterize the Gaussian model that we assumed to describe the data distribution (e.g., Eifler et al. 2009; Morrison & Schneider 2013; Krause & Eifler 2017; Blot et al. 2020). This correctly works for the number counts, where the assumption of Gaussian likelihood is well motivated (e.g., Payerne et al. 2023). A Gaussian likelihood is also a common choice for the analysis of the two-point correlation function, but there have been claims that the use of a cosmology-dependent covariance matrix may lead to an underestimation of the posteriors (Carron 2013). This approximation applied to the clustering of cluster statistic has been extensively tested in Euclid Collaboration (2024). The authors found that a Gaussian likelihood with a cosmology-dependent covariance matrix is statistically preferred over the fixed-cosmology case.

To evaluate the effective impact of the cosmology-dependent clustering covariance on the SDSS redMaPPer cluster catalog used in our analysis, we compare in Fig. 4 the cosmological posteriors obtained by using a cosmology-dependent covariance matrix and by fixing it at a single cosmology. The latter was chosen as the best-fit model obtained for the cosmology-dependent case. We performed this test for the three probe combinations that contain the cluster clustering. The most significant difference is observed for the CL+M_WL analysis, where the cosmology-dependent covariance matrix improves the constraints on σ₈ by approximately 100%. As already mentioned in Sect. 4.1, the additional information comes from the shot-noise term proportional to the inverse of the mean density of halos, that is, the integrated halo-mass function. On the other hand, the posterior distributions of the other parameters remain mostly unaffected by the cosmology-dependent covariance matrix. This suggests that the integrated halo-mass function only provides additional information on σ₈, while the dependence of the halo-mass function on the other parameters is smoothed out by integration over the entire mass spectrum.

Fig. 4. Comparison of Ωm − σ8 posteriors obtained using a cosmology-dependent (filled contours) and fixed-cosmology (empty contours) clustering covariance. The different colors represent different probe combinations (color-coded as Fig. 1). The fixed-cosmology covariance is computed at the best-fit parameters of the corresponding cosmology-dependent case.

Conversely, for the combination NC+CL, the two posteriors look very similar. This is expected because the cosmological information carried by the shot-noise term through the halo-mass function is already exhausted by the number count data, and thus the contribution of the cosmology-dependent covariance matrix becomes negligible. In this context, it might be argued that the information of the two observables is counted twice: As already stated, the clustering covariance includes a shot-noise contribution that contains part of the information of the number counts. Vice versa, the number count covariance contains part of the clustering information through the sample variance term. However, these terms do not directly affect the mean values, but contribute to modulate the fluctuations around the mean, thus characterizing the covariance matrices of the two statistics. It is noteworthy that if there were a correlation that is erroneously overlooked, the halo-mass function information would be counted twice and the cosmology-dependent contours of the NC+CL case (dark green lines) should tighten with respect to to the fixed-cosmology contours (light green contours), which is not the case. For the same reason, there is no difference between the cosmology-dependent and fixed-cosmology cases when all three probes are combined.

5. Discussion and conclusions

In this work, we have extended the analysis described in C19 by adding number counts and weak-lensing mass estimates to the cosmological information contained in the two-point correlation function of galaxy clusters. The twofold purpose of our analysis is to determine whether the inclusion of cluster clustering can help us, on the one hand, to improve cosmological constraints, and, on the other hand, to highlight possible tensions between the data sets. For this purpose, we analyzed the redMaPPer cluster catalog from the SDSS DR8 (Aihara et al. 2011), which contains 6964 clusters with richness λ^ob ≥ 20 and redshift z^ob ∈ [0.1, 0.3]. Assuming a flat ΛCDM model with massive neutrinos, we constrained the cosmological and mass-observable relation parameters using different combinations of the probes: cluster number counts, cluster clustering, and weak-lensing masses.

The combination of all the three probes has proved to be highly constraining for the parameters defining both the underlying cosmological model and the mass-richness relation. The inclusion of cluster clustering helps us to better constrain Ω_m and shifts its posterior toward higher values. Moreover, it allows us to obtain a precise measurement of the Hubble parameter ($h=0.{64}_{-0.04}^{+0.03}$), which is independent of and competitive with other probes. Instead, cluster clustering does not prove useful in constraining the neutrino mass any more than does the combination NC+M_WL. In all the cases, the posteriors on the neutrino mass are consistent with its prior.

Cluster clustering combined with number counts alone does not carry enough information on the scaling relation to compensate for the lack of dedicated mass-calibration data, and this comes at the expense of very weak constraints on σ₈. Instead, the different parameter degeneracies of clustering and weak-lensing masses allows us to combine them to constrain cosmology better than the standard combination of number counts and weak-lensing masses. It is interesting to note that while the amplitude of Ω_m − σ₈ contours is comparable with the three-probe case, the scaling relation parameters are much less constrained when number counts are not included in the fit. This indicates that the combination of clustering and weak-lensing masses without the full information of the number counts can tightly constrain cosmology even with a relatively limited knowledge of the scaling relation compared to the standard NC+M_WL analysis. However, even when number counts are not included in the analysis, part of the halo-mass function information is contained in the cosmology-dependent covariance matrix entering in the two-point correlation function likelihood. We showed that this information is useful to obtain better parameter constraints when the number counts are not directly included in the inference process.

In Fig. 5 we compare our contours with constraints obtained from other probes and from other cluster surveys. In the top panel, we consider the outcomes of the following analyses: Bocquet et al. (2019) analyzed the cluster count data from the 2500 deg² SPT-SZ catalog in combination with high-quality X-ray and lensing follow-up data for 121 systems; Chiu et al. (2023) examined cluster abundances from the X-ray sample identified in the eROSITA Final Equatorial Depth Survey (eFEDS) in combination with Hyper Supreme Camera weak-lensing data; Dark Energy Survey Collaboration (2022) combined cosmic shear, galaxy clustering, and galaxy-galaxy lensing from the third-year release of DES data (DES-Y3); and the baseline analysis of the CMB anisotropies with free neutrino mass conducted by Planck Collaboration VI (2020). In the bottom panel, we compare our results of the first-year release of the DES cluster sample data: To et al. (2021a) combined DES-Y1 cluster number counts and galaxy, shear, and cluster auto- and cross-correlation functions; Costanzi (2021) combined the DES-Y1 cluster count data with SPT-SZ multiwavelength data; Dark Energy Survey Collaboration (2020) combined DES-Y1 number counts and weak-lensing mass estimates. Our NC+M_WL+CL results in the σ₈ − Ω_m plane agree with all the other analyses we considered within one sigma except for the Ω_m value preferred by Planck and Dark Energy Survey Collaboration (2020), which are different at 1.22σ and 2.21σ from our results. Our results also show a small tension with those of Dark Energy Survey Collaboration (2022), corresponding to 1.35σ on Ω_m and 1.30σ on σ₈. The constraining power is similar to those of other recent cluster abundance studies. The comparison with the DES-Y1 NC+4x2pt results is particularly interesting because the two cluster samples have similar sizes, and the observables considered in this analysis mirror those analyzed by To et al. (2021a; in addition to the galaxy clustering and galaxy-cluster cross-correlations). It is worth mentioning that our contours fully match the DES results (at 0.45σ and 0.47σ level for Ω_m and σ₈, respectively), but the amplitude is more constrained. A possible explanation for the larger contours found in To et al. (2021a) are the different scales considered in the two-point correlation functions: The analysis is limited to a scale of 100 h⁻¹ Mpc, so that part of the information carried by the BAO peak is lost. In this work, the clustering was measured up to 130 h⁻¹ Mpc, allowing us to fully exploit the BAO feature (see Appendix C). At the lower end, they did not use data below 8 h⁻¹ Mpc and thus excluded the one-halo term information from the lensing analysis and significantly reduced the capability of constraining the cluster masses. Moreover, To et al. (2021a) considered a selection bias term, which modulates the amplitude of the cluster auto- and cross-correlation functions and further reduces the constraining power of the sample.

Fig. 5. Comparison of Ωm − σ8 posteriors obtained by different surveys. Top panel: 2500 deg2 SPT-SZ cluster survey (Bocquet et al. 2019, blue contours), eROSITA (eFEDS) cluster survey (Chiu et al. 2023, yellow contours), DES-Y3 NC+3x2pt (Dark Energy Survey Collaboration 2022, green contours), and Planck TT,TE,EE+lowE with free mν (Planck Collaboration VI 2020, black countours), compared to our results in the NC+MWL+CL (red contours) analyses. Bottom panel: DES-Y1 number counts and weak-lensing mass estimates (Dark Energy Survey Collaboration 2020, cyan contours), DES-Y1 cluster count data and SPT-SZ follow-up data (Costanzi 2021, orange contours), and DES-Y1 NC+4x2pt (To et al. 2021a, dark green contours), compared to our results from the NC+MWL (dotted gray contours) and NC+MWL+CL (red contours) results.

In Fig. 6 we compare the mean richness-mass relation derived from our NC+M_WL+CL analysis with those obtained from the analyses of the DES-Y1 cluster sample considered above. For the comparison, all the DES-Y1 relations were evolved to the mean redshift of the SDSS sample, z = 0.22, and were rescaled by the factor 0.93 to correct for the systematic richness offset between SDSS and DES-Y1 catalogs (see McClintock et al. 2019). All the scaling relations agree well over the relevant mass range, M ≳ 10^13.5 h⁻¹ M_⊙ except for the Dark Energy Survey Collaboration (2020) results, which exhibit a shallower slope and deviate from the other relations by 1σ at the low-mass end. As discussed in Dark Energy Survey Collaboration (2020), this difference is likely due to a flawed modeling of the stacked weak-lensing signal, and it drives the tension with other cosmological probes. It is interesting to note that despite the similar technique and modeling adopted to estimate the weak-lensing masses in this work and in Dark Energy Survey Collaboration (2020), the SDSS result does not seem to be affected by the same systematics. A possible explanation is the different scale cut adopted in the lensing analyses: Simet et al. (2017) considered only scales below R ≃ 3 h⁻¹ Mpc, while McClintock et al. (2019) included scales up to 30 h⁻¹ Mpc and also modeled the two-halo term. These large scales are most affected by the optical selection bias discussed in Sect. 4.2 (see e.g., Sunayama et al. 2020; Wu et al. 2022), which represents one of the main systematics in Dark Energy Survey Collaboration (2020). Conversely, To et al. (2021a) excluded the one-halo term from the analysis by considering R > 8 h⁻¹ Mpc, but included a selection bias correction term that was left free to vary in the analysis, which can explain the good agreement with our result. Finally, the weak-lensing follow-up data used in Costanzi (2021) to calibrate the scaling relations, being based on a SZ-selected cluster sample, are not affected by this selection bias.

Fig. 6. Comparison of the mean richness-mass relations derived in this work (red band) and those derived from different analyses of the DES-Y1 redMaPPer cluster sample: DES-Y1 number counts and weak-lensing mass estimates (Dark Energy Survey Collaboration 2020, cyan band), DES-Y1 cluster count data and SPT-SZ follow-up data (Costanzi 2021, orange-dotted band), and DES-Y1 cluster abundance with galaxy, lensing, and cluster auto- and cross-correlation function (To et al. 2021a, black-hatched band). The band width corresponds to a 1σ uncertainty of the mean relation. Lower panel: ratio of the DES-Y1 scaling relations with respect to the SDSS result.

Finally, it is worth mentioning the comparison with Sunayama et al. (2023), who also analyzed the SDSS redMaPPer cluster catalog in combination with the Hyper-Suprime Cam (HSC) Year3 shape catalog by combining cluster abundances, weak-lensing masses, and projected cluster clustering. We obtain consistent results, but their posterior amplitude is slightly larger. Again, this is mainly ascribed to the different scales used for cluster clustering: In their analysis, the clustering is measured up to 50 h⁻¹ Mpc, thus excluding the BAO peak that carries most of the information. The broadening of the contours generated by examining the range of separation r ∈ [20 − 60] h⁻¹ Mpc shown in Fig. C.1 (blue contours) is wide enough to be consistent with the amplitude of their posteriors. Moreover, they included a selection bias term that modulates both the lensing and the two-point correlation function amplitude.

In summary, in the analysis presented here, we showed some of the potential of cluster clustering, which has proven to be useful for constraining cosmology, although the catalog we analyzed still has rather limited statistics and a narrow redshift coverage. To better exploit the clustering information, a full redshift- and richness-dependent analysis should be performed. Moreover, the inclusion of redshift-space distortions in the modeling of the two-point function may further increase the constraining power. Finally, exploring cluster clustering analysis within the context of nonstandard cosmological models is an interesting avenue of research. This approach allowed us to investigate the impact of phenomena such as dark energy or modified gravity on the formation and evolution of the large-scale structure of the Universe. These nonstandard models have the potential to induce substantial alterations in the clustering patterns of cosmic structures, including that of galaxy clusters. Notably, in these nonstandard scenarios, the halo bias may acquire a scale-dependent character, offering a valuable way for probing and constraining the underlying cosmological framework through a cluster clustering analysis.

We expect that future large survey catalogs, such as or LSST, will make it possible to measure cluster clustering with growing accuracy, making its combined analysis with other cosmological probes increasingly powerful and competitive.

Acknowledgments

A.S., T.C., and S.B. are supported by the INFN INDARK PD51 grant. A.F. acknowledges support from Brookhaven National Laboratory. A.S. is also supported by the ERC ‘ClustersXCosmo’ grant agreement 716762. T.C. and A.S. are also supported by the FARE MIUR grant ‘ClustersXEuclid’ R165SBKTMA.

References

Appendix A: Cluster counts and clustering cross-correlation matrix

In this appendix, we discuss the cross-correlation between the number counts and clustering. In principle, these observables are correlated by density fluctuations on scales larger than the survey size, that is, the supersample covariance (see, e.g., Takada & Hu 2013; Krause & Eifler 2017). The effect of the supersample covariance only depends on the survey size and should become negligible for large volumes. Moreover, it mostly impacts small, nonlinear scales.

To assess the entity of this correlation in a survey with properties similar to those analyzed in this work, we used a set of 1000 light cones simulated with the PINOCCHIO code (Monaco et al. 2002; Munari et al. 2017), which are described in detail in Euclid Collaboration (2021). The light cones covered an area of about 13 000 deg², which is comparable with the area of the cluster catalog we analyzed in this work, and we considered the redshift range z ∈ [0.1, 0.3]. We measured number counts and two-point correlation function from each of the 1000 catalogs and used them to compute a numerical covariance matrix, whose normalized version is shown in Fig. A.1. The block-diagonal elements represent the cluster counts (upper left) and cluster clustering (lower right) covariances, while the off-diagonal elements are the cross-covariance between the two statistics. The cross correlation is entirely dominated by noise and consistent with zero signal.

Fig. A.1. Normalized covariance matrix between cluster counts and cluster clustering. The block-diagonal terms represent the autocorrelation matrix for the number counts (upper left) and clustering (lower right), and the off-diagonal blocks are the cross-correlation matrix.

This result suggests that the two observables can be considered to be independent. A similar outcome was obtained in Sunayama et al. (2023).

Appendix B: Validation tests

In this appendix, we present the results from tests with which we validated our analyses.

In Fig. B.1 we compare our results for the NC+M_WL combination with the same posterior distribution as C19. By comparing our results (solid black contours) with the C19 results (shaded blue area), we note a good agreement in general, with only a small difference in the posterior amplitudes. This difference can be ascribed to numerical issues that are due to the use of a different code. Although the likelihood and sampler we adopted are the same, it is unlikely that two different codes achieve the same numerical accuracy. This may produce sizable differences in the predictions of the observables and, consequently, in the cosmological constraints. However, we reiterate that the two results from our analysis and from that of C19 are fully consistent with each other. This consistency check ensures the absence of additional systematic effects in our analysis with respect to C19 and confirms the good agreement of the two analyses.

Fig. B.1. Comparison between the results of C19 and this work. The dashed gray lines are the best-fit parameters from C19.

We also assessed the impact of different halo-mass function models and of their interaction with the halo bias model. The halo bias was derived from the halo-mass function through the peak-background split formalism (Cole & Kaiser 1989; Mo & White 1996), and it is thus related to the mass function model from which it is obtained. We tested three parameterizations for the halo-mass function, that is, Tinker et al. (2008), Tinker et al. (2010), and Euclid Collaboration (2023), while we only assumed a single model for the halo bias, the model from Tinker et al. (2010).

In Fig. B.2 we show the posteriors resulting from the use of different halo-mass function parameterizations. The orange contours refer to the model of Tinker et al. (2008), the blue contours show the model of Tinker et al. (2010), and the red contours show the model of Euclid Collaboration (2023). We show here the constraints from the M_WL+CL analysis to better estimate the impact of using the same halo bias in combination with different halo-mass function models. The results from the other combinations of observables are consistent with those presented here. By comparing the results, we note that the halo-mass function model does not have any significant impact on the parameter constraints at this level of statistics. Neither does the use of a bias model that is not associated with the halo-mass function produce measurable effects on the contours. From this result, we conclude that the choice of the halo-mass function can be made without any particular constraint. Therefore, we chose to use the model of Euclid Collaboration (2023), which has been demonstrated to be accurate at the subpercent level when tested on numerical simulations.

Fig. B.2. Comparison between different halo-mass function models from the MWL+CL analysis. Tinker et al. (2008) is plotted in orange, Tinker et al. (2010) is shown in blue, and Euclid Collaboration (2023) is shown in red.

Appendix C: Clustering scales

We show in Fig. C.1 the cosmological constraints from the NC+M_WL+CL analysis that we obtained by varying the clustering scales. The standard separation range adopted throughout this paper is given by r ∈ [20 − 130] h⁻¹ Mpc (empty black contours). We split this range in order to quantify the impact of the BAO scales on the final parameter constraints. In this way, we obtained a low-scale interval r ∈ [20 − 60] h⁻¹ Mpc (filled blue contours) and a high-scales interval r ∈ [60 − 130] h⁻¹ Mpc (filled orange contours). The comparison highlights that the low-scale part of the (linear) 2PCF has a far lower predictive power than that measured around the BAO peak, which carries most of the information. A similar outcome is found by expanding the range beyond the BAO scales, that is, r ∈ [20 − 200] h⁻¹ Mpc: The contours remain almost unchanged, indicating that no further information can be extracted from the 3D cluster clustering. The latter test also assures us that our results do not include border effects when the 2PCF is modeled because no systematic errors due to the assumption of spherical shells (see Eq. 14) are added when the separation is increased to 200 h⁻¹ Mpc.

Fig. C.1. Contour plots at the 68 and 95 percent of confidence level, corresponding to different separation ranges in the cluster clustering modeling. In all the cases, the contours correspond to the NC+MWL+CL combination.

Appendix D: Richness and weak-lensing correlation

We attempted to increase our information of the clustering and to constrain a possible mass bias induced by the correlation of the richness and lensing signal (see Sect. 4.2). In particular, we assumed that the mean weak-lensing mass and richness follow a bivariate log-normal distribution. With this assumption, we corrected the mean-mass predictions by adding to Eq. (10) the bias term

$$\begin{array}{c}\hfill \rho {\sigma }_{M_{\mathrm{WL}}}\frac{ln{\lambda }^{\mathrm{tr}}-⟨ln{\lambda }^{\mathrm{tr}}|M⟩}{{\sigma }_{ln{\lambda }^{\mathrm{tr}}}},\end{array}$$(D.1)

where ρ is the correlation coefficient between the richness and weak-lensing mass, σ_MWL is the intrinsic scatter at a fixed mass of the latter, and ⟨lnλ^tr|M⟩ and σ_lnλ^tr are given by Eqs. (4) and (5). In addition to assuming a log-normal distribution for the richness-mass relation P(λ^tr|M, z^tr), we neglected for this exercise the observational scatter in the richness mass relation, P(λ^ob|λ^tr, z^tr) in Eq. (3), and instead assumed λ^ob = λ^tr. We verified that these two approximations vary the posteriors very little. In addition to the parameters listed in Table 1, we set ρ and σ_MWL to be free parameters with priors given by ρ = [ − 1, +1] and σ_MWL = 𝒩(0.2, 0.1). The three-probe combination NC+CL+M_WL was examined. The analysis is not conclusive because the two parameters remain basically unconstrained. The lack of meaningful constraints derived on these two parameters is in line with the result that tensions disappear in the analysis based on the combination of the three probes.

All Tables

All Figures

Fig. 1. Contour plots at 68 and 95% of confidence level for different observable combinations. The number counts and weak-lensing masses are plotted in blue, the number counts and clustering are shown in green, the clustering and weak-lensing masses are shown in orange, and all the three probes are shown in red. The posteriors are obtained by marginalizing over the other parameters listed in Table 1; σ8 is a derived parameter. We only show the parameter posteriors that are not dominated by a prior.

In the text

	Fig. 2. Contour plots at the 68 and 95% confidence level in the Ωm − h plane for different observable combinations.
In the text

Fig. 3. Observed number counts (left), two-point correlation function (middle), and weak-lensing masses (right), compared to the corresponding predicted quantities evaluated at the best-fit cosmology of the contours in Fig. 1 (same color code). The black dots are the observed points, and the colored lines and shaded areas represent the mean and standard deviation of each statistics, respectively, evaluated at 100 cosmologies that were randomly extracted from the posterior distribution, convoluted with the covariance matrix through a multivariate Gaussian distribution.

In the text

	Fig. 4. Comparison of Ωm − σ8 posteriors obtained using a cosmology-dependent (filled contours) and fixed-cosmology (empty contours) clustering covariance. The different colors represent different probe combinations (color-coded as Fig. 1). The fixed-cosmology covariance is computed at the best-fit parameters of the corresponding cosmology-dependent case.
In the text

Fig. 5. Comparison of Ωm − σ8 posteriors obtained by different surveys. Top panel: 2500 deg2 SPT-SZ cluster survey (Bocquet et al. 2019, blue contours), eROSITA (eFEDS) cluster survey (Chiu et al. 2023, yellow contours), DES-Y3 NC+3x2pt (Dark Energy Survey Collaboration 2022, green contours), and Planck TT,TE,EE+lowE with free mν (Planck Collaboration VI 2020, black countours), compared to our results in the NC+MWL+CL (red contours) analyses. Bottom panel: DES-Y1 number counts and weak-lensing mass estimates (Dark Energy Survey Collaboration 2020, cyan contours), DES-Y1 cluster count data and SPT-SZ follow-up data (Costanzi 2021, orange contours), and DES-Y1 NC+4x2pt (To et al. 2021a, dark green contours), compared to our results from the NC+MWL (dotted gray contours) and NC+MWL+CL (red contours) results.

In the text

Fig. 6. Comparison of the mean richness-mass relations derived in this work (red band) and those derived from different analyses of the DES-Y1 redMaPPer cluster sample: DES-Y1 number counts and weak-lensing mass estimates (Dark Energy Survey Collaboration 2020, cyan band), DES-Y1 cluster count data and SPT-SZ follow-up data (Costanzi 2021, orange-dotted band), and DES-Y1 cluster abundance with galaxy, lensing, and cluster auto- and cross-correlation function (To et al. 2021a, black-hatched band). The band width corresponds to a 1σ uncertainty of the mean relation. Lower panel: ratio of the DES-Y1 scaling relations with respect to the SDSS result.

In the text

	Fig. A.1. Normalized covariance matrix between cluster counts and cluster clustering. The block-diagonal terms represent the autocorrelation matrix for the number counts (upper left) and clustering (lower right), and the off-diagonal blocks are the cross-correlation matrix.
In the text

	Fig. B.1. Comparison between the results of C19 and this work. The dashed gray lines are the best-fit parameters from C19.
In the text

	Fig. B.2. Comparison between different halo-mass function models from the MWL+CL analysis. Tinker et al. (2008) is plotted in orange, Tinker et al. (2010) is shown in blue, and Euclid Collaboration (2023) is shown in red.
In the text

	Fig. C.1. Contour plots at the 68 and 95 percent of confidence level, corresponding to different separation ranges in the cluster clustering modeling. In all the cases, the contours correspond to the NC+MWL+CL combination.
In the text