Clustering dark energy and halo abundances

Within the standard paradigm, dark energy is taken as a homogeneous fluid that drives the accelerated expansion of the universe and does not contribute to the mass of collapsed objects such as galaxies and galaxy clusters. The abundance of galaxy clusters -- measured through a variety of channels -- has been extensively used to constrain the normalization of the power spectrum: it is an important probe as it allows us to test if the standard $\Lambda$CDM model can indeed accurately describe the evolution of structures across billions of years. It is then quite significant that the Planck satellite has detected, via the Sunyaev-Zel'dovich effect, less clusters than expected according to the primary CMB anisotropies. One of the simplest generalizations that could reconcile these observations is to consider models in which dark energy is allowed to cluster, i.e., allowing its sound speed to vary. In this case, however, the standard methods to compute the abundance of galaxy clusters need to be adapted to account for the contributions of dark energy. In particular, we examine the case of clustering dark energy -- a dark energy fluid with negligible sound speed -- with a redshift-dependent equation of state. We carefully study how the halo mass function is modified in this scenario, highlighting corrections that have not been considered before in the literature. We address modifications in the growth function, collapse threshold, virialization densities and also changes in the comoving scale of collapse and mass function normalization. Our results show that clustering dark energy can impact halo abundances at the level of 10\%--30\%, depending on the halo mass, and that cluster counts are modified by about 30\% at a redshift of unity.

However, there exists many DE models which can be described by a perfect fluid with varying sound speed. Models that show this behavior include the tachyon field [8,9] and, in general, the whole class of minimally coupled K-essence scalar fields [10,11], in which the sound speed can be freely chosen. Models with negligible sound speed can also be constructed with two scalar fields [12] and effective field theory [13].
When DE has negligible sound speed, its perturbations have no pressure support and can grow at the same pace as matter perturbations. Therefore, they can evolve nonlinearly and impact the formation of galaxy clusters. As cosmological simulations that include DE fluctuations have not yet been conducted, one has to resort to semi-analytical approaches such the spherical collapse model (SCM) [14] together with Press-Schechter (PS) or Sheth-Tormen (ST) mass functions [15,16] in order to explore the impact of clustering DE in the abundances of halos. In fact, the SCM was generalized to study clustering DE and it was shown that DE fluctuations can indeed become nonlinear, impact matter growth and the critical density threshold for collapse, δ c , and contribute to the mass of the clusters [17][18][19][20][21][22][23][24][25]. The impact of clustering DE was also studied in higher order perturbation theory [26][27][28].
In this work we carefully examine how the halo mass function is modified in the presence of clustering DE -a dark energy fluid with negligible sound speed -with a redshift-dependent equation of state, highlighting corrections that have not been considered in the literature before. We reanalyze some modifications that were already considered in the literature, such as the growth function, collapse threshold and virialization densities. Moreover, we also consider two new modifications: the comoving scale of collapse and mass function normalization. Section 2 is devoted to the study of the SCM in the case of clustering DE, while Section 3 discusses how the halo mass function is modified in order to take into account the contribution of DE fluctuations to the halo virialization. In Section 4 we show that clustering dark energy can impact halo abundances at the level of 10%-30%, depending on the halo mass, and that cluster counts are modified by about 30% at a redshift of unity. We draw our conclusions in Section 5.

