C³ Cluster Clustering Cosmology I. New Constraints on the Cosmic Growth Rate at z ∼ 0.3 from Redshift-space Clustering Anisotropies

The catalog analyzed in this work consists of optically selected clusters of galaxies that have been identified by Wen et al. (2012, hereafter WHL12) ¹⁰ from the SDSS-III, Data Release (DR) 8 (Aihara et al. 2011).

The WHL12 catalog lists 132,684 galaxy clusters on a sky area of ∼15,000 deg², spanning the redshift range of 0.05 < z < 0.8. The cluster identification is based on a friends-of-friends procedure (Huchra & Geller 1982). This approach has been already exploited to find groups and clusters using volume-limited spectroscopic samples of galaxies (see, e.g., Berlind et al. 2006; Tempel et al. 2014) at low redshifts (z < 0.2). The WHL12 cluster sample extends the technique on photometric redshift samples of galaxies, allowing the detection of galaxy overdensities around the BCGs at higher redshifts.

A candidate cluster is included in the catalog if N₂₀₀ ≥ 8 and ${R}_{{L}_{* }}\geqslant 12$, where N₂₀₀ is the number of member candidates within ${\tilde{r}}_{200}$, and ${R}_{{L}_{* }}$ is the optical richness defined as ${R}_{{L}_{* }}={\tilde{L}}_{200}/{L}_{* }$, where ${\tilde{L}}_{200}$ is the total r-band luminosity within an empirically determined radius ${\tilde{r}}_{200}$ and L_* is the evolved characteristic galaxy luminosity (Blanton et al. 2003). The subscript 200 denotes quantities measured in a sphere whose mean density is 200 times the critical density at the cluster redshift. The cluster photometric redshifts reported in the catalog are the median value of the photometric redshifts of the galaxy members. These selections have been applied to avoid contaminations by bright field galaxies with overestimated photometric redshifts.

The cluster masses are estimated from the weak-lensing cluster mass scaling relation calibrated in WHL12, with data from Wen et al. (2010):

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{M}_{200}}{{10}^{14}{M}_{\odot }}\right)=(-1.49\pm 0.05)+(1.17\pm 0.03)\mathrm{log}\left({R}_{{L}_{* }}\right).\end{eqnarray} \tag{ 1 }$$

To verify the robustness of our analysis, we consider also the scaling relation provided independently by Covone et al. (2014), finding consistent results.

According to WHL12, the catalog completeness, which is the fraction of the selected galaxy clusters over the full sample, is close to 1 in the redshift range of 0.1 < z < 0.42, for M₂₀₀ ≳ 2 × 10¹⁴ M_⊙, while the detection rate decreases down to ∼75%, including all clusters down to the minimum mass of the sample, M₂₀₀ = 6 × 10¹³ M_⊙. WHL12 also quantified the false cluster detections to be at the level of 6% for ${R}_{{L}_{* }}=12$, decreasing to <1% for cluster of richness ${R}_{{L}_{* }}\geqslant 23$. Possible effects on the properties of the WHL12 cluster sample caused by incompletenesses of SDSS-III's Baryon Oscillation Spectroscopic Survey (BOSS) galaxies at high stellar masses (Leauthaud et al. 2016; Saito et al. 2016) are neglected in this work. As we will detail in Section 4, our statistical analysis does not depend on the cluster catalog completeness.

An accurate and precise estimate of the redshift is crucial when reconstructing the statistical properties of the large-scale distribution of matter. Large redshift uncertainties, as in photometric redshift surveys, lead to severe distortion effects that reflect in the 2PCF measurement, complicating its analysis and cosmological interpretation (see, e.g., Marulli et al. 2012; Sereno et al. 2015; García-Farieta et al. 2020). In order to construct a spectroscopic cluster sample, we take advantage of the spectroscopic data from SDSS, focusing on the final spectroscopic DR12 from BOSS (Dawson et al. 2013; Alam et al. 2015, 2017), which is part of the SDSS-III program. This survey measured the spectra for millions of galaxies. We assign spectroscopic redshifts to WHL12 clusters by crossmatching with the spectroscopic galaxy sample. ¹¹ The total cluster catalog with spectroscopic information consists of 72,563 objects, spanning the redshift range 0 < z < 1. Following WHL12, we cut the sample in the redshift range 0.1 < z < 0.42 to minimize incompleteness uncertainties. The number of remaining clusters is 43,743, with a median redshift of z ∼ 0.3, covering an area of about 10,800 deg².

We make no distinction between galaxy clusters and BCGs in this analysis, since the coordinates of galaxy clusters are estimated as the coordinates of their BCGs, without further refinements. Clusters with no spectroscopic information for their BCGs are discarded, even if some of their member galaxies have a measured spectroscopic redshift. The rationale of this choice is to reduce contamination from nonlinear dynamics in hosting virialized dark matter haloes, thus minimizing the impact of theoretical uncertainties in the RSD modeling at small scales.

Figure 1 shows the angular distribution of the spectroscopic galaxy cluster catalog analyzed in this work, compared to the WHL12 photometric sample. The three north fields and the south one are shown with different colors. The early (E) north fields (E2, E3) have been included in DR12. They are characterized by a lower galaxy density and a different redshift distribution, relative to the north and south fields (Beutler et al. 2017), and will be treated differently when constructing the random catalog (see Section 3.1).

Zoom In Zoom Out Reset image size

