Cold dark energy constraints from the abundance of galaxy clusters

TL;DR

This paper investigates the constraints on cold dark energy with negligible sound speed using galaxy cluster abundance data, highlighting its clustering properties and potential to improve cosmological parameter estimates.

Contribution

It introduces a new mass function calibration for cold dark energy models and demonstrates their distinguishability from standard models using cluster data.

Findings

Cold dark energy models show up to 1σ shift from quasi-homogeneous models.

Cluster data combined with other probes constrains dark energy parameters effectively.

Forecasts indicate 50% tighter constraints on key parameters with DES data.

Abstract

We constrain cold dark energy of negligible sound speed using galaxy cluster abundance observations. In contrast to standard quasi-homogeneous dark energy, negligible sound speed implies clustering of the dark energy fluid at all scales, allowing us to measure the effects of dark energy perturbations at cluster scales. We compare those models and set the stage for using non-linear information from semi-analytical modelling in cluster growth data analyses. For this, we recalibrate the halo mass function with non-linear characteristic quantities, the spherical collapse threshold and virial overdensity, that account for model and redshift dependent behaviours, as well as an additional mass contribution for cold dark energy. We present the first constraints from this cold dark matter plus cold dark energy mass function using our cluster abundance likelihood, which self-consistently accounts for selection effects, covariances and systematic uncertainties. We combine cluster growth data with CMB, SNe Ia and BAO data, and find a shift between cold versus quasi-homogeneous dark energy of up to $1\sigma$. We make a Fisher matrix forecast of constraints attainable with cluster growth data from the on-going Dark Energy Survey (DES). For DES, we predict $\sim$50$\%$ tighter constraints on $\left(\Omega_\mathrm{m},w \right)$ for cold dark energy versus $w$CDM models, with the same free parameters. Overall, we show that cluster abundance analyses are sensitive to cold dark energy, an alternative, viable model that should be routinely investigated alongside the standard dark energy scenario.

Figures

Tables

Table 1. Table 1: Marginalized best-fitting values and 68.3 per cent confidence intervals for σ 8 subscript 𝜎 8 \sigma_{8} , Ω m subscript Ω m \Omega_{\mathrm{m}} and w 𝑤 w , for both cold dark energy with sound speed c s = 0 subscript 𝑐 s 0 c_{\mathrm{s}}=0 and quasi-homogeneous dark energy with c s = 1 subscript 𝑐 s 1 c_{\mathrm{s}}=1 . Results are shown for cluster+CMB+BAO+SNIa data as described in Section  4 . Note for example the shift to lower σ 8 subscript 𝜎 8 \sigma_{8} , higher Ω m subscript Ω m \Omega_{\mathrm{m}} and higher w 𝑤 w for c s = 0 subscript 𝑐 s 0 c_{\mathrm{s}}=0 as compared to c s = 1 subscript 𝑐 s 1 c_{\mathrm{s}}=1 .

Analysis	${\sigma }_8$	${\Omega }_m$	$w$
Clusters + CMB + BAO + SNIa, $c_s=0$	0.806 $\pm$ 0.014	0.302 $\pm$ 0.013	–1.14 $\pm$ 0.05
Clusters + CMB + BAO + SNIa, $c_s=1$	0.823 $\pm$ 0.017	0.296 $\pm$ 0.013	–1.19 $\pm$ 0.07

Table 2. Table 2: Survey specifications used for the Fisher matrix forecast of cosmological parameters with cluster number counts (see Section  6 ). Shown are the sky fraction f sky subscript 𝑓 sky f_{\mathrm{sky}} , the threshold mass M thr subscript 𝑀 thr M_{\mathrm{thr}} , the logarithmic mass binning Δ ​ log ⁡ ( M / M ⊙ ) Δ 𝑀 subscript 𝑀 direct-product \Delta\log\left(M/M_{\odot}\right) , the maximum survey redshift z max subscript 𝑧 max z_{\mathrm{max}} and the redshift binning Δ ​ z Δ 𝑧 \Delta z .

Survey	$f_{\mathrm{sky}}$ (per cent)	$M_{\mathrm{thr}}$ ($M_{\odot }$)	$\Delta \log \left(M/M_{\odot }\right)$	$z_{\max }$	$\Delta z$
DES	12	$1\times {10}^{14}$	0.2	1.0	0.05

Table 3. Table 3: Marginalized 68.3 per cent confidence intervals for DES, for the Fisher matrix forecasts of Section  6 for both c s = 0 subscript 𝑐 s 0 c_{\mathrm{s}}=0 and c s = 1 subscript 𝑐 s 1 c_{\mathrm{s}}=1 . We also show the FoM in the ( Ω m , w ) subscript Ω m 𝑤 \left(\Omega_{\rm m},w\right) -plane as defined in Section  6 . The fiducial parameter values are { Ω m , σ 8 , w } = { 0.287 , 0.87 , − 1.124 } subscript Ω m subscript 𝜎 8 𝑤 0.287 0.87 1.124 \left\{\Omega_{\mathrm{m}},\sigma_{8},w\right\}=\left\{0.287,0.87,-1.124\right\} . The marginalized confidence intervals for the offset in the mass calibration B M , 0 subscript 𝐵 𝑀 0 B_{M,0} , which was left to vary at the 3 % percent 3 3\% level, are Δ ​ B M , 0 = 1.6 % Δ subscript 𝐵 𝑀 0 percent 1.6 \Delta B_{M,0}=1.6\% for c s = 0 subscript 𝑐 s 0 c_{\mathrm{s}}=0 and Δ ​ B M , 0 = 0.2 % Δ subscript 𝐵 𝑀 0 percent 0.2 \Delta B_{M,0}=0.2\% for c s = 1 subscript 𝑐 s 1 c_{\mathrm{s}}=1 . The tighter constraints on this parameter for c s = 1 subscript 𝑐 s 1 c_{\mathrm{s}}=1 are due to the larger number of clusters detected in this case.

	$\Delta {\sigma }_8$ (10-3)	$\Delta {\Omega }_m$ (10-3)	$\Delta w$ (10-3)	FoM (103)
$c_s=0$	5.2	3.1	31.4	11.4
$c_s=1$	2.9	6.1	43.3	10.5

Equations

$$\frac{\partial {\rho }_i}{\partial t}+{\nabla }_{\mathbf{r}}\cdot \left({\rho }_i+p_i\right){\mathbf{v}}_i=0,$$

$$\frac{\partial {\mathbf{v}}_i}{\partial t}+\left({\mathbf{v}}_i\cdot {\nabla }_{\mathbf{r}}\right){\mathbf{v}}_i+\frac{{\nabla }_{\mathbf{r}}{p}_i+{\mathbf{v}}_i\dot{p_i}}{{\rho }_i+p_i}+{\nabla }_{\mathbf{r}}\Phi =0,$$

$${\dot{\delta }}_i$$

$${\dot{\theta }}_i$$

$${\nabla }^2\Phi =-4\pi G\sum _I\left(1+3c_{\mathrm{s},i}^2\right)a^2{\bar{\rho }}_i{\delta }_i,$$

$${\ddot{\delta }}_{\mathrm{m}}$$

$${\dot{\delta }}_{\mathrm{de}}$$

$${\delta }_{\mathrm{de},i}=\frac{1+w}{1-3w}{\delta }_{\mathrm{m},i},$$

$$\frac{\ddot{R}}{R}=-\frac{4\pi G}{3}\left({\rho }_{\mathrm{m}}+{\rho }_{\mathrm{de}}+3{\bar{p}}_{\mathrm{de}}\right),$$

$${\dot{\rho }}_i+3\frac{\dot{R}}{R}\left({\rho }_i+p_i\right)=0.$$

$${\Delta }_{\mathrm{vir}}=\left({\delta }_{\mathrm{NL},\mathrm{vir}}+1\right)=\left({\delta }_{\mathrm{ta}}+1\right){\left(\frac{a_{\mathrm{vir}}}{a_{\mathrm{ta}}}\right)}^3{\left(\frac{R_{\mathrm{ta}}}{R_{\mathrm{vir}}}\right)}^3,$$

$${\Delta }_{\mathrm{vir}}\left(z\right)=\frac{\left[18{\pi }^2+a\left(1-{\Omega }_{\mathrm{m}}\left(z\right)\right)+b{\left(1-{\Omega }_{\mathrm{m}}\left(z\right)\right)}^2\right]}{{\Omega }_{\mathrm{m}}\left(z\right)},$$