Spherical collapse in the fluid and radius approach
The SC model can be generalized to treat fluids other than pressureless matter using the Pseudo-Newtonian Cosmology approach [19][20][21]. This framework is particularly useful to study clustering DE (c s = 0), in which case the equations are: where θ is the divergence of the peculiar velocity, δ m = δρ m /ρ m and δ de = δρ de /ρ de . In this approach, the redshift of collapse, z c , is determined when δ m → ∞, which is numerically implemented as certain threshold that reproduces the EdS results. As we will show, although this method is widely used to study the nonlinear evolution, the determination of the critical threshold might not be the best approach when DE fluctuations are present.
In order to understand the issue, let us now consider the radial approach. The original SC model, which determines the evolution of a spherical shell of physical radius R, only considers pressureless matter and can be solved analytically. When DE with negligible c s is present, the dynamical equation for R can be written as We also need equations for the contrasts of matter and DE. For matter, we assume that the total mass inside the shell with radius R is conserved, then we havė Comparing with Eq. (2.1), we identify Note that deriving Eq. (2.6) with respect to time and using Eq. (2.4), we recover Eq. (2.3). The radial approach itself can not determine the evolution of DE fluctuations -note that the evolution of matter fluctuations is given by the assumption that the total matter within R is conserved, which is not valid for DE. It is possible to derive an equation for DE fluctuations using the relation between θ and R, and Eq. (2.2), which yieldṡ Hence, the system of Eqs. (2.4), (2.5) and (2.7) determines the evolution of the shell radius, matter and DE fluctuations. The dynamical evolution in the two approaches is identical, however, as we will show, the commonly used collapse criteria in the two approaches are not equivalent.
In the radius approach, the redshift of collapse, z c , is given when R → 0. In the EdS model, the RHS of Eq. (2.4) only has matter quantities and the critical density of collapse takes the standard value δ c = δ lin m (z c ) 1.686. In the presence of DE, however, Eq. (2.4) shows that DE fluctuations also contribute to the collapse, hence the density contrast that gives the threshold density contrast has to be modified to take this new contribution into account. Therefore, instead of determining the collapse threshold only with the matter contribution, the natural quantity to use is the weighted total fluctuation: At high-z, when DE is subdominant in the background, the contribution of DE fluctuations to δ tot is very small, regardless of the magnitude of δ de . Conversely, at low-z, DE dominates the background and even relatively small DE fluctuations impact δ tot . So, the definition of δ tot reflects the fact the gravitational potential, which drives the evolution of R, is sourced by densities fluctuations,ρδ. The same quantify was used to study higher order perturbation theory in clustering quintessence models [26].
The critical contrast is then given by the linearly evolved total fluctuation: where z c is the redshift at which δ tot is above a certain numerical threshold. In the same fashion, the normalized growth function is given by .
Another interesting point about the relation between these two approaches is the error in δ c when using the fluid approach, as reported in [29] (see also [30]). The authors show that the numerical collapse threshold is not the same for different redshifts, which in turn generates a spurious increase of δ c with z. This can be corrected by calibrating the numerical threshold as a function of the redshift in order to reproduce some known δ c and demanding that it approaches the EdS value at high-z. Once this calibration is done, we verified that both approaches give essentially the same results for the models under consideration in this work. However, the radial approach is more robust because it does not demand any redshiftdependent calibration and for this reason it will be used here.

Background evolution and growth function
In order to show the impact of DE fluctuations on the halo collapse, we adopt a DE with a linear parametrization of the equation of state, w = w 0 + (1 − a) w a . In this work we do not consider models with w < −1 because in such case, during the nonlinear regime, the dynamical evolution can generate δ de < −1, which implies a non-physical negative energy density. In fact, we observed this happens around virialization densities for 1 + w −0.1. Hence we prefer to focus on the less problematic region of parameter space of non-phantom equations of state in order to clearly show the modifications in the mass function.
The latter happens because in Eq. (2.2) the coupling between the density contrast and gravitational force is proportional to (1 + w + δ de ). Hence, already at the linear regime, matter overdensities will create DE underdensities if w < −1. Indeed, in the matter dominated era and for constant w one has [20,23]: In the nonlinear regime, δ de continues to decrease as δ m grows, but the coupling does not vanish as δ de → −1 negative energy densities can appear. To better understand this behavior, bear in mind that voids of matter always have δ m > −1 because the coupling between its density contrast and gravitational force is (1 + δ m ), so as δ m → −1 the coupling vanishes and the decrease in δ m is halted. In the phantom case, (1 + w + δ de ) is always negative inside matter halos, therefore δ de continues to decrease as δ m grows. This pathological behavior might be an indication of a flaw in phantom models with negligible sound speed. This issue is not present for w > −1 because in this case matter halos induce positive δ de and, even in the case of an extreme matter void with δ m −1, δ de is smaller in magnitude than δ m for w −1.
We choose the four sets of parameters shown in Table 1. As initial conditions, we assume that only the growing mode of matter is present and δ de is given by Eq. (2.11).
In Figure 1 we show, for each set of parameters, the ratios of the total growth function, G tot,cs=0 (left panel) and the usual matter only growth function, G m,cs=0 , (right panel) to the w 0 w a S1 -0.9 0.2 S2 -0.9 0.1 S3 -0.9 0 S4 -0.9 -0.1 Table 1. The four sets of parameters w 0 and w a used in this work.
as a function of redshift for each set of parameters given in Table 1. matter growth function for the case of negligible DE perturbations, G m,cs=1 , as a function of redshift. As we can see, the impact of DE perturbations on G tot is much larger. Given this fact, we already expect that our proposed growth function will strongly impact the abundance of galaxy clusters. Since all functions are normalized at z = 0, we can see that clustering DE enhances the growth of perturbations because they grow from smaller values in the past to reach the same value as the homogeneous case today. Also note that the greatest impact of DE perturbations appear for S1. This happens because S1 has the largest value of w at high-z, hence both δ de and Ω de (z) /Ω m (z) are enhanced. This kind of behavior is responsible for large deviations of the growth factor with respect to the so called GR value [31].