Figure 2 compares the redshift and mass distributions of the spectroscopic cluster sample analyzed in this work to the ones of the WHL12 photometric and spectroscopic cluster catalogs, and of the WHL12 sample restricted to the BOSS area. The shape of the mass distribution of the selected cluster sample is overall consistent with theoretical ΛCDM predictions by Tinker et al. (2008). However, we do not attempt to exploit the cluster mass distribution in this work, to avoid systematics due to possible inaccurate knowledge of the sample selection function. The estimated masses are used instead to set a prior on the linear bias of the selected cluster sample.

Zoom In Zoom Out Reset image size

To estimate the three-dimensional 2PCF of a sample of extragalactic sources, a geometric selection function is needed. The clustering estimator adopted in this work (see Section 3.2) requires this function to be provided as a catalog of objects randomly distributed in the same area of the real catalog, and with the same selection along the line of sight. As we will explain in Section 4, the likelihood we will use to extract cosmological constraints from the analyzed clustering data set is independent of the catalog completeness, which is the fraction of selected galaxy clusters over the full sample (see Marulli et al. 2018, for a similar analysis). Thus, our final results depend only on the geometric selection that enters the clustering estimator.

We construct the random catalog following the same methodology used in galaxy clustering analyses, in particular, the one used to measure the 2PCF of BOSS galaxies. In the assumption that the angular and redshift distributions of the selected galaxy clusters are independent, we assign angular coordinates, i.e., R.A., decl., and redshifts in two separate steps. The R.A.-decl. coordinates are extracted with MANGLE (Swanson et al. 2008), using publicly available survey footprints. ¹² The redshifts are then sampled from the true redshift distribution of the catalog, smoothed with a Gaussian kernel of σ_z = 0.02 in order not to introduce spurious clustering along the line of sight. As we verified, the impact of this assumption is negligible. Due to the different density and redshift distributions in the north, north E2, north E3, and south fields, the extraction of random coordinates and redshifts is performed in each field separately. The final random catalog, which is obtained by adding the random catalogs in the four fields, is constructed to be 50 times larger than the selected cluster catalog in order to minimize the impact of the shot noise. The redshift distribution of the random objects normalized to the number of the spectroscopic cluster sample is shown in Figure 2.

3.1.1. Weights

We take into account the systematic uncertainties in the angular selection function due to the position-dependent completeness of the cluster sample. To do this, we cross correlate the angular cluster counts with the maps of observational systematics provided by Leistedt & Peiris (2014). ¹³ We find that the cluster counts are anticorrelated with stellar densities, and less strongly, with the r-band extinction. We correct for this effect by weighting the objects in the random sample accordingly.

Moreover, we weight the clusters in the catalog to account for the spectroscopic target selection, following the same weighting scheme adopted by Reid et al. (2016) and Ross et al. (2017):

$$\begin{eqnarray}&&w={w}_{\mathrm{see}}{w}_{\mathrm{star}}({w}_{\mathrm{cp}}+{w}_{\mathrm{noz}}-1),\end{eqnarray} \tag{ 2 }$$

which considers the impact of seeing (w_see), star contamination (w_star), and presence of a close target (w_cp), as well as spectroscopic measurement failures (w_noz).

The cosmological analysis performed in this work is based on the redshift-space clustering wedges of the 2PCF of the spectroscopic cluster catalog presented in Section 2.2. The clustering wedges have been introduced by Kazin et al. (2010) as a convenient projection statistic, similar to the clustering multipoles, to compress the anisotropic 2PCF signal. The main advantage of this approach is to reduce the dimension of the data set to be analyzed, and the associated covariance matrix.

To measure the three-dimensional 2PCF, we first convert the observed coordinates of the galaxy clusters (R.A., decl., redshift) into comoving Cartesian coordinates, assuming Planck 18 cosmology.

The comoving distances, d_c, are related to the cosmological redshifts, z, as follows:

$$\begin{eqnarray}&&{d}_{{\rm{c}}}(z)=c{\int }_{0}^{z}\displaystyle \frac{{\rm{d}}z^{\prime} }{H(z^{\prime} )},\end{eqnarray} \tag{ 3 }$$

where c is the speed of light, and H is the Hubble parameter, which in a flat ΛCDM model reads as

$$\begin{eqnarray}&&H(z)={H}_{0}{\left[{{\rm{\Omega }}}_{{\rm{M}}}{\left(1+z\right)}^{3}+(1-{{\rm{\Omega }}}_{{\rm{M}}})\right]}^{1/2}.\end{eqnarray} \tag{ 4 }$$

Neglecting redshift uncertainties and second-order corrections, the observed redshift, z_obs, is related to the cosmological redshift, z, as follows:

$$\begin{eqnarray}&&{z}_{\mathrm{obs}}=z+\displaystyle \frac{{v}_{\parallel }}{c}(1+z),\end{eqnarray} \tag{ 5 }$$

where v_∥ is the peculiar velocity along the line of sight. In this analysis, the impact of cluster redshift uncertainties on the 2PCF is minor, especially on large scales, as we consider only the spectroscopic cluster sample. Assuming that the cluster spectroscopic redshift uncertainties follow a Gaussian distribution (e.g., Sereno et al. 2015), their effects on the 2PCF are degenerate with those of small-scale peculiar random motions, and do not require any additional parameters to be modeled (see Section 4).

Since the peculiar velocities of the analyzed cluster sample are unknown, we estimate the comoving distances by substituting z with z_obs in Equation (3), thus introducing the so-called RSD. Hereafter, the redshift-space spatial coordinates are indicated with s.