$$M_{\mathrm{vir},\mathrm{de}}=4\pi {\int }_0^{R_{\mathrm{vir}}}\mathrm{d}R{R}^2\delta {\rho }_{\mathrm{de}},$$

$$\frac{\mathrm{d}{n}_{\mathrm{cal}}}{\mathrm{d}M}\left(M,z\right)=\frac{\mathrm{d}{n}_{\mathrm{ST}}/\mathrm{d}M\left(M,z;c_{\mathrm{s}}=0\right)}{\mathrm{d}{n}_{\mathrm{ST}}/\mathrm{d}M\left(M,z;c_{\mathrm{s}}=1\right)}\times \frac{\mathrm{d}{n}_{\mathrm{T}}}{\mathrm{d}M}\left(M,z\right),$$

$$\frac{\mathrm{d}{n}_T}{\mathrm{d}M}\left(M,z\right)=f\left(\sigma \right)\frac{{\bar{\rho }}_{\mathrm{m}}}{M}\frac{\mathrm{d}\log {\sigma }^{-1}}{\mathrm{d}M}.$$

$$f\left(\sigma \right)=A\left[{\left(\frac{\sigma }{b}\right)}^{-a}+1\right]e^{-c/{\sigma }^2},$$

$${\sigma }^2\left(R\right)=\frac{1}{2{\pi }^2}\int P\left(k\right){\left|W\left(kR\right)\right|}^2{k}^2\mathrm{d}k.$$

$$\frac{\mathrm{d}{n}_{\mathrm{ST}}}{\mathrm{d}M}\left(M,z\right)=\nu f\left(\nu \right)\frac{{\bar{\rho }}_{\mathrm{m}}}{M^2}\frac{\mathrm{d}\log \nu }{\mathrm{d}\log M}.$$

$$\nu f\left(\nu \right)=A\sqrt{\frac{2a{\nu }^2}{\pi }}\left[1+{\left(a{\nu }^2\right)}^{-p}\right]\exp \left[-a{\nu }^2\right],$$

$$F_{ij}=\sum _{l,m}=\frac{\partial N_{l,m}}{\partial p_i}{C}_{l,m}^{-1}\frac{\partial N_{l,m}}{\partial p_j},$$

$$N_{l,m}=$$

$${\int }_{M_{l,m}^{\mathrm{ob}}}^{M_{l,m+1}^{\mathrm{ob}}}\frac{\mathrm{d}{M}^{\mathrm{ob}}}{M^{\mathrm{ob}}}{\int }_0^{+\infty }\mathrm{d}Mn\left(M,z\right)p\left(M^{\mathrm{ob}}|M;z\right),$$

$$\frac{\mathrm{d}V}{\mathrm{d}z}=\frac{4\pi c{\left(1+z\right)}^2}{H\left(z\right)}{d}_{\mathrm{A}}^2\left(z\right),$$

$$p\left(M^{\mathrm{ob}}|M;z\right)=\frac{\exp \left[-x^2\left(M^{\mathrm{ob}},M,z\right)\right]}{\sqrt{2{\sigma }_{\ln M}^2(z)}},$$

$$x\left(M^{\mathrm{ob}},M,z\right)=\frac{\ln M^{\mathrm{ob}}-\ln M_{\mathrm{bias}}(z)-\ln M}{\sqrt{2{\sigma }_{\ln M}^2(z)}}.$$

$$N_{l,m}=$$

$${\int }_0^{+\infty }\mathrm{d}Mn\left(M,z\right)\left[\mathrm{erfc}\left(x_m\right)-\mathrm{erfc}\left(x_{m+1}\right)\right],$$

$$\ln M_{\mathrm{bias}}\left(z\right)$$

$${\sigma }_{\ln M}^2\left(z\right)$$

$${\mathbf{p}}_{\mathbf{N}}=\left\{B_{M,0}=0,\alpha =0,{\sigma }_{\ln M,0}=0.4,\beta =0.0\right\}.$$

Full text

††pubyear: 2017††pagerange: Cold dark energy constraints from the abundance of galaxy clusters–Cold dark energy constraints from the abundance of galaxy clusters

1 Introduction

Cosmology has entered a phase where ever more precise measurements enable us to put increasingly tight constraints on model parameters. Since the first evidence was found for late-time accelerated expansion (Riess et al., 1998; Perlmutter et al., 1999), the possibilities of either a cosmological constant, dynamical dark energy or modifications of gravity has been a question central to cosmology (for reviews see e.g. Copeland et al., 2006; De Felice & Tsujikawa, 2010; Clifton et al., 2012; Joyce et al., 2015). To guarantee the accuracy of such precise constraints, our modelling toolkit needs to be extended. For this, it is key to explore a wider range of model and parameter spaces, while correctly translating model characteristics into quantities testable against data. The objective is to maximize the information gain and not overlook distinctive signatures.

Both in the linear and non-linear regimes, galaxy cluster surveys are competitive probes of cosmological models and fundamental physics. This has been shown in results ranging from estimating standard cosmological parameters (Mantz et al., 2008; Vikhlinin et al., 2009; Mantz et al., 2010a; Rozo et al., 2010; Allen et al., 2011; Benson et al., 2013; Mantz et al., 2014, 2015; Planck Collaboration et al., 2016; Applegate et al., 2016; de Haan et al., 2016) to constraining non-Gaussianities of primordial density fluctuations (Sartoris et al., 2010a; Shandera et al., 2013a; Mana et al., 2013) and testing predictions of General Relativity and modified gravity scenarios (Schmidt et al., 2009a; Rapetti et al., 2010; Lombriser et al., 2012; Rapetti et al., 2013; Cataneo et al., 2015). The high-mass end of the halo mass function (HMF), which can be constrained by observations of galaxy clusters, is particularly sensitive to cosmological models through both the background evolution and the linear and non-linear growth of structure. For a vanilla model with a cosmological constant and cold dark matter ($\Lambda$CDM), the HMF has been carefully modelled and calibrated (Sheth & Tormen, 1999; Tinker et al., 2008; Jenkins et al., 2001; Maggiore & Riotto, 2010; Tinker et al., 2010; Corasaniti & Achitouv, 2011; Despali et al., 2016; Bocquet et al., 2016). Also, the calibration of mass function parameters with observations has been demonstrated (Castro et al., 2016; Mehrabi et al., 2017). Continued efforts have been directed to the modelling of the mass function for extended theories (Schmidt et al., 2009b; Bhattacharya et al., 2011; Cui et al., 2012; Barreira et al., 2013; Kopp et al., 2013; Cataneo et al., 2016; Nazari-Pooya et al., 2016).

For many models beyond $\Lambda$CDM, dark energy can be effectively parametrized as a fluid with an equation of state, $w$, and a sound speed $c_{\mathrm{s}}^{2}$. When the sound speed is non-adiabatic, i.e., when dealing with imperfect fluids, $w$ and $c_{\mathrm{s}}$ are not necessarily directly related to each other as they are in the adiabatic case. For dynamical dark energy models, a sound speed different from one (speed of light), even time- and scale-dependent, is quite natural. One of the simplest examples for dark energy models with a varying sound speed is Quintessence with non-canonical kinetic terms, known as K-essence (Armendariz-Picon et al., 1999, 2000). Beyond K-essence, in principle any model with an additional propagating scalar degree of freedom exhibits sound speeds different from one, for example the Horndeski class of scalar-tensor theories (Horndeski, 1974; Kobayashi et al., 2011; De Felice et al., 2011).

First attempts to constrain the sound speed of dark energy at the linear level came from Weller & Lewis (2003) using cosmic microwave background (CMB), large-scale structure (LSS) and supernova data, Bean & Doré (2004) using CMB and CMB LSS cross-correlation data, and Hannestad (2005) using CMB, galaxy clustering and weak-lensing data. Later, Abramo et al. (2009) forecasted sound speed constraints from galaxy cluster number counts, Appleby et al. (2013) considered constraints possible from number counts and the CMB, and  Hojjati & Linder (2016) studied the potential use of CMB lensing data for this purpose.  Creminelli et al. (2010b) and Batista & Pace (2013) extended the analysis of dark energy models with negligible sound speed to the non-linear level of structure formation utilizing the spherical collapse model.  Basse et al. (2011) and Basse et al. (2014) followed a similar approach but for dark energy models with an arbitrary speed of sound, and also forecasted parameter constraints.