Virialization
Now we consider the nonlinear evolution of the fluctuations. First we compute the virialization time, which will be used to define the critical density threshold for collapse, the virial mass and virial overdensity.
If c s = 0, DE fluctuations have no pressure support and follow matter trajectories. Their order of magnitude can be estimated using Eq. (2.11), which also shows that overdensities in matter create overdensities in DE for w > −1. Moreover, as DE fluctuations grow in the nonlinear regime, its local equation, tends to zero for w > −1, thus indicating that its properties become more similar to matter, see [17,32].
Therefore, the total mass, M tot , is expected to have an extra contribution due to DE: The mass associated with matter is and is conserved. It is very debatable how to include DE contribution from a phenomenological point of view. Since only DE fluctuations with negligible sound speed are subjected to the same forces that act on matter particles, we define which is not conserved throughout the evolution. Note that we are not including the contribution from the homogeneous energy densityρ de (that we did include in the case of M m ). This is choice is also helpful to connect to the case of homogeneous dark energy, whose uniform contribution is not usually taken into account when calculating the halo mass.
In order to determine the virialization time, we use the method described in Ref. [33], which takes into account the non-conservation of M tot caused by the inclusion M de . In this approach the virialization time is defined when the moment of inertia of a sphere of non-relativistic particles is in steady state, which yields the equation In EdS, R vir = 1 2 R ta , and δ m (z v ) 145.8, see Ref. [34] for a discussion about this quantity and its impact on the mass function.
In Fig. 2 we show R v /R ta for each set of parameters and for c s = 0 and c s = 1. DE fluctuations affect both the time variation of M tot and R. The general modification is that time derivatives of M tot delay the moment of virialization, so R v /R ta is smaller when DE fluctuations are present.
In Fig. 3 (left panel) we show as a function of virialization redshift. This quantity is important because, as we will show, it is related to the normalization of mass functions and also represents the fraction of DE mass in the halo. It is clear that the contribution of DE increases as DE dominates the background evolution. Since M de gives rise to the non-conservation of total mass, the largest impact of DE fluctuations on quantities computed at virialization is present at very low redshifts. It is important to note that, the larger is 1 + w at high-z, the larger is . As we can see, DE fluctuations can account for a few % of the total mass of the halo. As gives us the fractional contribution of DE with respect to matter, it should be linked to the contribution of the DE perturbation to G tot . It turns out that the function  Table 1.  where G m is the matter growth function in the presence of clustering DE, is very similar to G tot . As one can see from Fig. 3 (right panel), these two quantities differ at most by 2.6%   Table 1. at very low redshifts. This is a important consistency check between linear and nonlinear quantities that we have defined in this work. Expressing the total mass of the halo in the form M tot = 4π 3 R 3 virρ c ∆ v , the virialization overdensity is given by (2.18) In Fig. 4 we show the evolution of ∆ v as a function of redshift for the four different EoS evolution studied. In all cases ∆ v becomes smaller at low-z, as Ω m decreases. The presence of DE fluctuations increases ∆ v , which is a consequence of the delayed virialization that we discussed in the case of R v /R ta .