The anisotropic 2PCF in redshift-space is computed with the Landy & Szalay (1993) estimator:

$$\begin{eqnarray}&&\hat{\xi }(s,\mu )=\frac{{N}_{\mathrm{RR}}}{{N}_{\mathrm{CC}}}\frac{\mathrm{CC}(s,\mu )}{\mathrm{RR}(s,\mu )}-2\frac{{N}_{\mathrm{RR}}}{{N}_{\mathrm{CR}}}\frac{\mathrm{CR}(s,\mu )}{\mathrm{RR}(s,\mu )}+1,\end{eqnarray} \tag{ 6 }$$

where μ is the cosine of the angle between the line of sight and the comoving separation s, CC(s, μ), RR(s, μ), and CR(s, μ) are the numbers of cluster-cluster, random-random, and cluster-random pairs in bins of s and μ, i.e. in s ± Δs and μ ± Δμ, N_C, and N_R are the total numbers of clusters and random objects, and N_CC = N_C(N_C − 1)/2, N_RR = N_R(N_R − 1)/2, and N_CR = N_C N_R are the total numbers of cluster-cluster, random-random, and cluster-random pairs, respectively. The Landy & Szalay (1993) estimator of the 2PCF is widely used as it provides the minimum variance when ∣ξ∣ ≪ 1, and it is unbiased in the limit of an infinitely large random sample (Keihänen et al. 2019). We estimate the comoving separation associated with each bin as the average cluster pair separation inside the bin (e.g., Zehavi et al. 2011).

Lastly, to efficiently compress the information contained in the clustering signal, we estimate the so-called wedges of the 2PCF (Kazin et al. 2012), that consist in the integrals of ξ(s, μ) over wide bins of μ:

$$\begin{eqnarray}&&{\xi }_{w}(s)\equiv \displaystyle \frac{1}{{\rm{\Delta }}\mu }{\int }_{{\mu }_{1}}^{{\mu }_{2}}d\mu \,\xi (s,\mu ),\end{eqnarray} \tag{ 7 }$$

where Δμ = μ₂ − μ₁ is the wedge width. Here, we set Δμ = 0.5, which leads to two clustering wedges, i.e., the transverse wedge, ${\xi }_{\perp }(s)\equiv {\xi }_{1/2}({\mu }_{\min }=0,s)$, and the radial wedge, ${\xi }_{\parallel }(s)\equiv {\xi }_{1/2}({\mu }_{\min }=0.5,s)$, computed in the ranges of 0 ≤ μ < 0.5 and 0.5 ≤ μ ≤ 1, respectively.

We model the redshift-space transverse and radial wedges of the 2PCF of our cluster catalog with an extended version of the Taruya et al. (2010) model, which includes the nonlinear biasing model by McDonald & Roy (2009). Following Beutler et al. (2014), we will refer to this as the extended Taruya, Nishimichi & Saito (eTNS) model.

The redshift-space power spectrum of galaxy clusters in the eTNS model is approximated as follows:

$$\begin{eqnarray}&&\begin{array}{l}{P}^{s}(k,\mu )=D(k,\mu ,f,{\sigma }_{v})\left[{P}_{{\rm{c}},\delta \delta }(k)+2f{\mu }^{2}{P}_{{\rm{c}},\delta \theta }(k)\right.\\ \left.+\,{f}^{2}{\mu }^{4}{P}_{\theta \theta }(k)+{b}_{1}^{3}{C}_{A}(k,\mu ,f,{b}_{1})+{b}_{1}^{4}{C}_{B}(k,\mu ,f,{b}_{1})\right],\end{array}\end{eqnarray} \tag{ 8 }$$

where f is the linear growth rate, b₁ is the linear bias, θ( k ) ≡ [−i k · v ( k )]/[af(a)H(a)] is the velocity divergence, P_{c,δ δ} (k) is the real-space density cluster power spectrum, P_{c,δ θ} (k) and P_{θ θ} (k) are the real-space density-velocity divergence cross spectrum and the real-space velocity divergence auto-spectrum of clusters, respectively, assuming no velocity bias, i.e., P_{c,θ θ} (k) = P_{θ θ} (k), D(k, f, μ, σ_v ) is a damping factor used to model the random peculiar motions at small scales, and C_A and C_B are two additional terms to correct for systematics at small scales. The cluster power spectra are computed with the nonlinear biasing model by McDonald & Roy (2009) as follows:

$$\begin{eqnarray}\begin{array}{rcl}{P}_{{\rm{c}},\delta \delta }(k) & = & {b}_{1}^{2}{P}_{\delta \delta }(k)+2{b}_{2}{b}_{1}{P}_{b2,\delta }(k)+2{b}_{{s}^{2}}{b}_{1}{P}_{{bs}2,\delta }(k)\\ & & +\,2{b}_{3\mathrm{nl}}{b}_{1}{\sigma }_{3}^{2}(k){P}_{{\rm{m}}}^{\mathrm{lin}}(k)+{b}_{2}^{2}{P}_{b22}(k)\\ & & +\,2{b}_{2}{b}_{{s}^{2}}{P}_{b2s2}(k)+{b}_{{s}^{2}}^{2}{P}_{{bs}22}(k)+N,\end{array}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\begin{array}{rcl}{P}_{{\rm{c}},\delta \theta }(k) & = & {b}_{1}{P}_{\delta \theta }(k)+{b}_{2}{P}_{b2,\theta }(k)+{b}_{{s}^{2}}{P}_{{bs}2,\theta }(k)\\ & & +\,{b}_{3\mathrm{nl}}{\sigma }_{3}^{2}(k){P}_{{\rm{m}}}^{\mathrm{lin}}(k),\end{array}\end{eqnarray} \tag{ 10 }$$