The goal of this paper is to capture the rich non-linear information of structure formation imprinted on cluster growth data. For this, we implement a semi-analytical framework that incorporates the dominant effects of cold dark energy into the HMF. Using this framework, we present the first observational constraints on cold (clustering) dark energy ($c_{\mathrm{s}}^{2}=0$). We also compare these results with those for quasi-homogeneous dark energy ($c_{\mathrm{s}}^{2}=1$) to investigate the impact of the common assumption of $c_{\mathrm{s}}^{2}=1$ when constraining $w$ and the other relevant standard parameters.

This paper is organized as follows. Section 2 describes the semi-analytical framework we use to calculate the non-linear model characteristics of interest for both cold and quasi-homogeneous dark energy, and illustrates their behaviour. In Section 3 we proceed to recalibrate the HMF with these non-linear quantities. Section 4 briefly describes the data we use. In Section 5 we present our results on standard cosmological parameters for both cold and quasi-homogeneous dark energy, and in Section 6 forecast such constraints for DES. We summarize our findings and discuss the broader implications in Section 7.

2 Non-linear characteristics

In this section we review the effects of assuming cold dark energy on non-linear quantities such as the cosmology-dependent linear threshold of collapse, the virial overdensity and the cold dark energy mass contribution to virialized objects. We use the spherical collapse model to calculate the perturbations stemming from dark energy’s being clustering instead of quasi-homogeneous. This enables us to recalibrate the cluster mass function for cold dark energy by incorporating into it these non-linear quantities.

Cold dark energy designates a dark energy fluid whose sound speed is extremely low, i.e. approaching the limit of zero sound speed. A varying, and possibly low, sound speed is natural in scalar-tensor theories, such as for the simple example K-essence, where it avoids pathologies like ghost degrees of freedom or instabilities during the so-called phantom crossing (a crossing of the $w=-1$ line during the evolution of the K-essence field to negative values; Creminelli et al., 2009). When analysing data, it might therefore be advisable to test models with sound speed zero.

In the following we will assume a dark energy fluid with an effective sound speed $c_{\mathrm{s}}^{2}$ and equation of state $w$, at the limits where $c_{\mathrm{s}}^{2}\rightarrow 0$ for cold dark energy and $c_{\mathrm{s}}^{2}=1$, i.e., the speed of light, for quasi-homogeneous dark energy. The effective sound speed is defined here as $c_{\mathrm{s}}^{2}=\delta p_{\mathrm{e}}/\delta\rho_{e}$, with $\delta p_{\mathrm{e}}$ and $\delta\rho_{e}$ being the pressure and density perturbations of the dark energy fluid, respectively. This relation is more general than that of the adiabatic sound speed $c_{\mathrm{a}}^{2}$, where $c_{\mathrm{s}}^{2}=c_{\mathrm{a}}^{2}$, with $c_{\mathrm{a}}^{2}=w$, only for perfect fluids, but not necessarily for more general theories of dark energy. We restrict our current analysis to an effective constant equation of state $w=p_{\mathrm{e}}/\rho_{e}$, where $p_{\mathrm{e}}$ is the pressure and $\rho_{e}$ the energy density of the fluid. Assuming $w$ constant certainly is an approximative choice for models of sound speed different from one; leaving $w$ to vary together with the sound speed will be interesting for future studies. Here we aim at a first analysis exploring the effects of using the two most extreme possibilities for the sound speed. Also, judging from current observational constraints, we expect only a moderate variation of $w$ away from –1. In addition, we assume negligible anisotropic stress and a flat Friedmann–Lemaître–Robertson–Walker (FLRW) background.

2.1 Fluid equations and spherical collapse

Here, we describe the set of equations we employ within the spherical collapse model to derive the non-linear quantities that describe the cosmology-dependent formation of bound structures, which are needed to recalibrate the HMF (Section 3). These are the density threshold of collapse, the overdensity at virialization and a cold dark energy mass contribution to the bound objects. Their behaviour will be discussed in detail in the following Sections 2.2–2.4.

The background evolution is governed by the first Friedmann equation in the presence of dark matter and dark energy, $H^{2}\left(a\right)=H_{0}^{2}\left[\Omega_{\mathrm{m,0}}a^{-3}+\Omega_{\mathrm{de,0}}e^{-3\left(1+w\right)}\right]$, with $\Omega_{\mathrm{m,0}}$ and $\Omega_{\mathrm{de,0}}$ being the present-day mean matter and dark energy densities, and $H_{0}$ the Hubble constant. We neglect the impact of radiation at the late times relevant for this study, and treat $\Omega_{\mathrm{m}}$ as being made up of CDM + baryons, neglecting baryonic pressure for the collapsing objects. This corresponds to the standard treatment of the fluid equations when calculating collapse threshold and virial overdensity, as they are input for HMFs calibrated with dark-matter-only N-body simulations.

At late times, all scales relevant for structure formation are well within the horizon and we can safely take the Newtonian limit for the fluid equations. The continuity and Euler equation then read

[TABLE]

with density $\rho_{i}$, three-velocity $\boldsymbol{v}_{i}$ and pressure $p_{i}$ for each species $i$, and $\Phi$ denoting the Newtonian potential.

To investigate the non-linear evolution of density fluctuations we expand the fluid quantities, both for dark matter and dark energy (Pace et al., 2010, 2014). Fluctuations in density $\delta_{i}$ and pressure $\delta p_{i}$, as well as the peculiar velocity $\boldsymbol{u}_{i}$, are defined for each species $i$ through $\rho_{i}=\bar{\rho}_{i}\left(1+\delta_{i}\right)$, $p_{i}=\bar{p}_{i}+\delta p_{i}$, and $\boldsymbol{v}_{i}=a\left[H\left(a\right)\boldsymbol{x}+\boldsymbol{u}_{i}\right]$, respectively, where $\boldsymbol{x}$ is the comoving coordinate and overbars denote the corresponding background quantities. For fluids with constant sound speed $c_{\mathrm{s},i}$ and constant equation of state $w_{i}$, the corresponding equations for density perturbations and velocity divergence, $\theta_{i}=\left(\boldsymbol{\nabla}_{\boldsymbol{x}}\cdot\boldsymbol{u}\right)_{i}$, read

[TABLE]

The Poisson equation describes how the potential $\Phi$ is sourced by the density and pressure perturbations, that is

[TABLE]

where the sum runs over each species considered, here dark matter and dark energy. In the case of quasi-homogeneous dark energy with $c_{\mathrm{s}}=1$, the non-linear equations and their linearized counterparts are taken in the limit of negligible dark energy perturbations $\delta_{\mathrm{de}}\rightarrow 0$. This amounts to neglecting large-scale modes. For cold dark energy with $c_{\mathrm{s}}\rightarrow 0$ we have negligible dark energy pressure perturbations, $\delta p_{\mathrm{de}}\ll\delta\rho_{\mathrm{de}}$, since $\delta p_{\mathrm{de}}=c_{\mathrm{s}}^{2}\delta\rho_{\mathrm{de}}$.

For cold dark energy, the combination of the linearized version of equations (3)–(5) for dark matter and dark energy gives

[TABLE]

Initial conditions are chosen during matter domination when $\delta_{\mathrm{m}}\propto a$ holds. By solving equation (7) in this case, the initial dark energy contrast $\delta_{\mathrm{de},i}$ is obtained in terms of the dark matter density contrast $\delta_{\mathrm{m},i}$ as

[TABLE]

with $\dot{\delta}_{\mathrm{m,}i}=H\left(a_{i}\right)\delta_{\mathrm{m,}i}$ at initial scalefactor $a_{i}$, which we fix to $10^{-5}$. The initial dark matter density contrast $\delta_{\mathrm{m,}i}$ is adjusted around $\delta_{\mathrm{m},i}\sim 10^{-7}$ such that the point of collapse takes place at the low redshifts of interest, as we need to calculate the collapse threshold at the redshifts of our galaxy cluster data. This choice is degenerate with the initial scalefactor $a_{i}$ chosen. The point of collapse is defined as the singularity where the non-linear matter density perturbation diverges when evolving the set of coupled non-linear equations (3)–(5).

Following the non-linear density evolution corresponds to tracking the evolution of a spherical homogeneous tophat overdensity of radius $R$ until its radius reaches zero, i.e. the point of collapse (Gunn & Gott, 1972). We also show the radius evolution here, as we employ it to calculate the turnaround radius needed to connect to the radius at virialization through energy conservation (see Section 2.3). From Birkhoff’s theorem, following the radius evolution is equivalent to following the evolution of a separate closed FLRW universe where the scalefactor $a$ is replaced by a distinct scalefactor $R$ within the overdensity. The radius evolution of the spherical overdensity is obtained from the isotropic and homogeneous solution of the Euler equation (2) in the Hubble flow, $\boldsymbol{v}=H\boldsymbol{x}$, which corresponds to $\ddot{a}/a=-\boldsymbol{\nabla}\Phi$, with the scalefactor $a$ replaced by the radius $R$ and the gradient potential $\boldsymbol{\nabla}\Phi=\left(4\pi G/3\right)\sum\left(\rho_{i}+3p_{i}\right)\boldsymbol{x}$ inserted. Hence, the evolution of the spherical overdensity in the presence of cold dark energy is described by (Creminelli et al., 2010a)