Critical density threshold
Usually PS and ST mass functions make use of δ c as the critical density threshold that defines collapsed regions. However, when DE perturbations are important, mass is not conserved and the sensible choice is to define the halo properties at the virialization time. Consequently, it seems inappropriate to use δ c for the collapse threshold as in this case the mass function would depend on properties determined at the two different redshifts z c and z v . Therefore, we believe that the most natural quantity to use in the presence of DE fluctuations is the  Table 1. See Section 2.3 for more details. extrapolated linear contrast δ lin tot at the moment of virialization, which we define as δ v in order to avoid confusion with the usual δ c . In Sec. 3 we show that this change can be absorbed into a redefinition of one parameter of the ST mass function.
In Fig. 5 we show the evolution of δ v as a function of z. As for the other quantities, the largest impact of DE fluctuations occur for S1. Note that, for both homogeneous and clustering DE, δ v grows as function of z, which is opposite to the usual behavior of δ c (see, e.g., [21][22][23]35]). The impact of DE fluctuations in δ v can be as large as 7% for S1, whereas the impact in the usual δ c is at most 1% [23]. These results are valid only for DE with negligible c s . For some range of small but non-negligible c s , DE fluctuations are not negligible and can be treated linearly, [22]. In this case, the impact of DE fluctuations is smaller, but induces mass dependence, δ c (z) → δ c (z, M ), which is not considered in this work.