where ${P}_{{\rm{m}}}^{\mathrm{lin}}(k)$ is the linear power spectrum. The adopted biasing model has the following four parameters, besides the linear bias term, b₁: the second-order local and nonlocal bias parameters, b₂ and b_s ², the third-order nonlocal bias parameter, b_3nl, and the constant stochasticity term, N. The latter parameter affects only the smallest comoving separations, which are not considered in our analysis. As we verified, its impact on the cosmological outcomes of this work is in fact negligible. In the local Lagrangian framework the nonlocal bias terms can be written as a function of b₁ as follows (Chan et al. 2012; Saito et al. 2014):

$$\begin{eqnarray}&&{b}_{{s}^{2}}=-\displaystyle \frac{4}{7}({b}_{1}-1),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{b}_{3\mathrm{nl}}=\displaystyle \frac{32}{315}({b}_{1}-1).\end{eqnarray} \tag{ 12 }$$

The P_{δ δ} (k), P_{δ θ} (k), and P_{θ θ} (k) terms are estimated in the standard perturbation theory (SPT), which consists of expanding the statistics as a sum of infinite terms, corresponding to the n-loop corrections (see, e.g., Gil-Marín et al. 2012). Considering corrections up to the first loop order, the matter power spectrum can be modeled as follows:

$$\begin{eqnarray}&&{P}^{\mathrm{SPT}}(k)={P}_{{\rm{m}}}^{\mathrm{lin}}(k)+2{P}_{13}(k)+{P}_{22}(k),\end{eqnarray} \tag{ 13 }$$

where the one-loop correction terms are computed with the CPT Library ¹⁴ (Taruya & Hiramatsu 2008; Zhao et al. 2021). To test the model accuracy, we compared the outcomes of our reference analysis with the ones obtained by adopting the Bel et al. (2019) universal fitting functions for P_{δ θ} (k) and P_{θ θ} (k), finding negligible differences (see also Pezzotta et al. 2017; de la Torre et al. 2017). The other power spectrum terms in Equations (9) and (10), i.e., P_b2,δ (k), P_bs2,δ (k), P_b22(k), P_b2s2(k), P_bs22(k), P_b2,θ (k), and P_bs2,θ (k), are computed as a function of ${P}_{{\rm{m}}}^{\mathrm{lin}}(k)$, as prescribed in, e.g., Beutler et al. (2014) and Gil-Marín et al. (2014).

The damping term is assumed to be Lorentzian in Fourier space (see, e.g., de la Torre et al. 2017):

$$\begin{eqnarray}&&D(k,f,\mu ,{\sigma }_{v})=\displaystyle \frac{1}{1+{k}^{2}{f}^{2}{\mu }^{2}{\sigma }_{v}^{2}},\end{eqnarray} \tag{ 14 }$$

where σ_v is a nuisance parameter to marginalize over (Davis & Peebles 1983; Fisher et al. 1994; Zurek et al. 1994).

Finally, we compute the correction terms C_A and C_B in Equation (8) in SPT as follows (Taruya et al. 2010; de la Torre & Guzzo 2012):

$$\begin{eqnarray}\begin{array}{rcl}{C}_{A}(k,\mu ,f,{b}_{1}) & = & (k\mu f)\displaystyle \int \displaystyle \frac{{d}^{3}{\boldsymbol{p}}}{{\left(2\pi \right)}^{3}}\displaystyle \frac{{p}_{z}}{{p}^{2}}\\ & & \times \,\left[{B}_{\sigma }({\boldsymbol{p}},{\boldsymbol{k}}-{\boldsymbol{p}},-{\boldsymbol{k}})-{B}_{\sigma }({\boldsymbol{p}},{\boldsymbol{k}},-{\boldsymbol{k}}-{\boldsymbol{p}})\right],\end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{C}_{B}(k,\mu ,f,{b}_{1})={\left(k\mu f\right)}^{2}\int \displaystyle \frac{{d}^{3}{\boldsymbol{p}}}{{\left(2\pi \right)}^{3}}F({\boldsymbol{p}})F({\boldsymbol{k}}-{\boldsymbol{p}}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&F({\boldsymbol{p}})=\displaystyle \frac{{p}_{z}}{{p}^{2}}\left[{P}_{\delta \theta }(p)+f\displaystyle \frac{{p}_{z}^{2}}{{p}^{2}}{P}_{\theta \theta }(p)\right],\end{eqnarray} \tag{ 17 }$$

where B_σ is the cross bispectrum. The C_A and C_B terms are proportional to ${b}_{1}^{3}$ and ${b}_{1}^{4}$, respectively, and can be expressed as a power series expansion of b₁, f, and μ (see, e.g., García-Farieta et al. 2019, 2020, for more details).

We note that the eTNS model given by Equation (8) reduces to the Taruya et al. (2010) model if all the nonlinear bias terms are neglected, to the Scoccimarro (2004) model if also the C_A and C_B terms are neglected, and to the so-called dispersion model if both P_{δ θ} (k) and P_{θ θ} (k) are approximated as P_{δ δ} (k), which is valid in the linear regime (Kaiser 1987; Peacock & Dodds 1996).

The power spectrum multipoles can be estimated from Equation (8), as follows:

$$\begin{eqnarray}&&{P}_{l}(k)=\displaystyle \frac{2l+1}{2{\alpha }_{\perp }^{2}{\alpha }_{\parallel }}{\int }_{-1}^{1}{\rm{d}}\mu \,{P}^{s}({k}^{{\prime} },{\mu }^{{\prime} }){L}_{l}(\mu ),\end{eqnarray} \tag{ 18 }$$

where the Alcock & Paczynski (AP) (1979) geometric distortions, caused by a possibly incorrect assumption of the background cosmology used to convert cluster redshifts into comoving distances in Equation (3), are modeled by rescaling the wave numbers as follows (Beutler et al. 2014):

$$\begin{eqnarray}&&{k}^{{\prime} }=\displaystyle \frac{k}{{\alpha }_{\perp }}{\left[1+{\mu }^{2}\left(\displaystyle \frac{{\alpha }_{\perp }^{2}}{{\alpha }_{\parallel }^{2}}-1\right)\right]}^{1/2},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\mu }^{{\prime} }=\mu \displaystyle \frac{{\alpha }_{\perp }}{{\alpha }_{\parallel }}{\left[1+{\mu }^{2}\left(\displaystyle \frac{{\alpha }_{\perp }^{2}}{{\alpha }_{\parallel }^{2}}-1\right)\right]}^{-1/2},\end{eqnarray} \tag{ 20 }$$