[TABLE]

where we have used $\delta p_{\mathrm{de}}\approx 0$. For quasi-homogeneous dark energy with $c_{\mathrm{s}}=1$, we have $\rho_{\mathrm{de}}=\bar{\rho}_{de}$. For the evolution of the dark matter and dark energy densities within a spherical overdensity of radius $R$, the continuity equation (1) gives

[TABLE]

In order to solve for collapse, equations (9) and (10) are evolved until a singularity is reached. At the initial time $t_{i}$, we set the radius to $R_{i}=1$, the expansion rate in the linear regime during matter domination to $\left.\mathrm{d}R/\mathrm{d}t\right|_{i}=2\left(1-\delta_{\mathrm{m},i}/3\right)/(3t_{i})$, and employ equation (8) for the initial dark energy density contrast.

We checked that solving for the radius of the spherical overdensity using equations (9) and (10) gives the same time of collapse as obtained from the non-linear set of equations (3)–(5). In fact, these two approaches are equivalent. This can for example be shown for dynamical mutation of dark energy, where the effective equation of state within an overdensity or collapsed region is different from the background equation of state (Abramo et al., 2008). Thus, for the collapse overdensity we have $p_{\mathrm{c}}=w_{\mathrm{c}}\rho_{\mathrm{c}}$ and $\rho_{\mathrm{c}}=\bar{\rho}\left(1+\delta\right)$, for which the continuity equation (10) reads $\dot{\delta}+3\left(1+\delta\right)\left[h\left(1+w_{\mathrm{c}}\right)-H\left(1+w\right)\right]$, with $h=\dot{R}/R$. This is the same as equation (3) when expressing the equation of state inside the overdensity, $w_{\mathrm{c}}$, in terms of the background equation of state via $w_{\mathrm{c}}=w+\left(c_{\mathrm{s}}^{2}-w\right)\delta/\left(1+\delta\right)$.

Note that for dark energy sound speeds different from one or zero, the tophat profile evolution for density perturbations, with which equations (3)–(5) comply, does not hold. The absence of a sharp tophat profile leads to a scale- (or mass-) dependence in the perturbations, which propagates to derived quantities like the density threshold of collapse, so that either an interpolation down to the well-behaved case of sound speed zero or an averaging of the derived quantities, e.g. the threshold of collapse, over a tophat profile are necessary (Basse et al., 2011, 2014). This scale- (or mass-) dependence of the perturbations that leads to a breakdown of the tophat approximation can be translated into a Jeans mass, depending on the sound speed of the dark energy fluid. This mass corresponds to a characteristic scale where the effects of dark energy clustering become important. For example, sound speeds of the order of $10^{-4}$ and $10^{-5}$ correspond to masses of the order of $10^{14}M_{\odot}$ and $10^{15}M_{\odot}$, respectively, which are typical masses for galaxy clusters. In this work we opt for comparing cold dark energy of negligible sound speed against the standard quasi-homogeneous dark energy. The advantage is that these two limiting cases are fully consistent with the semi-analytical treatment described above, in which the tophat evolution of the spherical overdensity is physically motivated.

2.2 Collapse threshold

Here we discuss the cosmology-dependent linear density threshold of collapse for a range of values of the present-day matter density and dark energy equation of state. We include examples with $w$ values differing from -1, that are disfavoured by combined constraints with the CMB but allowed by cluster-only data, to display the qualitative trends of the characteristic quantities of collapsed structures.

To solve for the evolution of a spherical overdensity we evolve equations (3)–(5) until the non-linear density perturbation diverges and the point of collapse is reached. The linear density contrast of matter at the time of collapse is the so-called spherical collapse threshold $\delta_{\mathrm{c}}$. In an Einstein-de Sitter (EdS) universe (i.e. $\Omega_{\rm m}=1$), $\delta_{\mathrm{c}}=1.686$, independent of redshift and initial overdensity (see e.g. Peebles, 1980; Bertschinger & Jain, 1994). The inclusion of dark energy modifies the dynamics of the spherical collapse, introducing a redshift-dependence in the threshold of collapse. This redshift dependency of the collapse threshold is displayed in the top left panel of Fig. 1 for $w\in[-1.4,-0.6]$ (in steps of $0.2$ from top to bottom) for $c_{\mathrm{s}}^{2}=1$ and $c_{\mathrm{s}}^{2}=0$. Note that for $\Lambda$CDM, i.e. $w=-1$, these two cases coincide, as perturbations vanish for a cosmological constant. For both $c_{\mathrm{s}}=1$ and $c_{\mathrm{s}}=0$, $w>-1$ lowers $\delta_{\mathrm{c}}$ with respect to that of $\Lambda$CDM, while dark energy with $w<-1$ has the opposite effect. But, for cold dark energy the collapse threshold is lowered for $w<-1$ and enhanced for $w>-1$ as compared to $\Lambda$CDM, i.e., for both $w>-1$ and $w<-1$ the curves for $c_{\mathrm{s}}=0$ in comparison to $c_{\mathrm{s}}=1$ are closer to the $\Lambda$CDM case of $w=-1$.

As can be seen from the radius evolution of the spherical overdensity in the middle panels of Fig. 1, the change of collapse dynamics with $w$ translates into a slightly delayed collapse in a cold dark energy scenario as opposed to quasi-homogeneous dark energy for $w<-1$. This is because dark energy becomes important earlier, hindering the collapse, with $\rho_{\mathrm{de}}\sim\exp\left[-H\left(1+w\right)\right]$ growing for $w<-1$ with a positive argument in the exponential. The opposite is true for $w>-1$, where spherical overdensities collapse earlier, as dark energy starts to dominate later. In all cases, cold dark energy tends to bring radius evolution and collapse threshold closer to those of $\Lambda$CDM with respect to quasi-homogeneous dark energy. Also, the modification of the spherical collapse dynamics by the inclusion of dark energy becomes more apparent when further away from an EdS scenario, i.e. the lower $\Omega_{\mathrm{m}}$ the lower $\delta_{\mathrm{c}}$. Note that an increased $\Omega_{\mathrm{m}}$ is partly degenerate with shifting from $c_{\mathrm{s}}=1$ to $c_{\mathrm{s}}=0$ for $w>-1$ and vice versa for $w<-1$, which is to be expected for clustering dark energy behaving like some kind of extra matter contribution. This behaviour is shown in the top right panel of Fig. 1 for different values of $\Omega_{\mathrm{m}}<1$, with $w=-1$ fixed.

2.3 Virial overdensity

In the context of the Press–Schechter formalism (Press & Schechter, 1974) and its extensions (see e.g. Zentner (2007) for a review), a halo mass is defined as the mass enclosed by the virial radius $R_{\rm vir}$, within which there is an interior mean overdensity $\Delta_{\rm vir}$ with respect to a reference background density $\bar{\rho}_{\rm ref}$. In general, the virial overdensity $\Delta_{\rm vir}$ depends on redshift and cosmology. In contrast, observational masses are usually associated with more convenient fixed overdensities, with typical values such that the interior average density is $\rho_{\rm int}=300\bar{\rho}_{\rm m}$ or $\rho_{\rm int}=500\bar{\rho}_{\rm crit}$, where $\rho_{\rm crit}=3H_{0}^{2}/\left(8\pi G\right)$ is the critical density of the Universe. Thus, to convert between virial and observed masses we need to derive $\Delta_{\rm vir}$ explicitly for each cosmology.

In the spherical collapse model, the redshift of virialization is obtained by relating it to the turnaround of the collapsing sphere via the virial theorem. The time of turnaround is reached when the radius of the sphere is maximal, and equivalently the quantity $\log\left(\delta_{\rm NL}+1\right)/a^{3}$, with non-linear density contrast $\delta_{\rm NL}$, is minimal before diverging at the time of collapse. Radius and overdensity at turnaround can then be connected to radius and overdensity at virialization by making use of energy conservation. We use the connection between virial and turnaround radius given by Lahav et al. (1991). The virial radius gives the scalefactor of virialization, where we tracked the radius evolution of the sphere as described in Section 2. Having the time and radius of virialization, the virial overdensity $\Delta_{\mathrm{vir}}$ is due to mass conservation given by (Basilakos et al., 2009; Lee & Ng, 2010; Meyer et al., 2012)