Halo mass function
The halo mass function f (M, z) gives the fraction of the total mass in halos of mass M at redshift z. It is related to the (comoving) number density n(M, z) by: where ρ cc is the comoving average of the energy density that makes up the collapsed halos.
In the usual case of homogeneous DE ρ cc (z) = ρ mc ≡ a 3ρ m , where the latter is the constant matter density in a comoving volume and a is the background scale factor.
The halo mass function acquires an approximate universality when expressed with respect to the variance of the mass fluctuations on a comoving scale r at a given redshift z, ∆(r, z). The latter is given by: where W (kr) is the Fourier transform of the top-hat window function and∆ 2 (k, z) is the dimensionless power spectrum extrapolated using linear theory to the redshift z: In the previous equation P (k, z) is the power spectrum, n s is the spectral index, δ H0 is the amplitude of perturbations on the horizon scale today, G tot (z) is the normalized linear growth function given in Eq. (2.10) and T (k) is the transfer function [36]. Note that the use of G tot (z) accounts for the contribution of DE to the power spectrum. Since we are considering only the cases of c s = 1 and c s = 0, the transfer function shape should remain the same because DE perturbations are either negligible or have the same scale dependence of matter perturbations, thus changing only the overall normalization which is fixed via where σ 8 takes the fiducial value given earlier. Note that we are normalizing all the models to the same present-day value of σ 8 . For the impact of general c s in the power spectrum see [37]. We model the mass function according to the functional form proposed by Sheth and Tormen (ST) [38]: In the previous equation a and p are free parameters while A is fixed by the normalization condition´f dM = 1, according to which: We adopt here the original values from [38], which are a = 0.707 and p = 0.3 and are expected to be approximately universal (see, however, [39]). In the case of clustering DE this is justified by the fact that for c s = 0 the shape of the power spectrum is not changed. However, we should point out that, to our knowledge, a numerical simulation that includes clustering DE has not yet been run and, therefore, the ST parameters a and p may depend on the DE properties. A possible change will likely increase the impact of inhomogeneous DE on the halo abundances. Consequently, by neglecting this effect our results will be more conservative. We adopted the ST mass function as it depends on the peak height parameter ν = δ c /∆ and not on ∆ alone, as in the case of the precise mass function proposed by [1]. This makes the ST function more sensitive to the physics that goes into the halo collapse that we are studying in this paper. In order to obtain a more precise mass function, the usual approach -adopted, for example, in the context of non-Gaussianities [40] and baryonic feedback [41] is to multiply the more precise but less cosmology-dependent mass functions by a correcting factor, which in our case is: We shall then focus on this correcting factor when presenting our results.
As previously discussed, we define halos' properties at the moment of virialization. The ST mass function, however, uses the critical contrast at collapse, δ c . Therefore, we need to adjust the parameter a -which is degenerated with δ c -in order to account for the different critical contrast used:ã where within the EdS mode it is δ v 1.583 and δ c 1.686. The ST mass function becomes: Note that, because of the change of variable, f ST is related to our original definition of f by: In the latter, we have ∆(M, z) = ∆(r(M, z), z), where the function r(M, z) relates the comoving scale that collapsed to form the halo with the actual halo mass. Within homogeneous DE this relation is simply: that is, thanks to mass conservation, one can take the mass corresponding to the collapsing sphere at early times when the perturbation is negligible and ρ m ≈ρ m . Therefore, one get a "background level" relation r = r(M ) and it is not necessary to know the virialization time and the corresponding overdensity and halo radius. As discussed in Section 2.2, when dark energy perturbations are present the total mass M tot is not conserved. Hence, we need to define the actual halo mass and the radius at the virialization and the relation r(M, z) is given by:  Figure 6. Plot of the ratio b v,cs=0 /b v,cs=1 as function of redshift for the four sets of w 0 and w a given in Table 1. See Section 3 for more details.
where the physical halo radius is is independent of the comoving scale r and so of the halo mass M .
If only matter contributes to the mass, M tot = M m (z i ) = M m (z v ) is conserved and we can define the conserved b v according to 14) which is equal to the general case of Eq. (3.13) if dark energy fluctuations are negligible. We note that the common practice of neglecting δ m (z i ) in (3.11) induces an error ∼ 0.1% in b cons v . In Fig. 6 we show the evolution of the ratio b v,cs=0 /b v,cs=1 for each set of parameterizations. In all cases, the modifications with respect to b cons v grow with redshift. The largest differences occur for S3 and S4, a direct consequence of the impact of DE fluctuations on R v as shown in Fig. 2. Fig. 7 shows how different is the non-standard r = r(M, z) relation with respect to the background relation r = r(M ) (left panel) and the impact of the non-standard r = r(M, z) relation on the logarithmic derivative of (3.10). Here we only show the results for the case S1, which suffers the largest modifications. As one can see, the impact of clustering DE in the mass-scale relation is around 1% at z = 0 (left panel of Fig. 7), however the actual impact on the mass function is one order of magnitude smaller because it depends on a logarithmic derivative (right panel of Fig. 7) and can be safely neglected. . These plots refer to the set S1 given in Table 1.
Finally, the last quantity that we need to evaluate is ρ cc . The total average mass density in halos at redshift z can be obtained from the following equation: (3.15) where it has been assumed that, at redshift z, all halos have the same ratio (z) given in Eq. (2.16). This neglects the fact that large halos are produced by the mergers of smaller halos that have virialized earlier and so with a different virialization overdensity and DE contribution. As we can see, the impact on the mass function is linear in (left panel of Fig. 3), so it increases abundances as much as 5% on all mass scales.