with

$$\begin{eqnarray}&&{\alpha }_{\parallel }=\frac{{{\rm{H}}}^{\mathrm{fid}}(z){r}_{s}^{\mathrm{fid}}({z}_{d})}{{\rm{H}}(z){r}_{s}({z}_{d})},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\alpha }_{\perp }=\frac{{D}_{A}(z){r}_{s}^{\mathrm{fid}}({z}_{d})}{{D}_{A}^{\mathrm{fid}}(z){r}_{s}({z}_{d})},\end{eqnarray} \tag{ 22 }$$

where H^fid(z) and ${D}_{A}^{\mathrm{fid}}(z)$ are the fiducial values for the Hubble constant and angular diameter distance, respectively, and ${r}_{s}^{\mathrm{fid}}({z}_{d})$ is the fiducial sound horizon at the drag redshift assumed in the power spectrum template.

The corresponding 2PCF multipoles in configuration space read as

$$\begin{eqnarray}&&{\xi }_{l}(s)={i}^{l}{\int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}k}{2{\pi }^{2}}{k}^{2}{P}_{l}(k){j}_{l}({ks}),\end{eqnarray} \tag{ 23 }$$

where L_l (μ) are the Legendre polynomials and j_l are the spherical Bessel functions of order l.

Lastly, we assess the redshift-space wedges from the multipole moments through the following relation:

$$\begin{eqnarray}&&{\xi }_{w}(r)=\sum _{l}{\xi }_{l}(s){\bar{L}}_{l},\end{eqnarray} \tag{ 24 }$$

where ${\bar{L}}_{l}$ is the average value of the Legendre polynomials over the interval [μ₁, μ₂]. In particular, neglecting minor contributions from multipoles with l > 2 and considering the wedge width Δμ = 0.5, Equation (24) can be written as follows (Kazin et al. 2012):

$$\begin{eqnarray}&&\left(\begin{array}{l}{\xi }_{\parallel }\\ {\xi }_{\perp }\end{array}\right)=\left(\begin{array}{cc}1 & \tfrac{3}{8}\\ 1 & -\tfrac{3}{8}\end{array}\right)\left(\begin{array}{c}{\xi }_{0}\\ {\xi }_{2}\end{array}\right).\end{eqnarray} \tag{ 25 }$$

A full validation of the implemented likelihood algorithms on simulated galaxy and cluster catalogs will be presented in a forthcoming paper (García-Farieta et al. 2021, in preparation).

As a common practice, and to directly compare to previous similar analyses performed on galaxy samples, we parameterize the model as a function of σ_v and the three parameter products $\left[f{\sigma }_{8},{b}_{1}{\sigma }_{8},{b}_{2}{\sigma }_{8}\right]$ (but see Sanchez 2020), fixing the other parameters to Planck 18 cosmology.

In this work we focus on scales smaller than the BAO ones, which does not allow us to put strong enough constraints on the geometric distortions that are degenerate with RSD (e.g., Taruya et al. 2011). Nevertheless, to marginalize over the AP distortion parameters, we allow them to vary, considering Gaussian priors with a standard deviation of 0.01 (see, e.g., de la Torre et al. 2017, for a similar approach).

The posterior distribution constraints on these parameters are assessed through a Markov Chain Monte Carlo (MCMC) statistical analysis, assuming a standard Gaussian likelihood:

$$\begin{eqnarray}&&-2\mathrm{ln}{ \mathcal L }=\sum _{i,j=1}^{N}\left[{\xi }_{k}^{D}({s}_{i})-{\xi }_{k}^{M}({s}_{i})\right]{C}_{k}^{-1}({s}_{i},{s}_{j})\left[{\xi }_{k}^{D}({s}_{j})-{\xi }_{k}^{M}({s}_{j})\right],\end{eqnarray} \tag{ 26 }$$