[TABLE]

where $a_{\mathrm{ta}}$ and $R_{\mathrm{ta}}$ are the turnaround scalefactor and radius, and $a_{\rm vir}$ and $\delta_{\mathrm{NL,vir}}$ are the scalefactor and non-linear density contrast, respectively, at the time of virialization.

The virial overdensity needs to be calculated for every redshift and cosmological parameter set of interest. In order to speed up the calculations, we fit the virial density threshold on a grid of $\Omega_{\mathrm{m}}$ and $w$, aiming at subpercent accuracy. The fitting formula is an expansion around the EdS case of $\Omega_{\mathrm{m}}=1$ and constant $\Delta_{\mathrm{vir}}=18\pi^{2}$, with $\Omega_{\mathrm{m}}\left(z\right)$ order of unity at relevant redshifts. It is normalized for $\Delta_{\mathrm{vir}}$ to be inversely proportional to $\Omega_{\mathrm{m}}$. The fitting formula reads

[TABLE]

with parameters $a$, $b$ and the fractional matter density parameter $\Omega_{\mathrm{m}}\left(z\right)$. For our HMF calculations in Section 3, we interpolate the values of the virial density threshold on the fitted grid to convert observed cluster masses to virial masses.

The bottom panels of Fig. 1 show the virial overdensity for different values of $w$ and $\Omega_{\mathrm{m}}$. As was the case in Section 2.2 for the collapse threshold, the presence of dark energy tends to increase the virial density threshold compared to the constant EdS case. The effect is bigger the earlier dark energy becomes important. When comparing cold to quasi-homogeneous dark energy in the bottom left panel of Fig. 1, the virial overdensity is larger for $w<-1$, as collapse is hindered, and lower for $w>-1$, with cold dark energy helping the collapse in this case. For a cosmological constant, the bottom right panel of Fig. 1 makes clear again that lowering $\Omega_{\mathrm{m}}$ increases the virial overdensity, in accordance with the fitting formula used.

2.4 Dark energy mass contribution

In the case of cold dark energy, where dark energy of negligible sound speed is comoving with dark matter, dark energy perturbations can contribute to the total mass of the object. We therefore include the contributions from dark energy perturbations to the total mass by altering the predicted HMF. To estimate the extra dark energy contribution to the total mass, we calculate the dark energy mass at virialization as shown in Creminelli et al. (2010a). We assume a tophat profile and calculate the virial radius $R_{\mathrm{vir}}$ as in Section 2.3. We define the dark energy mass contribution at virialization as

[TABLE]

with dark energy overdensity $\delta\rho_{\mathrm{de}}$. The total halo mass $M$ is then rescaled by the dark energy mass contribution as $M\rightarrow M\left(1+\epsilon_{vir,de}\right)$, with the ratio between dark energy and dark matter mass at virialization $\epsilon_{\mathrm{vir}}\left(z\right)=M_{\mathrm{vir,de}}/M_{\mathrm{vir,m}}$, shown in the left-hand panel of Fig. 2 for different values of $w$. As expected, the extra dark energy mass contribution grows as $w$ deviates from the $\Lambda$CDM value of $w=-1$, i.e. with increasing inhomogeneity, and changes sign from a positive to a negative contribution when going from the non-phantom $w>-1$ to the phantom $w<-1$ side. The negative mass contribution for $w<-1$ corresponds to dark energy density escaping from dark matter potential wells, so that effectively the total density perturbations are smaller than without dark energy (Weller & Lewis, 2003). The right-hand panel of Fig. 2 shows the ratio $\epsilon_{\mathrm{vir}}\left(z\right)$ for different values of $\Omega_{\mathrm{m}}$ and fixed $w=-0.5$. For these models, the positive mass correction is higher the lower $\Omega_{\mathrm{m}}$, as the dark energy mass contribution increases compared to that of dark matter. Generally, the dark energy contribution also becomes more significant at lower redshifts for the same reason.

An alternative to integrating only over the perturbed dark energy density $\delta_{\rm de}$ is to rescale the total mass by integrating over a dark energy accretion rate for the halo, see again Creminelli et al. (2010a). As a conservative, lower estimate, throughout this paper we will use only the dark energy mass contribution from the dark energy overdensity as evolved at second order in perturbations.

3 Recalibrated Halo Mass Function

In this section, we show the impact of cold dark energy on the HMF, i.e., the comoving number density of haloes as a function of mass and redshift, by recalibrating the standard HMF for CDM to account for non-linear perturbative effects stemming from dark energy’s being cold (clustering). To do so we include the non-linear characteristic quantities derived in Section 2 into our mass function recalibration, as described below. For this, it is crucial to evaluate these quantities, i.e., $\delta_{\mathrm{c}}$, $\Delta_{\mathrm{vir}}$ and $\epsilon_{\mathrm{vir}}$, for each set of cosmological parameters when estimating and forecasting constraints in Sections 5 and 6, respectively.

In practice, we calculate both the Tinker HMF (Tinker et al., 2008) and the Sheth–Tormen HMF (Sheth & Tormen, 1999) for each set of cosmological parameters. We then form the ratio of the Sheth–Tormen HMF of cold dark energy to the Sheth–Tormen HMF of quasi-homogeneous dark energy and multiply it by the standard Tinker HMF. The idea here is that the ratio of Sheth–Tormen HMFs encapsulates the difference between cold and quasi-homogeneous dark energy, i.e., the impact of cold dark energy on the HMF, while the Tinker mass function provides us with a precise estimate for the standard cosmological mass function, informed by N-body simulations. The ratio of Sheth–Tormen HMFs then includes both effects at the linear level (which is also the case for the Tinker HMF) as well as second-order effects introduced via spherical collapse quantities, which might, due to their scale-dependence, enable us to distinguish the two models.

The recalibrated HMF for the expected number of haloes of mass $M$ at redshift $z$ in a cold dark energy scenario reads

[TABLE]

with ‘ST’ designating a Sheth–Tormen and ‘T’ a Tinker HMF. A similar ratio has been employed in Sartoris et al. (2010b), see also Shandera et al. (2013b) and Cataneo et al. (2015) in order to constrain primordial non-Gaussianity, in the former, as first prescribed in LoVerde et al. (2008), and to distinguish $f\left(R\right)$ modified gravity theories from GR+$\Lambda$CDM, in the latter.

As a reminder, the functional form of an HMF under spheroidal collapse with mass-independent collapse threshold reads

[TABLE]

The Tinker HMF is defined by the multiplicity function $f(\sigma)$ (Tinker et al., 2008), which describes non-linear clustering behaviour as a function of matter variance,

[TABLE]

with the parameters $A$, $a$, $b$ and $c$ fitted using CDM simulations. These parameters depend on the definition of cluster mass as multiples of the overdensity $\Delta$ with respect to $\rho_{\mathrm{crit}}$. The variance $\sigma\left(R\right)$ of the density field smoothed at radius $R$ is defined as

[TABLE]

This integrates the matter power spectrum $P\left(k\right)$ over the tophat window function $W\left(kR\right)$ for a sphere of radius $R$ in Fourier space. The theoretical uncertainty of the Tinker HMF from simulations is $\leq 5\%$. We include this uncertainty in our cluster analysis by increasing the width of the Gaussian priors on the parameters of the HMF to $10\%$, as well as accounting for their covariance.

The Sheth–Tormen HMF accounts for ellipsoidal instead of spherical collapse, giving an improved fit to numerical simulations. The functional form for an HMF under ellipsoidal collapse with mass-dependent peak height $\nu$ states

[TABLE]

$f\left(\nu\right)$ depends now on the peak height $\nu=\delta_{\mathrm{c}}/\sigma$, with collapse threshold $\delta_{\mathrm{c}}$, and not solely on $\sigma$ as in the case of spheroidal collapse. As can be seen in equation (18), this HMF depends on the cosmological background via the mean matter density $\bar{\rho}_{m}$, as well as on linear density perturbations through $\sigma$ and the linear growth function $D\left(z\right)$. The multiplicity function for the Sheth-Tormen HMF reads (Sheth & Tormen, 1999)

[TABLE]

with parameters $a$ and $p$ fitted to simulations, where $a$ is determined by the number of massive haloes, $p$ by the shape of the mass function at the low-mass end and $A$ is the normalization ensuring the integral of $f\left(\nu\right)$ over all $\nu$ gives unity.