Results
We show in Fig. 8 the quantity n(M,z)| cs=0 n(M,z)| cs=1 for the four sets of w 0 and w a given in Table 1 and at three different redshifts. As one can see, DE perturbations can impact halo abundances at the level of 10%-30%, depending on the halo mass and redshift. Bear in mind that all models have the same σ 8 . In this case, clustering DE always diminishes the abundance of massive clusters M > 10 14 M . Interestingly, the number of small halos, 10 10 M < M < 10 12 M , is enhanced at z = 0 and z = 1. The main impact of clustering DE for large masses comes from the modifications in the quantity δ v /G tot , which is always larger than in the homogeneous case and strongly reduces the number of halos around the exponential tail of the mass function. For small masses, the main impact is related to the density normalization of mass function, Eq. (3.15).
Had we normalized the amplitude of fluctuations at the redshift of CMB decoupling, clustering DE would have generated more massive clusters in comparison to homogeneous DE, therefore worsening the Planck Clusters tension. Clustering DE models alleviate this tension only if w < −1, in which case the impact of DE fluctuations is opposite to the case w > −1. However, as we discussed in Sect. 2.1, phantom clustering DE can generate negative energy densities, δ de < −1. We will address this problem in a forthcoming paper [6]. Perturbations are normalized today so that all models have the same σ 8 . The four sets S1-4 of w 0 and w a are given in Table 1.
There are numerous ongoing and planned surveys (e.g., DES [42,43], J-PAS [44,45], Euclid [46][47][48]) which will count the number of high-mass systems on the sky in an effort to better constrain the cosmological parameters that control the growth rate of structure. As we have seen the impact of DE perturbations on the mass function is greatest at large masses. Therefore, it is interesting to quantify how much cluster counts are affected by such perturbations.
Here, we neglect the effect of the mass-observable relation and we do not bin in mass. We will carry out survey-specific analyses in a forthcoming work. We define a cluster as an object with mass M ≥ M thr ≡ 10 14 h −1 M , a figure in agreement with the expected sensitivity from Euclid. The total number of clusters at the redshift z and within the redshift bin ∆z is then: where the quantity dV /dz is the cosmology-dependent comoving volume element per unit redshift interval which is given by: where d A is the angular diameter distance and H(z) is the Hubble rate at redshift z. . Fractional change in cluster number counts due to dark energy perturbations. N (z) is computed using (4.1) and N 0 (z) refers to dark energy with sound speed c s = 0, and similarly for N 1 (z). Perturbations are normalized today so that all models have the same σ 8 . The four sets S1-4 of w 0 and w a are given in Table 1. Figure 9 shows that cluster counts are modified by about 30% at a redshift of unity. Therefore, cluster counts is a very promising observable if one is to constrain models of dark energy that feature negligible sound speeds.

Conclusions
In this paper we have discussed in detail how the halo mass function is modified when dark energy has a negligible sound speed that causes its perturbations to collapse. We have developed a consistent framework according to which the halo properties are defined at the moment of virialization and that takes into account the non-conservation of the mass relative to the dark-energy fluctuations.
To summarize, the corrections that directly impact the mass function due to the presence of DE fluctuations are : • Growth function: based on the spherical collapse model in the radius approach, we have argued that the total perturbation density should be used, Eqs. (2.8) and (2.10). This new function differs from matter growth function in the homogeneous case by 3% to 8% for the equations of state we have analyzed.
• Threshold density for collapse: by demanding consistency in the definition of halo mass and the mass function quantities, we have argued that most natural threshold density to be used is the total linear perturbation density at the moment of virialization, δ v , Eq. (2.19). The impact of clustering DE on this quantity is similar to the one on the growth function, but has a distinct redshift evolution in comparison with the usual δ c .
• Mass-scale relation: since the total mass of a halo with some fraction of DE fluctuation is not conserved, one can not use the background mass-scale relation for matter only.
The corrected relation is given by Eqs. (3.12) and (3.13). This, however, is a subpercent correction and can be safely neglected.
• Normalization density: also because DE fluctuations induce non-conservation of total mass, the mean energy density that normalizes the mass function has to be corrected, Eq. (3.15). This is a linear modification of the order of the DE mass fraction and so it changes the mass function by a few % on all mass scales.
We have then computed the impact of clustering dark energy considering all the mentioned modifications on the mass function and found that halo abundances can be altered at the level of 10%-30%, depending on the halo mass and redshift. As the change is largest at the high masses of galaxy clusters, we have then computed the impact on cluster number counts and obtained that they are modified by about 30% at a redshift of unity, as big an effect as the one that comes from including baryons within cosmological simulations [41]. A comprehensive analysis that takes into account the nuisances from the scaling relations will be the subject of forthcoming work [6].