where N is the number of bins at which the wedges are computed, and the superscripts D and M refer to data and model, respectively.

The covariance matrix, C_k , which measures the variance and correlation between 2PCF wedge bins, is defined as follows:

$$\begin{eqnarray}&&{C}_{w}\left({s}_{i},{s}_{j}\right)=\displaystyle \frac{1}{{N}_{R}-1}\sum _{n=1}^{{N}_{R}}\left[{\xi }_{w}^{n}\left({s}_{i}\right)-{\overline{\xi }}_{w}\left({s}_{i}\right)\right]\left[{\xi }_{w}^{n}\left({s}_{j}\right)-{\overline{\xi }}_{w}\left({s}_{j}\right)\right].\end{eqnarray} \tag{ 27 }$$

The indices i and j run over the 2PCF wedge bins, while w = 0, 0.5 refers to each clustering wedge. In both cases, ${\bar{\xi }}_{w}\,=1/{N}_{R}{\sum }_{n=1}^{{N}_{R}}{\xi }_{w}^{n}$ is the average wedge of the 2PCF, and N_R = 100 is the number of realizations obtained by resampling the catalogs with the bootstrap method. We correct the inverse covariance matrix estimator to account for the finite number of realizations as in Hartlap et al. (2007), while the parameter uncertainties are corrected to take into account the uncertainties in the covariance estimate as in Percival et al. (2014).

Similarly to RSD analyses of galaxy clustering, we adopt large flat priors on fσ₈, b₂ σ₈, and σ_v , specifically $f{\sigma }_{8}=\left[0,1\right]$, ${b}_{2}{\sigma }_{8}=\left[-10,10\right]$, and ${\sigma }_{v}=\left[0,100\right]{\rm{Mpc}}\,{h}^{-1}$, respectively. While σ₈ could be constrained directly from the cluster mass function, this would have required an accurate knowledge of the cluster selection function to avoid systematic uncertainties. To provide conservative linear growth constraints, we prefer to focus the current analysis on cluster clustering, setting all the cosmological parameters, including σ₈, to Planck 18 values. The constraint on the linear growth rate we will derive in this paper has to be considered in this respect, though to compare to previous analyses we will express our results in terms of fσ₈.

Differently from galaxy clustering analyses, we can set a strong prior on b₁ σ₈ thanks to the knowledge of galaxy cluster masses inferred from weak-lensing scaling relations. In particular, the b₁ σ₈ prior is centered on the effective linear bias of the cluster sample, which is estimated as in Marulli et al. (2018):

$$\begin{eqnarray}&&{b}_{\mathrm{eff}}^{2}=\langle b({\tilde{M}}_{i},{z}_{i})b({\tilde{M}}_{j},{z}_{j})\rangle ,\end{eqnarray} \tag{ 28 }$$

where the linear bias of each cluster, b, is computed with the Tinker et al. (2010) model, while $\tilde{{M}_{i}}$ and $\tilde{{M}_{j}}$ are the masses of the two clusters of each pair, at redshifts z_i and z_j , respectively, estimated from the weak-lensing cluster mass scaling relation given by Equation (1). This represents the key difference with respect to analogous RSD analyses of galaxy samples, as in those cases the masses of the dark matter haloes hosting the galaxies are unknown and no priors can be reliably assumed on the bias of the catalog. We will discuss the impact of this assumption in Section 5.

Drawing a set of mass samples from the scaling relation, we computed the average bias and variance of our cluster catalog, which correspond to b₁ σ₈ = 1.722 ± 0.002. We consider a Gaussian prior on b₁ σ₈ centered on the latter value. To provide conservative constraints, we adopt a prior width 5 times larger than the estimated standard deviation, i.e., 0.01, to include possible systematic uncertainties in the adopted bias model and scaling relation.

As described in Section 3, the clustering wedges are computed by integrating the redshift-space 2PCF, ξ(s, μ), over two bins of μ. In Figure 3 we present the redshift-space 2PCF of the selected cluster sample in two coordinate systems, i.e., as a function of perpendicular (s_⊥) and parallel (s_∥) separations to the line of sight (Cartesian coordinates), and as a function of distance modulus (s) and cosine of the angle (μ) between the line of sight (polar coordinates). In real space the contour lines of the former statistics would be circular, while the ones of the latter statistics would be straight. RSD introduces anisotropies in the derived map that warp these contour lines, an effect that depends directly on the value of the linear growth rate of cosmic structures.

Zoom In Zoom Out Reset image size

The shape of the Cartesian 2PCF of the selected clusters shown in the left panel of Figure 3 appears similar to the one of galaxies, as expected (e.g., Alam et al. 2017; Marulli et al. 2017). In fact, as described in Section 2.2, the 2PCF of galaxy clusters we measure in this work coincides with the 2PCF of BCGs, by construction. Nevertheless, the cosmological analysis of this data set provides a clear advantage as, differently from galaxy clustering studies, we can infer in this case the linear bias of the sample from the richness-mass scaling relation of the galaxy clusters hosting the selected BCGs, as explained in Section 4.2.

The Fingers-of-God distortions at small scales due to incoherent peculiar motions are not completely negligible, though much less strong relative to the case of the 2PCF of lower biased tracers (see the discussion in Marulli et al. 2017). The polar 2PCF shown in the right panel of Figure 3 is the statistics we integrate along the μ direction to compute the wedges.