When comparing the HMF to observations, the non-linear information on the properties of collapsed structures can be included via $\delta_{\mathrm{c}}$, $\Delta_{\mathrm{vir}}$ and $\epsilon_{\mathrm{vir}}$. Generally, these quantities are redshift-dependent, and change with cosmology, and therefore depend on the dark energy sound speed. However, when comparing structure formation data to theoretical expectations, the collapse threshold and virial overdensity are sometimes taken at their constant EdS values, or as fitting functions for a $\Lambda$CDM cosmology (see e.g. Kitayama & Suto (1996) and Nakamura & Suto (1997)).

We go beyond this practice and advocate a more accurate calibration of the model-dependent non-linearities in the HMF by using the results obtained within the spherical collapse model of Sections 2.2–2.4. We insert $\delta_{\mathrm{c}}\left(z\right)$ for a given cosmology into the Sheth-Tormen HMF of equation (18), both for quasi-homogeneous and cold dark energy, and form the ratio in equation (14). We convert between observed cluster and virial masses, defined using $\Delta_{\mathrm{vir}}$ from Section 2.3, assuming a Navarro–Frenk–White (NFW) profile (Navarro et al., 1996; Hu & Kravtsov, 2003). We also include the mass correction at virialization due to dark energy calculated as in Section 2.4. This accounts for differences in virialization, as well as for the mass shift introduced in the HMF due to cold dark energy.

The standard use of the NFW profile is motivated by N-body simulations incorporating quasi-homogeneous dark energy in the background, and CDM particles. We expect this profile to also hold reasonably well in the case of cold dark energy, since for negligible sound speed dark energy will comove with dark matter and therefore should not significantly alter the radial profile, while adding more CDM-like substance to the flow. However, this assumption should be checked in detail with future N-body simulations including cold dark energy. As we account for here, cold dark energy has an impact on collapsing objects, for example changing the typical collapse overdensities, affecting the number of objects that collapsed in the first place.

In the top left panel of Fig. 3, we show the recalibrated mass function of equation (14) as a function of virial mass for both $c_{\mathrm{s}}=1$ and $c_{\mathrm{s}}=0$, and $w=-1$ ($\Lambda$CDM), $w<-1$ and $w>-1$, with $\Omega_{\mathrm{m}}=0.3$ fixed for all curves. The predicted cluster abundances are systematically lower for cold compared to quasi-homogeneous dark energy when $w<-1$, and systematically higher for $w>-1$. For $w=-1$ both curves coincide, as dark energy perturbations vanish. This behaviour is to be expected, since the contribution from dark energy perturbations suppresses structure formation for $w<-1$ and enhances it for $w>-1$ with respect to the quasi-homogeneous case. The bottom left panel of Fig. 3 underlines these sizeable deviations, showing the ratio of the Sheth–Tormen cluster mass functions of $c_{\mathrm{s}}=0$ to $c_{\mathrm{s}}=1$ dark energy for different values of $w$ and $\Omega_{\mathrm{m}}$. We recalibrate the HMF with this ratio.

The top right panel of Fig. 3 shows the ratio of the recalibrated cluster mass functions for cold and quasi-homogeneous dark energy to the mass function in a $\Lambda$CDM universe for the same set of parameters as in the top left panel. For a cold dark energy scenario, this ratio indicates deviations from $\Lambda$CDM by up to a number of times, depending on parameters chosen. The effect is especially pronounced at the high-mass end where massive clusters of galaxies reside.

The bottom right panel of Fig. 3 shows the dependence on $\Omega_{\mathrm{m}}$ of the recalibrated cluster mass function, for both $c_{\mathrm{s}}=1$ and $c_{\mathrm{s}}=0$ dark energy. We display curves for $\Omega_{\mathrm{m}}=0.2,0.3,0.4$ at fixed $w=-0.7$. As expected, the predicted cluster abundances are systematically higher for cold compared to quasi-homogeneous dark energy as it was in the corresponding case in the bottom left panel for a $w>-1$. Note that the difference between cold and quasi-homogeneous dark energy is more pronounced for low matter densities, as the relative importance of dark energy grows.

Even though here we focused on displaying the parameter dependence of the mass function recalibration on $\Omega_{\mathrm{m}}$, $w$ and $c_{\mathrm{s}}$, we stress that the mass function estimate to be obtained with our Markov Chain Monte Carlo (MCMC) data analysis depends on a multitude of parameters, such as those of the mass-observable relations fitted simultaneously with cosmology. Dependencies on astrophysical parameters of course interplay with the ones on cosmological parameters displayed in this section, while our recalibration remains important in order to incorporate the effect of dark energy clustering. In Section 5 we will use the recalibrated mass function presented here to constrain cosmological parameters, while marginalizing over astrophysical parameters, with measurements of cluster number counts.

4 Data

For the parameter constraints in Section 5 we tested both a cluster-only data set and a combination of cluster, CMB, baryon acoustic oscillation (BAO) and type Ia supernova (SN Ia) data. The cluster data used here consist of three X-ray samples of clusters taken from the ROSAT All-Sky Survey (RASS; Truemper, 1993), the Brightest Cluster Sample (BCS; Ebeling et al., 1998), the ROSAT-ESO Flux Limited X-ray sample (REFLEX; Boehringer et al., 2004) and the bright sample of the Massive Cluster Survey (MACS; Ebeling et al., 2010), together with follow-up X-ray and optical imaging data.111Note that we do not include measurements of the gas mass fraction $f_{\rm gas}$ (Mantz et al., 2014), as the relation between X-ray gas mass, total cluster mass and baryonic fraction has not been investigated yet for cold dark energy cosmology. To select for high-mass clusters and avoid incompleteness, a flux limit of 0.1–2.4 keV for luminosities $>2.5\times 10^{44}h^{-2}$ erg s*-1* was applied; furthermore systems with X-ray emission dominated by active galactic nuclei have been removed. After cuts, the sample contains 224 clusters. For 94 clusters of the sample, follow-up ROSAT and/or Chandra data were used to derive X-ray luminosities, temperatures and gas masses in order to constrain cluster scaling relations.222X-ray luminosities, temperatures and gas masses are direct measurements of intrinsic quantities that we utilize to fit cluster scaling relations simultaneously with cosmology, as described in detail in Mantz et al. (2010b). The absolute cluster mass scale is calibrated using weak-lensing measurements for 50 massive clusters taken from the Weighing the Giants program introduced in von der Linden et al. (2014), Applegate et al. (2014) and Kelly et al. (2014).

When using our cluster data set alone, consisting of RASS catalogues, and follow-up X-ray and lensing data, we also include Gaussian priors on the cosmic baryon density $100\,\Omega_{\mathrm{b}}h^{2}=2.202\pm 0.045$  (Cooke et al., 2014) and the present-day Hubble parameter $h=0.738\pm 0.024$ (Riess et al., 2011). For the full combination of clusters+CMB+BAO+SNIa, priors on $h$ and $\Omega_{\mathrm{b}}$ are not required. For the CMB data, we combine Planck (1-year release plus WMAP polarisation data; Planck Collaboration et al., 2014a), together with Atacama Cosmology Telescope (ACT; Das et al., 2014) and South Pole Telescope (SPT; Keisler et al., 2011; Reichardt et al., 2012; Story et al., 2013) measurements at higher multipoles. We also use the Union2.1 compilation of SNIa data (Suzuki et al., 2012), as well as BAO data from the 6-degree Field galaxy Survey (6dF; Beutler et al., 2011) and the Sloan Digital Sky Survey (SDSS release 7 and 9, see Padmanabhan et al., 2012; Anderson et al., 2014).

In the following, we employ our cluster-only data set and its combination with the above complementary data sets to constrain cosmological parameters for both cold and quasi-homogeneous dark energy, using the recalibrated mass function introduced in the previous section.

5 Parameter estimation

Here, we present the first observational constraints on a dark energy model with constant $w$ and $c_{\mathrm{s}}=0$. For this, we employ our cluster growth data and additional complementary data sets, as described in the previous section. Our likelihood analysis and parameter estimation via an MCMC sampler includes theoretical and experimental scatter in the HMF and the mass-observable scaling relations, and also incorporates modified characteristic quantities for structure formation: collapse threshold (Section 2.2), virial overdensity (Section 2.3) and cold dark energy mass contribution at virialization (Section 2.4). These quantities are interpolated on a grid for computational speed, aiming at per cent accuracy.