Figure 4 shows the redshift-space transverse and radial wedges of the cluster 2PCF, defined by Equation (7). The horizontal and vertical error bars are the standard deviation around the mean pair separation in each bin, and the diagonal values of the bootstrap covariance matrix, respectively. In real space, the radial and transverse wedges would be statistically equal, for isotropy. On the other hand, redshift-space anisotropies make these two statistics significantly different, as shown in Figure 4.

Zoom In Zoom Out Reset image size

The correlation matrix, which is ${C}_{i,j}/\sqrt{{C}_{i,j}{C}_{j,i}}$, of the redshift-space transverse and radial wedges is shown in Figure 5. The algorithms to estimate the covariance matrix have been highly validated in previous works on both simulations and real cluster catalogs (e.g., Veropalumbo et al. 2014, 2016; Marulli et al. 2017; García-Farieta et al. 2020). To further check the results of the current analysis, we compare the bootstrap error estimates with the ones obtained with either the jackknife method or the analytic Gaussian model provided by Grieb et al. (2016) (for the theoretical modeling of non-Gaussian contributions to the covariance, which are neglected in the current analysis, see Sugiyama et al. 2020). The diagonal values of the bootstrap, jackknife, and analytic matrices are compared in Figure 6. The estimated bootstrap clustering uncertainties are statistically consistent with the analytic ones at scales larger than about 20 h⁻¹ Mpc , while the jackknife uncertainties appear slightly larger.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bootstrap method allows us to draw a greater number of realizations relative to the jackknife method, providing a smoother covariance matrix, whose inverse is less affected by numerical noise. Moreover, it does not depend on free parameters, differently from the analytic covariance matrix, which depends on the sample bias and on the effective area of the survey, whose values are inferred within uncertainties. For the above reasons, in this work we rely on the bootstrap covariance uncertainties.

Following the method described in Section 4, we perform a joint statistical analysis of the redshift-space radial and transverse wedges of the selected spectroscopic cluster catalog, in the standard ΛCDM framework. The best-fit eTNS model obtained from the median of the MCMC posterior distribution is reported in Figure 4, together with its 68% uncertainty region. The fit is performed in the comoving scale range $10\lt s\left[{\rm{Mpc}}\,{h}^{-1}\right]\lt 80$. The model appears statistically consistent with the measurements in the scale range considered. We note only a minor mismatch in the radial wedge at $40\lt s\left[{\rm{Mpc}}\,{h}^{-1}\right]\lt 60$, though it is not statistically significant considering the covariance in the measurements. The reduced χ² at the best-fit eTNS model is ${\tilde{\chi }}^{2}=-2\mathrm{ln}{ \mathcal L }/{\rm{d.o.f.}}=0.71$ (where d.o.f. are the degrees of freedom of the data sample).

The marginalized posterior distributions on fσ₈, b₁ σ₈, b₂ σ₈, and σ_v are reported in Figure 7, together with the 68% and 95% posterior confidence regions. We assess the best-fit values and marginalized constraints from the median and percentile values of the posterior distribution. We note in particular that the fσ₈ posterior is consistent with a Gaussian distribution. We get fσ₈ = 0.44 ± 0.05, at the mean pair redshift z = 0.275. The relative statistical uncertainty is about 10%, which is a remarkable result considering the sparsity of the spectroscopic cluster sample considered. This is caused by the narrow Gaussian prior on the effective bias of the sample, assessed through the cluster mass-richness scaling relation, as described in Section 4. In fact, the posterior distribution we get on b₁ σ₈ is statistically indistinguishable from the assumed Gaussian prior distribution.

Zoom In Zoom Out Reset image size

We tested the impact of this assumption by running the statistical analysis assuming prior standard deviations of 0.002, 0.05, and 0.1. In particular, the latter prior width accounts for possible bias model uncertainties of about 5%. We obtained fσ₈ = 0.44 ± 0.05, f σ₈ = 0.42 ± 0.06, and f σ₈ = 0.41 ± 0.06, respectively, which are all statistically consistent, considering current measurement uncertainties. Assuming instead a flat prior distribution on b₁ σ₈ in $\left[0.1,5\right]$, we get f σ₈ = 0.40 ± 0.07, which is still statistically consistent with our fiducial result, though with a larger relative error of about 18%.

All the results obtained for the different prior assumptions considered are statistically consistent, showing that the accuracy in the mass estimates and in the linear bias model is high enough for the current analysis. Improving the mass and bias modeling will become crucial instead for cluster clustering analyses of next-generation surveys.

The reference value of the quadratic bias factor, b₂ σ₈, reported in Figure 7 is computed with the polynomial relation provided by Lazeyras et al. (2016):

$$\begin{eqnarray}&&{b}_{2}({b}_{1})=0.412-2.143\,{b}_{1}+0.929\,{b}_{1}^{2}+0.008\,{b}_{1}^{3},\end{eqnarray} \tag{ 29 }$$

while the reference value of σ_v is estimated in linear theory as follows (Taruya et al. 2010):

$$\begin{eqnarray}&&{\sigma }_{{\rm{v}},\mathrm{lin}}^{2}=\displaystyle \frac{1}{3}\int \displaystyle \frac{{d}^{3}{\boldsymbol{q}}}{{\left(2\pi \right)}^{3}}\displaystyle \frac{{P}_{{\rm{m}}}^{\mathrm{lin}}(q,z)}{{q}^{2}}.\end{eqnarray} \tag{ 30 }$$

To test the robustness of our results, we repeated the analysis fitting the wedges in narrower scale ranges, considering either jackknife or analytic clustering uncertainties instead of boostrap, and modeling small-scale random motions with a Gaussian damping term instead of a Lorentzian one (Marulli et al. 2012; Sridhar et al. 2017; García-Farieta et al. 2020). Overall we found consistent results, within the 68% confidence region.

As discussed in Section 4.1, we considered tight Gaussian priors on the AP distortion parameters with standard deviation of 0.01. To investigate the impact of this assumption, we run our analysis also for standard deviation values of 0, 0.05, and 0.1, obtaining fσ₈ = 0.44 ± 0.05, f σ₈ = 0.4 ± 0.1, and f σ₈ = 0.3 ± 0.2, respectively. As expected, the impact of this prior is significant. Since the current analysis focuses on scales below the BAO peak, we cannot break the degeneracy between RSD and geometric distortions. Constraints from a joint RSD+BAO analysis of both the two-point and three-point correlation functions of this cluster sample are presented in A. Veropalumbo et al. (2021, in preparation).

In Figure 8 we compare the fσ₈ constraint obtained in this work with a large collection of measurements at different redshifts from galaxy, quasar, and cosmic void samples and other tracers. The data shown provide the key observables to test the gravity theory on the largest cosmological scales. ¹⁵ Figure 8 summarizes our current understanding of the cosmological evolution of the linear growth rate of cosmic structures. Specifically, the data reported are from 6dF Galaxy Survey (6dFGS) (Beutler et al. 2012; Huterer et al. 2017; Adams & Blake 2017); First Amendment set of supernovae (SNe) peculiar velocities (Turnbull et al. 2012); Two Micron All Sky Survey (2MASS; Davis et al. 2011); Galaxy And Mass Assembly (GAMA) (Blake et al. 2013); WiggleZ (Blake et al. 2012); VIMOS-VLT Deep Survey (VVDS) (Guzzo et al. 2008); VIMOS Public Extragalactic Redshift Survey (VIPERS) (de la Torre et al. 2013, 2017; Pezzotta et al. 2017; Hawken et al. 2017; Mohammad et al. 2018); FastSound (Okumura et al. 2016); SDSS+BOSS+eBOSS (Howlett et al. 2015; Alam et al. 2017; Bautista et al. 2021; Gil-Marín et al. 2020; Tamone et al. 2020; de Mattia et al. 2021; Hou et al. 2021; Neveux et al. 2020; Nadathur et al. 2019, 2020; Alam et al. 2021). The latest, most constraining measurements from SDSS+BOSS+eBOSS are highlighted in green.

Zoom In Zoom Out Reset image size

The data are compared to the standard ΛCDM + general relativity (GR) predictions, which is f σ₈ = Ω_M (z)^0.545 σ₈(z), where σ₈ and Ω_M are computed assuming the Planck 18 cosmological parameters. The fσ₈ constraint obtained in this work appears fully consistent with the other data, with a competitive statistical uncertainty. By comparison, we also plot the model with $\gamma ={0.65}_{-0.04}^{+0.05}$, which provides a better fit to the H(z)/H₀ − f(z)σ₈(z) diagram, as found by Moresco & Marulli (2017).

The exploitation of growth rate measurements at different redshifts to discriminate among alternative dark energy and modified gravity models would require a detailed study which is outside the scope of this work. Nevertheless, to highlight the constraining power of current cluster clustering measurements, we compare in Figure 9 the fσ₈ constraint obtained in this work with the predictions of three popular alternative models, which is the f(R) model (e.g., De Felice & Tsujikawa 2010), the coupled dark energy (cDE) model (Wetterich 1995; Amendola 2000) and the Dvali–Gabadaze−Porrati (DGP) model (Dvali et al. 2000).

Zoom In Zoom Out Reset image size

The linear growth rate in f(R) and DGP models can be expressed with the so-called γ parameterization:

$$\begin{eqnarray}&&f={{\rm{\Omega }}}_{M}{\left(z\right)}^{\gamma \left(z\right)},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&\gamma (z)={\gamma }_{0}+{\gamma }_{1}\displaystyle \frac{z}{1+z}.\end{eqnarray} \tag{ 32 }$$

In particular, we consider the f(R) model by Hu & Sawicki (2007), setting the two free model parameters to n = 2 and λ = 3 following Di Porto et al. (2012), which corresponds to γ₁ = 0.43 and γ₂ = −0.2 in the limit of small, still linear, scales (Gannouji et al. 2009). For the DGP model we consider the flat-space case in which γ₁ = 0.633 and γ₂ = 0.041 (Maartens & Majerotto 2006; Fu et al. 2009). Finally, we model the linear growth rate in the cDE model with the so-called η-parameterization (di Porto & Amendola 2008):

$$\begin{eqnarray}&&f={{\rm{\Omega }}}_{M}{\left(z\right)}^{\gamma }\left(1+\eta \right),\end{eqnarray} \tag{ 33 }$$

where $\eta =c{\beta }_{c}^{2}$ quantifies the coupling strength. In particular, we consider the case with γ = 0.56, c = 2.1, and β_c = 0.16, which implies η = 0.056 (di Porto & Amendola 2008; Di Porto et al. 2012).

As shown in Figure 9, the current measurement uncertainties are not low enough to discriminate between these alternative cosmological scenarios at a sufficient statistical level. Next-generation experiments will have instead the required accuracy to achieve this key scientific task (e.g., Amendola et al. 2018).