We proceed as follows. For each trial cosmology, the ratio of the Sheth–Tormen HMF for cold versus quasi-homogeneous dark energy (equation 14) is formed, depending both on mass and redshift. In this way, the cosmology-dependent collapse threshold is included into the Sheth–Tormen HMF, on top of the background and linear level calculations taken from CAMB (Lewis et al., 2000; Howlett et al., 2012).333Version March 2013; http://camb.info/ In addition, the cold dark energy mass contribution is taken into account, resulting in a shift of the mass scale, i.e., adding mass for $w>-1$ and reducing mass for $w<-1$ (see Section 2.4). Equation (14) is then used as our recalibrated cold dark matter plus cold dark energy mass function to compare with observed cluster number counts. We convert between different mass definitions at differing overdensities assuming an NFW profile (see appendix in Hu & Kravtsov, 2003). Virial masses are converted to other mass definitions using $\Delta_{\mathrm{vir}}$, as derived in Section 2.3, which accounts for differences in virialization between the cold and quasi-homogeneous dark energy models.

We use a modified version of CosmoMC (Lewis & Bridle, 2002; Lewis, 2013),444Version October 2013; http://cosmologist.info/cosmomc/ which incorporates a module that evaluates our cluster growth likelihood (Mantz et al., 2015) using the data described in Section 4. We also employ uniform priors on $\Omega_{\mathrm{m}}h^{2}\in\left[0.025,0.3\right]$, $\Omega_{\mathrm{b}}h^{2}\in\left[0.005,0.1\right]$ and $w\in\left[-1.5,-0.5\right]$ in order to keep the spherical collapse calculations valid and numerically stable. Throughout our analysis, we assume a spatially flat cosmology, the standard effective number of relativistic species $N_{\mathrm{eff}}=3.046$ and the minimal species-summed neutrino mass $\sum m_{\mathrm{\nu}}=0.056$.555Note that the chosen assumption on the sum of neutrino masses can affect the high-mass end of the cluster mass function. In order to separate out the effect of dark energy clustering, we relied for now on fixing the effective sum of neutrino masses.

Fig. 4 shows our constraints on $\Omega_{\mathrm{m}}$, $w$ and $\sigma_{8}$, the rms variation of the matter density field within 8 $h^{-1}$Mpc spheres, using the combination of cluster growth data with CMB, BAO and SN Ia data, as described in Section 4. The blue contours show the results for the standard mass function analysis of quasi-homogeneous dark energy with sound speed $c_{\mathrm{s}}=1$, employing the Tinker HMF, and the magenta contours those for cold dark energy with sound speed $c_{\mathrm{s}}=0$. For the latter, we employ the recalibrated cluster mass function of Section 3 that we implemented into our cluster growth likelihood analysis alongside the standard analysis, which includes simultaneous fitting of cluster scaling relations. For the combination of data sets, a slight shift between the confidence contours for cold versus quasi-homogeneous dark energy is observed, suggesting that this effect will be important to consider with upcoming, more precise measurements. Note this shift puts the $w$ value for $c_{\mathrm{s}}=0$ closer to the $\Lambda$CDM case, as can be expected from the predicted HMFs relative to $\Lambda$CDM, for example in the top right panel of Fig. 3.

In Table 1, we show the corresponding marginalized best-fitting values and 68.3 per cent confidence intervals for $\sigma_{8}$, $\Omega_{\mathrm{m}}$ and $w$. We find that for the full combination of cluster+CMB+BAO+SNIa data, differences in the marginalized best-fit values at the 1–2$\,\sigma$-level appear between cold and quasi-homogeneous dark energy. The slight preference for low $w$ values is comparable to that found for other similar data combinations, like Planck+SNIa+BAO, with best-fitting values of $w\sim-1.1$ for $w$CDM.666https://wiki.cosmos.esa.int/planckpla/index.php/

Cosmological$\_$Parameters (Planck Collaboration XVI 2014) Using cluster growth data alone, we find that the derived best-fitting values for cold and quasi-homogeneous dark energy agree within their respective 68.3 per cent confidence intervals, and the estimates obtained for $w$ agree with $\Lambda$CDM at the 1$\sigma$ level.

6 Fisher forecast

In this section, we use a Fisher matrix formalism to forecast constraints on cosmological parameters from future cluster number counts for the case of both cold and quasi-homogeneous dark energy, in order to test for the impact of our model choice on parameter inference. For this analysis we use the cosmological parameter vector ${\bf p}=\left\{\Omega_{\mathrm{m}},\Omega_{\mathrm{b}},h,n_{\mathrm{s}},\sigma_{\mathrm{8}},w\right\}$, of which we vary only three parameters: $\Omega_{\mathrm{m}}$, $\sigma_{\mathrm{8}}$ and $w$. As fiducial values, we take those of the Planck best fit for the base case $w$CDM lowl+highL+SNLS in Planck Collaboration et al. (2014b),777PLA/base$\_$w$\_$planck$\_$lowl$\_$lowLike$\_$highL$\_$post$\_$SNLS*,*888https://wiki.cosmos.esa.int/planckpla/index.php/File:Grid$\_$limit68.pdf being $\Omega_{\mathrm{m}}=0.28$, $\Omega_{\mathrm{b}}=0.044$, $h=0.709$, $n_{\mathrm{s}}=0.9581$, $\sigma_{\mathrm{8}}=0.87$ and $w=-1.124$. Note also that we do not include a galaxy cluster power spectrum analysis.

The Fisher matrix for cluster number counts, with $N_{l,m}$ number of clusters in the $l$th redshift bin and $m$th observed mass bin, reads (Holder et al., 2001)

[TABLE]

where the inverse covariance is given by $C_{l,m}^{-1}=1/N_{l,m}$. The $N_{l,m}$ expected for a survey with sky coverage $\Delta\Omega$ is (Majumdar & Mohr, 2003)

[TABLE]

with comoving volume element per unit redshift interval $\mathrm{d}V/\mathrm{d}z$, HMF $n\left(M,z\right)$ and probability to assign an observed mass $M^{\mathrm{ob}}$ to a cluster of true mass $M$, $p\left(M^{\mathrm{ob}}|M\right)$. The cosmology-dependent comoving volume element is given by

[TABLE]

with angular diameter distance $d_{\mathrm{A}}\left(z\right)$ and Hubble rate $H\left(z\right)$. For quasi-homogeneous dark energy with sound speed $c_{\mathrm{s}}=1$ we use the HMF $n\left(M,z\right)$ from Tinker et al. (2008) at an overdensity of $\Delta=200$ with respect to the background density. For cold dark energy with sound speed $c_{\mathrm{s}}=0$, we use instead the recalibrated HMF of equation (14) with cosmology-dependent $\delta_{\mathrm{c}}$ (Section 2.2) and $\epsilon_{\mathrm{vir}}$ (Section 2.4). It is worth noting here that $\epsilon_{\mathrm{vir}}$ effectively introduces a mass bias between true and observed mass. This is due to the true mass $M$ being effectively shifted to $M\left(1+\epsilon_{\mathrm{vir}}\right)$ by the presence of clustering dark energy.

Following Lima & Hu (2005) and assuming a lognormal distribution with variance $\sigma^{2}_{\ln M}$ we have

[TABLE]

where $x\left(M^{\mathrm{ob}},M,z\right)$ relates true and observed mass as

[TABLE]

Inserting equation (23) into equation (21) and integrating over $\mathrm{d}M^{\mathrm{ob}}$, we have

[TABLE]

with $x_{m}=x\left(M^{\mathrm{ob}}_{l,m},M,z_{l}\right)$ and the complementary error function $\mathrm{erfc}\left(x\right)$.999The complementary error function is defined as

$\mathrm{erfc}\left(x\right)=1-\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-t^{2}}\mathrm{d}t$.

As in Sartoris et al. (2010b) and Sartoris et al. (2016), we model the mass-observable relation using the bias between observed and true masses $\ln M_{\mathrm{bias}}$ and the intrinsic scatter $\sigma^{2}_{\ln M}$ as

[TABLE]

with fiducial parameters

[TABLE]

For the quasi-homogeneous dark energy model, this choice of fiducial parameters corresponds to no bias in the mass-observable scaling relation. For cold dark energy, the shift between true and observed mass introduced by $\epsilon_{\mathrm{vir}}$ translates into a mass- and redshift-dependent mass bias.

We select this form for the mass-observable scaling relation in equation (27), including a mass scatter of $\sigma_{\ln M}=0.4$, to be consistent with the observed mass-richness relation for massive, X-ray-selected clusters at $z<0.5$ (Mantz et al., 2016). The richness measurements for that relation were drawn from SDSS data using the redMaPPer algorithm (Rykoff et al., 2014; Rozo et al., 2015), which is a red sequence cluster finder designed to handle an arbitrary photometric galaxy catalogue with excellent photometric redshift performance, completeness and purity, and which is currently in use for Dark Energy Survey (DES) cluster studies.101010http://www.darkenergysurvey.org We caution, however, that the scatter chosen may be optimistic when applied to optically selected clusters, as in some cases projection effects can boost the measured richness. Also, the redshift evolution of the scatter is not well known and the details will depend on the selection algorithm used. We therefore assume that it remains constant, even though the validity of this assumption may deteriorate at lower masses and higher redshifts, because there the impact of projection effects is expected to increase.

In our Fisher matrix forecast we fix the nuisance parameters $\alpha$, $\sigma_{\ln M,0}$ and $\beta$ to their fiducial values and account for an expected uncertainty in the absolute mass calibration at $\pm 3\%$. The cluster mass threshold of detection $M_{\mathrm{thr}}$ is approximated as independent of redshift, although in practice $M_{\mathrm{thr}}$ becomes higher at higher redshifts depending on the survey-specific selection function. For DES, we use a sky coverage of $5000\,$deg2 from $z=0.1$ out to $z=1.0$ with $\Delta z=0.05$ to match the expected photo-$z$ error for masses larger than $M_{\mathrm{thr}}\approx 10^{14}M_{\odot}$ (Abbott et al., 2016; Rykoff et al., 2016). We bin the mass in steps of $\Delta\log\left(M^{\mathrm{ob}}/M_{\odot}\right)=0.2$ from the threshold mass $M_{\mathrm{thr}}$ up to a maximum of $\log\left(M^{\mathrm{ob}}/M_{\odot}\right)\leq 16$. The survey specifications used for the Fisher matrix forecasts performed in this section are summarized in Table 2. For this analysis, we have implemented our HMF recalibration of Section 3 and the Fisher matrix calculation for cluster number counts sketched above into the publicly available CosmoSIS code (Zuntz et al., 2015).111111Version 1.4; https://bitbucket.org/joezuntz/cosmosis/wiki Within CosmoSIS we ran the provided Fisher matrix sampler together with CAMB,121212Version January 2015. making use of the existing routine for the calculation of the Tinker and Sheth-Tormen HMFs.131313http://wwwmpa.mpa-garching.mpg.de/$\sim$komatsu/crl/

Fig. 5 shows the redshift-binned number density of clusters to be expected for DES in the redshift range $0.1<z<1.0$ for both cold (magenta) and quasi-homogeneous (blue) dark energy. The inset panel shows the corresponding ratio of $c_{s}=1$ to $c_{s}=0$, with about $3\%$ more clusters detected in the $z=0.1-0.2$ bin in the $c_{s}=1$ model. In Fig. 6, we show the resulting marginalized forecasted constraints for cold (dashed, magenta contours) versus quasi-homogeneous (solid, blue contours) dark energy. We find that constraints on $\Omega_{\mathrm{m}}$ and $w$ are tighter for cold dark energy. This happens though at the expense of those on $\sigma_{8}$, with less degeneracy between this parameter and $\left(\Omega_{\mathrm{m}},w\right)$ and the degeneracy direction tilting accordingly. For physical intuition, compare the ratio between Sheth–Tormen HMFs for the two sound speeds in the bottom left panel of Fig. 3. This ratio exhibits a strong scale-dependence due to dark energy being clustering, and this dependency is sensitive to both $\Omega_{\mathrm{m}}$ and $w$, which results in tighter constraints on these parameters. In general, however, the HMF in the $c_{\mathrm{s}}=0$ case varies less e.g. when $w$ is varied with respect to $\Lambda$CDM, see for example the top left panel of Fig. 3. This leads to cold dark energy being less constrained at the linear level than quasi-homogeneous dark energy and therefore the uncertainty in $\sigma_{8}$ being larger, and the degeneracy of $\sigma_{8}$ with other parameters being smaller. The differences in marginalized confidence intervals are summed up in Table 3, together with a Figure of Merit (FoM), defined similarly to that of the Dark Energy Task Force (DETF – Albrecht et al., 2006) as the square root of the determinant of the inverse Fisher matrix, $\mathrm{FoM}=1/\sqrt{|F_{i}^{-1}|}$, in our case in the $\left(\Omega_{m},\,w\right)$-plane. We obtain constraints of the order of $\approx 10^{-3}$ for $\Omega_{m}$ and $\sigma_{8}$, and $\approx 10^{-2}$ for $w$, as well as an FoM $\approx 10^{4}$, with constraints on $\sigma_{8}$ being about half and those on $\Omega_{\mathrm{m}}$ and $w$ being about a factor of 1.5–2 tighter for cold dark energy.

Note that for Euclid141414http://www.euclid-ec.org and Large Synoptic Survey Telescope (LSST)151515https://www.lsst.org/lsst$\_$home.shtml data, with higher sky coverages and a larger redshift range, as well as better mass calibration, constraints on cosmological parameters are expected to be stronger. Since both LSST and Euclid extend to higher redshifts than DES, investigating the impact of cold dark energy with a varying dark energy equation of state should be particularly interesting. Due to the uncertain observable-mass relations at higher redshifts, we postpone such analyses for later studies.

7 Conclusions

Cold dark energy of sound speed zero presents an interesting phenomenology as compared to standard quasi-homogeneous dark energy with sound speed one. It for example adds on top of the clustering of matter an extra clustering component due to dark energy perturbations that renders the model potentially distinguishable. Within the semi-analytical spherical collapse model, in this paper we derived the non-linear characteristic quantities required for the recalibration of the cluster mass function for cold dark energy. These are the collapse threshold, the virial overdensity and a dark energy mass contribution for cold dark energy. We incorporated these quantities into the HMF, encoding non-linear cosmological model information, to obtain a new cold dark matter plus cold dark energy mass function.

We use this new recalibrated mass function to constrain cosmological parameters for both quasi-homogeneous dark energy and cold dark energy with current state-of-the-art cluster growth data, together with a combination of standard cosmological data sets. For this combined data set a shift in the best-fitting parameter values of up to $1\sigma$ can be detected, with for example $\sigma_{8}=0.806\pm 0.014$ for sound speed zero and $\sigma_{8}=0.823\pm 0.017$ for sound speed one. These results and the comparison of our recalibrated mass functions for both models makes clear that including further non-linear model information has the potential to distinguish cold dark energy from the standard quasi-homogeneous case, as well as to ease parameter tensions present in the literature (see e.g. Köhlinger et al. (2017)).

In order to predict the ability to distinguish cold versus quasi-homogeneous dark energy with upcoming cluster growth data, we made a Fisher matrix forecast of cosmological parameter constraints for the ongoing DES cluster survey, and found $\Omega_{m}$ and $\sigma_{8}$ constraints of the order of $10^{-3}$, and of $10^{-2}$ for $w$. Our results also show tighter constraints on $\Omega_{\mathrm{m}}$ and $w$, increasing them by about 50$\%$ to 100$\%$ for cold dark energy, while the contrary is true for $\sigma_{8}$, as well as significant differences in the parameter degeneracies between the cold and quasi-homogeneous dark energy models.

More and better data, as well as combinations with other data, should enhance the differences in the estimated parameters for cold versus quasi-homogeneous dark energy. An interesting direction for further studies would be a more realistic treatment that either allows both the dark energy sound speed and equation of state to vary as free parameters, or employs models for which these parameters naturally evolve with redshift, for example within a Parametrized Post-Friedmann framework (Hu & Sawicki, 2007; Hu, 2008; Fang et al., 2008). Here, using structure formation data in the non-linear regime sensitive to scale-dependent behaviour, which occurs for sound speeds different from zero or one, will prove valuable.

Acknowledgements

The computational analysis was performed using the High Performance Computing facility at the University of Copenhagen (HPC@UCPH), and the Gardar supercomputer of the Nordic HPC project. This work was partially supported by the SFB-Transregio TR33 ‘The Dark Universe’ and the DNRF. DR is supported by a NASA Postdoctoral Program Senior Fellowship at the NASA Ames Research Center, administered by the Universities Space Research Association under contract with NASA. MC acknowledges support from the European Research Council under grant number 647112. We acknowledge support from the U.S. Department of Energy under contract number DE-AC02-76SF00515, and from the National Aeronautics and Space Administration under grant number NNX15AE12G issued through the ROSES 2014 Astrophysics Data Analysis Program.