A short review on clustering dark energy

We review dark energy models which can present non-negligible fluctuations on scales smaller than Hubble radius. Both linear and nonlinear evolutions of dark energy fluctuations are discussed. The linear evolution has a well-established framework, based on linear perturbation theory in General Relativity, and is well studied and implemented in numerical codes. We highlight the main results from linear theory to explain how dark energy perturbations become important on the scales of interest for structure formation. Next, we review some attempts to understand the impact of clustering dark energy models in the nonlinear regime, usually based on generalizations of the Spherical Collapse Model. We critically discuss the proposed generalizations of the Spherical Collapse Model that can treat clustering dark energy models and their shortcomings. Proposed implementations of clustering dark energy models in halo mass functions are reviewed. We also discuss some recent numerical simulations capable of treating dark energy fluctuations. Finally, we make an overview of the observational predictions based on these models.


Introduction
Since the discovery of the accelerated expansion of the Universe Riess et al. [1], Perlmutter et al. [2], a great variety of explanations have been proposed. The most simple and well-studied proposal is that the accelerated expansion is caused by the Cosmological Constant, Λ, which is constant in space and time and naturally possesses no fluctuations. In this model, Dark Energy (DE) only modifies the background cosmological evolution, then the modifications in linear and nonlinear evolution of cosmological perturbations are straightforward to implement. Together with the Cold Dark Matter (CDM), which accounts for roughly 25% of the Universe energy density, the ΛCDM model provides an outstanding description of cosmological data obtained so far, e.g., Aghanim et al. [3], Abbott et al. [4], Asgari et al. [5] Although ΛCDM is very successful in describing almost all cosmological observations, on theoretical grounds, it is challenged by the Cosmological Constant Problem Weinberg [6], Carroll [7] and the Cosmic Coincidence Problem Zlatev et al. [8]. More recently, direct measurements of the Hubble constant, H 0 , have also been challenging the ΛCDM model, showing a disagreement of about 5σ with respect to the inferred value from Cosmic Microwave Background (CMB) data Verde et al. [9], Di Valentino et al. [10]. A less significant tension, about 3σ, in the ΛCDM model is related to the predicted normalization of the matter power spectrum, σ 8 , the S 8 = σ 8 √ Ω m0 /0.3 parameter Perivolaropoulos and Skara [11], Di Valentino et al. [12]. A recent analysis of possible solutions can be found in Schöneberg et al. [13].
Given these difficulties with Λ, a profusion of alternative models to explain the cosmic acceleration were proposed. One of the first and most popular alternatives to Λ is the quintessence class of models. In these models, a new scalar field minimally coupled to gravity and with no direct interactions to other types of matter plays the role of DE. This kind of model was studied even before the discovery of accelerated expansion by Peebles and Ratra [14], Wetterich [15]. In quintessence models, DE is time-dependent and its Equation of State arXiv:2204.12341v1 [astro-ph.CO] 26 Apr 2022 (EoS) parameter track the background, possibly alleviating the Coincidence Problem Zlatev et al. [8], Caldwell et al. [16], Steinhardt et al. [17].
Since the quintessence field is dynamical, it necessarily has fluctuations. These fluctuations, however, are relevant only on scales of the order of Hubble radius Ma et al. [18], Brax et al. [19], DeDeo et al. [20]. On the small scales of interest for structure formation, quintessence perturbations are much smaller than matter (dark matter plus baryons) perturbations and thus are usually neglected. This tiny amount of DE perturbations on small scales is a consequence of the sound speed of the field perturbations, which has a constant value c s = 1. This value of the sound speed implies that the sound-horizon scale of quintessence, a scale below which the perturbations are strongly suppressed by pressure support, is of the order of Hubble radius, c s /H 0 .
One of the first proposed models beyond the quintessence that can possible present relevant perturbations on small scales is the tachyon scalar field Sen [21], Padmanabhan [22], Bagla et al. [23]. The main difference with respect to quintessence is that the tachyon field has a time-varying speed of sound, which can be smaller than unity, thus allowing the field to cluster more effectively on smaller scales.
It was also observed that even more general scalar field models could be constructed. The so-called k-essence models were initially proposed in the context of Inflation Armendariz-Picon et al. [24], Garriga and Mukhanov [25]. In such models, one has the freedom to choose both the kinetic term and the potential of the scalar field, which translates into liberty to choose the EoS parameter and c s . Therefore, in this class of models, DE can have an arbitrarily low speed of sound, thus its perturbations can be the same order of magnitude of matter perturbations. In this scenario, DE has the potential to impact structure formation beyond the background level.
More recently, the class of Horndesky theories, Horndeski [26], was rediscovered and it was shown that both quintessence and k-essence models are sub-classes of Horndesky. At the linear perturbation level, these sub-classes are parameterized by the α K parameter Bellini and Sawicki [27]. In the context of Horndesky theories, many models beyond k-essence exist, including non-minimally coupled scalar fields, which also modify gravitational interaction. The various types of models that can explain the accelerated expansion can also be described in the Effective Field Theory framework Gubitosi et al. [28]. Many of these proposals are discussed in Amendola et al. [29].
In this review, we will focus on DE models described by minimally coupled to gravity scalar fields. We also restrict our attention to fields with no direct interaction with other matter fields. The k-essence class of models can be parametrized as perfect fluids defined by two time-dependent functions, w(t) and c s (t) and this description will suffice to analyze how large DE fluctuations can be and estimate their observational impact.
The first step to understand DE fluctuations is to study them at the linear perturbative level, which is described by the well-established theory of linear cosmological perturbations, e.g., Kodama and Sasaki [30], Mukhanov et al. [31], Ma and Bertschinger [32]. In this context, the study of DE perturbations is straightforward, but limited to large scales that did not develop nonlinear matter fluctuations. Difficulties arise when trying to study DE fluctuations in the nonlinear regime. Historically, structure formation was studied using the Newtonian theory, which can not deal with relativistic fluids, such as DE. The obvious approach of using full-blown General Relativity to include DE fluctuations in structure formation studies is undoubtedly challenging. In fact, relativistic studies of structure formation with Clustering Dark Energy (CDE) were developed quite recently Dakin et al. [33], Hassani et al. [34,35].
The first attempt to study Clustering DE (CDE) models in the nonlinear regime used an integration between the Newtonian Spherical Collapse Model (SCM) with a "local" Klein-Gordon equation, modified to permit the clustering of quintessence Mota and van de Bruck [36]. This study was able to show important new effects due to DE fluctuations, which were later confirmed by more general and well-justified models. For instance, depending on the evolution of its EoS, DE fluctuations can become nonlinear, impact the nonlinear evolution of matter fluctuations and change the virialization state of matter halos. Hence, it became clear that CDE can impact structure formation. Later on, some authors constructed a more formal and general framework to include CDE in the SCM, Abramo et al. [37], Creminelli et al. [38], Basse et al. [39].
The main focus of this review is to describe and discuss the applicability of the generalizations of the SCM capable of treating CDE models. We also pay special attention to the corresponding modifications on the Halo Mass Functions (HMF), which, up to now, have not been explored by numerical simulations. Moreover, we review the impact of CDE on cosmological observables and prospects for its detection.
The plan for this review is the following. Section 2 presents the essentials of relativistic perturbation theory that describe DE perturbations and their scale dependence. In Section 3, we review quintessence and k-essence models, highlighting the conditions under which relevant DE perturbations can be present. Section 4 discusses the SCM and its generalizations to study homogeneous and inhomogeneous DE models. Section 5 presents the Pseudo-Newtonian description of the SCM, which permits direct contact with linear relativistic perturbations on small scales and provides a clear and general framework to study CDE models in the nonlinear regime. We discuss the impact of DE fluctuations on the critical density threshold and HMF in Sections 6 and 7, respectively. The observational impact of CDE is reviewed in Section 8. We discuss the main results and perspectives in Section 9.

Linear perturbations
To understand the basic behaviour of DE perturbations, it is enough to consider scalar perturbations in the absence of anisotropic stresses. In this case, the perturbed line element in the Newtonian gauge is given by Ma and Bertschinger [32] where η is the conformal time and Φ the gravitational field. The energy momentum of a perfect fluid with energy density ρ, pressure p and four-velocity u µ is given by which include background (overbar quantities) plus perturbed quantities: We also restrict the analysis to the matter dominated era and the actual DE dominated phase. Then, at background level, we have the Friedman equations including matter (barions plus dark matter, indicated by the subscript m) and dark energy (indicated by the subscript de) given by where w =p de /ρ de is the EoS parameter for DE and the prime indicates derivative with respect to the conformal time, η. Density parameters of matter and DE are defined by and where ρ c =ρ m +ρ de is the critial density.
In this review, we present examples using the CPL parametrization of the EoS Chevallier and Polarski [40], Linder [41] where w 0 and w a are constants. Since Λ gives a good description for the background evolution, we fix w 0 = −1 and will vary w a to show the impact of DE fluctuations in a scenario that is very similar to ΛCDM at low redshift. In Fourier space, the (00) component of Einstein equations is given by where δ = δρ /ρ is the density contrast for a given component identified by the subscript . The conservation equations for perturbations of each fluid component, ∇ ν δT µ ν = 0, can be written as where θ = ik j v j is the divergence of the fluid peculiar velocity and c 2 a =p /ρ is the adiabatic sound speed of the fluid. The pressure perturbation is given by Bean and Dore [42] where c 2 s = (δp/δρ) rest is the sound speed of the fluid in its rest frame. An adiabatic or barotropic fluid is defined by c s = c a . The fluid is said to possess entropy perturbation if c s = c a . In general, DE models have entropy perturbations, but it is also possible to construct adiabatic models, e.g., Linder and Scherrer [43], Unnikrishnan and Sriramkumar [44] Let us analyze the behaviour of perturbations in an adiabatic DE model, which has simpler equations -see Ballesteros and Lesgourgues [45] for a more general treatment. In this case, the last term in (10) vanishes. Assuming w and c s as constants, from (8)-(10) we can determine a second order equation for DE density contrast Let us focus on small scales, k 2 H 2 , H and k 2 Φ HΦ , Φ . During the matter dominated era, we have the well-known solution δ m ∝ a and Φ = constant, and equation (11) simplifies to For non-negligible c 2 s , equation (11) has a constant solution Since, from (7), δ m ∼ k 2 Φ, DE perturbations are usually negligible on small scales when compared to the matter perturbations.
For negligible c s , equation (11) has the following solution Abramo et al. [46], Sapone et al. [47], Creminelli et al. [48] which is a good approximation even for Early DE models Batista and Pace [49]. As can be seen, the magnitude of DE perturbations is determined by both w and c s . If w −1, as should be at low-z, DE perturbations are strongly suppressed regardless of c s . For DE models in which w deviates from −1 at higher redshift and with low or negligible c s , DE perturbations can be as large as matter perturbations in the matter-dominated era, as indicated by (14). As we will see, in this case, DE fluctuations can become nonlinear and impact structure formation.
The solution (14) also indicates a correlation between matter and DE perturbations. Positive matter perturbations, which will be associated with the formation of halos in the nonlinear regime, will induce DE overdensities in the case of 1 + w > 0 and underdensities if 1 + w < 0. As we will see later, in the nonlinear regime, this correlation can induce a pathological behaviour for phantom DE models (w < −1), namely, δ de < −1.
The comoving scale bellow which DE perturbations are strongly suppressed by its pressure support is given by the sound horizon here H =ȧ/a is the usual Hubble parameter and the dot indicates derivative with respect to time. In regions below this scale, the pressure support halts the growth of DE perturbations. As we will see, quintessence models have c s = 1, then r s ∼ 1/H 0 and their perturbations are very small on scales below the Hubble radius. Let us estimate the value of c s that induce large DE perturbations on small scales. For instance, c s = 10 −3 gives a sound horizon r s 14Mpc (assuming best fit Planck18 cosmology Aghanim et al. [3]). This value is associated with a mass scale M = 4π/3ρ m0 (r s /2) 3 2 × 10 14 M , where ρ m0 is the matter density now. The value c s = 10 −3 was indeed reported to impact the formation of halos of such mass by Basse et al. [39]. As we will discuss later, in this work the authors considered that DE fluctuations do not break linear approximation. In the limit of c s → 0, however, DE fluctuations can become nonlinear, depending on the value of w. Nevertheless, we can consider that impact of DE fluctuations to become important for nonlinear structure formation for c s < 10 −3 .

Dark Energy Models
Now that we have understood under which circumstances DE perturbations can be significant on small scales, let us discuss some DE models that can present such large perturbations. The variety of DE models is enormous, with several reviews, books and extensive analysis about them, e.g., Copeland et al. [50], Amendola and Tsujikawa [51], Yoo and Watanabe [52], Tsujikawa [53], Ade et al. [54]. We will focus on the k-essence class of models, which provide enough generality for w and c s to permit large DE perturbations.
Models of k-essence were initially introduced in the context of inflation Armendariz-Picon et al. [24], Garriga and Mukhanov [25] and shortly after applied to describe DE Armendariz-Picon et al. [55], Erickson et al. [56]. A cosmological model with k-essence is described by the following action where κ 2 = 8πG, R is the Ricci scalar, X = − 1 2 g µν ∂ µ ϕ∂ µ ϕ, S m the action for matter fields (photons, neutrinos, baryons, dark matter) and L(X, ϕ) the Lagrangian for the k-essence field, ϕ. The k-essence energy-momentum tensor is given by where subscript , X represents the derivative with respect to X. Comparing it to the perfect fluid energy momentum tensor (2), we find that The sound speed of k-essence perturbations is given by Garriga and Mukhanov [25] In this framework it is clear that one has the freedom to choose the functions w and c s . Next, let us have a look at some specific popular realizations of DE models.

Quintessence
In the k-essence language, quintessence is defined by L = X − V(ϕ), then its EoS parameter, defined by p/ρ given by (18), and sound speed, given by (19), read If the kinetic energy of the field dominates over its potential energy, we have w 1. In the opposite case, w −1. Its sound speed, however, is always the same. Hence, although w can be far from −1, the sound speed of quintessence does not allow for large perturbations on small scales. This fact can also be understood via the Klein-Gordon equation for quintessence perturbations, which can be written as chan Hwang and Noh [57] δφ + 3Hδφ + c 2 where we explicitly introduced the the parameter c s . Quintessence models must have V ∼ H 2 0 in order to accelerate the Universe expansion recently, where H 0 = H(t 0 ) is the Hubble constant. Thus, for scales well bellow the Hubble radius, c 2 s k 2 dϕ 2 , and the field perturbations are strongly suppressed on these scales.
Since quintessence was initially the most popular alternative model to Λ, initially, many cosmologists considered DE perturbations were effectively negligible on small scales. This picture started to change with the appearance of new DE models which can have c s < 1.

Tachyon
A well-known model for DE that makes use of non-canonical scalar field is the tachyon field. It was introduced in the context of string theory Sen [21,58], but can also be understood as a generalization of the relativistic particle lagrangian, Padmanabhan and Choudhury [59]. The tachyon model is defined by the following lagrangian Comological studies of this model include Padmanabhan [22], Bagla et al. [23], Padmanabhan and Choudhury [59], Abramo and Finelli [60]. The EoS parameter and its sound speed are given by w = 2X − 1 and c 2 s = −w (23) The accelerated expansion occurs when X 1, thus w −1 and c 2 s 1. In this case, the behaviour is very similar to quintessence. However, earlier in the cosmological history, the kinetic term can be larger, which yields a smaller c s . In this case, both w and c s provide the conditions for large DE perturbations.
Considering a power-law potential, V(ϕ) ∝ ϕ −α , the actual impact of tachyon perturbations on CMB is small because the accelerated expansion requires α 0, which, in turn, suppresses the field perturbations Abramo et al. [61]. Nevertheless, this model clearly shows that DE perturbations are not necessarily small on scales below the Hubble radius.
For a constant potential, the tachyon model is equivalent to the Chaplygin gas Kamenshchik et al. [62], with equation of state where A is a constant. This model and its generalized version was studied by many authors, e.g., Bento

Clustering DE
It is possible to construct DE models with negligible sound speed. In the context of k-essence, a model with a constant arbitrary c s is given by Kunz et al. [72] where M is a constant mass scale. In the context of quartessence, Scherrer [73] proposed a model with very low c s . A possible issue when building such models is that distinct lagrangians can have the same w and c s , Unnikrishnan [74].
In effective field theory, models with c 2 s ∼ −10 −30 were described by Creminelli et al. [48]. Although negative c 2 s yields gradient instabilities, this value is so small that no relevant effect on cosmological scales is expected. This work also concludes that no pathological behaviour is present for phantom models, w < −1. However, as we will show later on, in the nonlinear regime, it is possible that phantom DE with negligible sound speed can present δ de < −1, i.e., a pathological situation with negative energy density, ρ de =ρ de (1 + δ de ).
Models with identically zero sound speed using two scalar fields were proposed by Lim et al. [75]. In Horndeski theories, there are even more possible realizations of DE with low sound speed, which can be chosen as a function of four physical parameters Bellini and Sawicki [27].
Given the variety of possible CDE models, the phenomenological implementation of parametrizations of w and c s is valuable to explore the possible impact of DE perturbations on structure formation. As already mentioned, we will show examples using w = w 0 + w a (1 − a) and constant c s .

The spherical collapse model
Once we know the basic behaviour of DE linear perturbations, we may ask how they impact the structure formation on small scales. The SCM is an analytical approach first proposed to study the nonlinear evolution of matter perturbations in the Universe. It can be used in connection to analytic and semi-analytic halo mass functions to determine the abundance of matter halos, which we will describe in section 7. In this section, we will review the classic formulation of the SCM and the main quantities of cosmological interest. We also discuss some generalizations capable of treating homogeneous DE and the first proposed model that considers the nonlinear evolution of DE fluctuations. Later on, we will present a more general framework for the SCM and a detailed discussion about threshold density, which is an important quantity that determines the abundance of halos.

Einstein-de-Sitter Universe
The SCM was first proposed by Gunn and Gott [76]. The important conclusions for our porpouses are described in Padmanabhan [77] and Sahni and Coles [78]. This model describes the dynamics of spherical shell of radius R that encloses a mass M m . The dynamical equation for the shell radius isR where the dot indicate derivative with respect to time. The mass is given by where δ m is assumed to be the averaged density contrast inside the sphere of radius R, which can also be understood as a top-hat profile for the matter fluctuation. The total mass is supposed to be conserved, dM m /dt = 0, which yields the following equation for the matter contrastδ which has the solution where the subscripts i indicate the initial values. The first integral of (26) is given by where E is a constant of integration identified as the total energy of the shell -the kinetic energy (K =Ṙ 2 /2) plus the potential energy per unit of mass (U = −GM m /R). Note that, in the derivation of (30), the fact that M m is constant was used. As we will see later on, in the presence of DE, the effective mass inside the shell is not conserved anymore.
Assuming that the shell initially expands with the background, R ∝ a, we haveṘ i = H i R i . Using (27), the total energy can be expressed as where Ω mi = Ω m (a i ). For 1 + δ mi > 1/Ω mi , the total energy of the shell is negative. Thus it is expected that the initially expanding dust cloud will grow to a maximum radius, contract and eventually form a bounded structure. In EdS Universe, Ω m = 1, and any region with positive δ mi will eventually collapse. As we will see next, this is indeed the behaviour given by solving for the time evolution of R, which will describe the formation of a matter halo. On the other hand, if 1 + δ mi < 1/Ω mi , the total energy of the shell is positive and the cloud will expand forever, forming a void with −1 < δ m < 0. The maximum shell radius is given byṘ = 0 and is called the turn-around radius, given by This is indeed a maximum radius because equation (26) indicates thatR < 0. This fact also implies that no solution of a minimal radius exists in the usual SCM, and consequently, the model can not describe an equilibrium configuration. This shortcoming is usually circumvented by assuming that the system achieves virial equilibrium at some point of its evolution and that this state represents the final bounded structure. There were efforts in implementing virialization in the dynamics of the SCM, see Engineer et al. [79], Shaw and Mota [80]. Equation (26) can be solved analytically with the following parametrization where A = R ta /2, B = t ta /π. Substituting these expressions in (26) one gets the relation As can be seen from (33), the shell radius is maximum at θ = π (turn-around) and then begins to shrink. Formally, as x → 2π, R → 0 and, from (29), δ m → ∞. This is identified as the moment of collapse. In reality, these values are not achieved because, for a given shell approaching collapse, other inner shells containing collisionless matter have already crossed each other. Thus, the total mass within a specific radius is not conserved anymore, signalling the break down of the model. Nevertheless, the collapse time is used to define quantities used to estimate the abundance of halos.
Using conservation of mass, given by equation (27), and the background density evolution in EdSρ m = 6πGt 2 −1 , the nonlinear evolution of matter fluctuations is given by where the superscript NL indicates the nonlinear value of the density constrast. Expanding for small x, we get the linear evolution where the superscript L indicates the linear value of the density contrast. At turn-around, we have and As can be seen, at this time, the linear evolution gives a density contrast of order one, indicating the transition from linear to nonlinear regime.
The collapse threshold (or critical density contrast) is the extrapolated linear overdensity value above which a bound structure is considered formed, a halo. Usually, it is defined at the collapse time As we will see later, the value δ c is of great importance in structure formation studies that make use of analytic or semi-analytic mass functions. This quantity is modified in the presence of homogeneous and inhomogeneous DE. It is also important to describe the state of equilibrium of the halo, which is assumed to obey the virial theorem for non-relativistic particles. Virialization occurs when U( From this, we can compute the virialization time (x = 3π/2) The collapse time ( With these quantities, two definitions of virialization overdensity are usually found in the literature. One possibility is define the virialization overdensity with respect to the matter background density Another common choice uses the critical density as reference Of course, ∆ cm and ∆ cc give essentially the same results in a matter dominated Universe, but they differ when DE becomes important, see figure 1. The overdensity ∆ ≈ 178 became a reference value used in N-body simulations in order to identify halos, but many other definitions can be used, see Despali et al. [81] for an analysis of several such choices. Interestingly, this paper has shown that the use of the quantity ∆ cc yields more universal mass functions, in the sense of its dependence on redshift and cosmological parameters. We will return to this discussion when dealing with the impact of CDE on halo mass functions.
Yet another possibility is to compute both the total and background densities at the virialization time, as suggested by Lee and Ng [82]. In this case the virialization overdensity is given by The value of ∆ v is slightly smaller than ∆ cc . Later on, will show that in the presence of DE these two quantities decay with redshift, whereas ∆ cm grows (see figure 1). Analogously, we can define the threshold density at virialization

Spherical collapse model with homogeneous dark energy
In the presence of homogeneous DE, the SCM has to be modified to include the effects of this new component on the dynamics of the shell radius. Now we havë where is the effective mass associated with homogeneous DE. This definition can be understood via the Poisson equation in the presence of a relativistic fluid The first study of the equation (46) was done by Lahav et al. [83] considering Λ as a DE component. It was found that the ratio of virial to turn-around radius is smaller in the presence of Λ, given by the following expression where η = Λ/4πGρ m (t ta ). Thus, in the presence of Λ, the halo has to contract more with respect to EdS to achieve the virial equilibrium. This indicates that halos formed in the late Universe have a distinct structure than those formed at high-z, when DE was very subdominant.
Other authors have studied the SCM in the presence of Λ, e.g., Lacey and Cole [84], Kitayama and Suto [85]. Although analytical solutions for δ c and ∆ cc were found, the following fits presented by Kitayama and Suto [85] became popular: and where w f (z) = 1/Ω m (z) − 1. Expression (50) shows that the the influence of Λ on the critical threshold is very small, giving δ c (z = 0) 1.676 for Ω m0 = 0.3 (only 0.64% different than the EdS value). However, the change in virialization overdensity is much larger. We have ∆ cm (z = 0) 334.2 for Ω m0 = 0.3, about two times EdS value. Studies of DE with constant w were done by Wang and Steinhardt [86], Weinberg and Kamionkowski [87]. The fitting functions proposed in the latter, shows variations of δ c (z) bellow 0.56% for −1 < w < −0.8 with respect to the ΛCDM fit. For ∆ cm the differences are lower than 17%. The differences can be slightly larger in some quintessence models Mainini et al. [88] (42), (43) and (44) for the ΛCDM model with Ω m0 = 0.3. In this plot, ∆ v is determined using the method summarized by equation (78), while ∆ cm and ∆ cc by the proper numerical method discussed in section (6.1) In figure 1, we show a plot of these three different definitions of virial overdensity. As discussed, the variations in δ c due to homogeneous DE are very small and will be shown together with the case of CDE in figure 2.

Spherical collapse model with inhomogeneous dark energy
The first study of SCM in the presence of CDE was done by Mota and van de Bruck [36]. DE was modeled as quintessence field, which equation of motion in background is given bÿ Remind that, in the context of quintessence models, the fluctuations of such field must be tiny on small scales because c s = 1. The authors have considered, however, that the field can cluster on small scales, assuming that, inside the collapsing region, the it obeys the following equationφ where is a quantity that describes the flux of DE in the collapsing region. If α = 1 DE evolves in the same way inside and outside the collapsing region, then DE is homogeneous, and ϕ c = ϕ. On the other hand, if α = 0 the field can evolve differently than in the background, allowing for DE fluctuations. The impact on δ c was computed in Nunes and Mota [89], where differences of a few per cent with respect to ΛCDM and redshift dependent features associated with the evolution of w were found. Although this model was a breakthrough regarding the possibility of nonlinear fluctuations of DE, the implementation of equation (53) is ad doc and formally inconsistent with the Klein-Gordon equation, (21). Equation (53) can be interpreted as (21) with c s = 0, where the the gravitational coupling is encoded in the dynamics of R. Later it was shown that the canonical scalar field indeed has negligible perturbations on small scales Mota et al. [90], Wang and Fan [91]. However, its perturbations can be more important on voids, which are much larger than halos.
Despite this inconsistency, this model presented some of the important impacts of DE fluctuations on the nonlinear regime, which were later confirmed by other studies. In particular, the following findings can be highlighted: 1.
The amount DE fluctuations strongly depends on the evolution of w.

2.
DE fluctuations impact the nonlinear evolution of δ m and virialization of halos. 3.
The  [108]. In the next section we will present a framework that can describe the nonlinear evolution in the CDE model.

Spherical collapse model in the Pseudo-Newtonian Cosmology
General Relativity is the standard theory to study gravitational phenomena of fluids with relativistic pressure, p ∼ c 2 ρ. In the context of structure formation, the gravitational fields are small, and the weak field limit can be used. Therefore, one can consider using the Newtonian theory, but taking into account the effects of relativistic pressure. This approach was indeed used in McCrea [109], Harrison [110] to study the background evolution of the Universe. It was shown that Newtonian equations with the correction for relativistic pressure reproduce the Friedman equations.
However, the application of these modified Newtonian equations for the study of cosmological perturbations is more controversial. Sachs and Wolfe [111] showed that perturbations in fluids with p = 0 do not agree with the relativistic treatment. For about three decades, the use of Newtonian equation for cosmology was halted until the origin of the discrepancy reported in 1967 by Sachs and Wolfe was found and treated in Lima et al. [112].
According to Lima et al. [112], the system of equations (the subscript identifies the different fluids under consideration), which we call Pseudo-Newtonian Cosmology (PNC), provides the correct growing modes for linear perturbations of fluids with relativistic pressure.
The PNC can also be used to generalize the SCM Abramo et al. [37,46]. The perturbed equations in comoving coordinates with background are given bẏ where v is the peculiar velocity of the fluid. Clearly, these equations are more general then those of the usual SCM. Let us consider the simplifying assumptions that will yield the correspondence between them. Remind that the SCM assumes a top-hat profile, thus we need ∇δ = 0 inside the shell. The other quantities must be consistent with this assumption. Equation (58) must depend only on time, thus ∇ · v = θ(t). Taking the divergence of (59) we geṫ Now let us turn our attention to the parameter c s . At first glance, assuming that ∇δ = 0 and ∇ · v = θ(t) seems compatible to any choice c s = c s (t). However, considering two distinct fluids, say pressureless matter (w m = c s m = 0) and DE (w de = w de (t) and c s de = c s de (t)), implies that each fluid has its own dynamical equation, (61). This indicates that there is no unique spherical shell radius, because the two fluids can flow with distinct velocities.
Another problem about using a generic function c s (t) arises when extending the analysis to regions outside the shell. For a top-hat profile, δ is discontinuous at the edge of the shell, then ∇δ is ill-defined at this point. A more realistic realization is to assume a smooth decay of δ around the edge. In this case, its gradient is well-defined and there exists a non-null pressure gradient, c 2 s ∇δ, around the shell radius. This, again, would make the two fluids flow with distinct velocities and, more drastically, disrupt the original homogeneous top-hat-like profile.
Therefore, regarding the value of c s , the equivalence of equations (58)-(60) with the usual SCM is achieved only for c s = 0. In this case, all fluids share the same dynamical equation and the evolution of the system is described bẏ For pressureless matter, equation (58) yieldṡ Comparing it with equation (28), we identify that the divergence of the peculiar velocity is given by Inserting this relation in (63), one obtains For the EdS model, we recover the usual spherical collapse equation, (26), (here as a function of the density contrast)R In the presence of CDE, we geẗ Particularizing to homogeneous DE models, δ de = 0, we recover the SC equation in Wang and Steinhardt [86].
The system of equations (62) and (63) was also derived in Creminelli et al. [38] using Fermi coordinates in a relativistic framework. The consistency of equation (63) with k-essence equation of motion was also verified in this work.
Summarizing, the SCM can be generalized to include other fluids with zero pressure perturbation, but allowing for non-zero background pressure. This is just the case of CDE with EoS parameter w. The equations governing the nonlinear evolution of the fluids are then given by:δ Note that, generic values of c s can be considered in the system of equations (58)-(60). However, the differential equations become partial in this case, and the correspondence with the orginal SCM is lost. We stress that the system (69)-(71) is valid only in the limit c s → 0.
A slightly different approach to treat DE perturbations in the SCM was proposed in Basse et al. [39]. In this work, DE perturbations were described in the linear regime, but allowing for non-null sound speed. Then, as discussed previously, DE perturbations can not maintain the top-hat profile. The new idea in Basse et al. [39] was to consider the spherical region with a DE perturbation given by where W(kR) is the top-hat window function in Fourier space, given by (83), and δ L de (k, t) obeys the linearized equations (58) and (59) in Fourier space. This approach has the advantage to allow the use of arbitrary sound speed up to values that do not break the linear approximation for DE perturbations. However, it can not be used in the full clustering regime, c s → 0, when DE fluctuations can become nonlinear. Nevertheless, the approach of Basse et al. [39] is essential to understand the impact of c s in the nonlinear evolution of matter fluctuations, showing that δ c and ∆ v become mass-dependent, which can be an observational signature of DE fluctuations in the abundance of galaxy clusters.

Density threshold definitions
Now that we have described the governing equations for the nonlinear evolution of matter and DE fluctuations, let us turn our attention to the determination of the threshold density that will be used in the Halo Mass Functions (HMF) in section (7). Although this determination is analytic in the usual in the ΛCDM model, in the presence of DE with general w and its possible fluctuations, a numerical computation is necessary, which may introduce some issues.
We will first discuss the calculation of the usual collapse threshold, δ c , and then the alternative virialization threshold, δ v . While the former is historically the most used one, the latter seems to be more consistent with the actual contribution of DE fluctuations in the collapsing region and in better accordance with results from simulations.

Collapse threshold, δ c
In order to determine δ c for CDE models, one has to solve numerically the system of equations (69)-(71) and its linearized version. The initial conditions for matter assume the EdS linear solution, δ m ∝ a and the corresponding value for θ; for CDE the solution (14) can be used.
Having solved the system, one has to determine a criterion to find the moment of collapse. In the standard SCM, the collapse threshold can be found analytically computing the value of the linearly evolved density contrast, (35) at the time of collapse, R → 0 or δ NL m → ∞. This suggests that, in general, the determination of δ c can be done by defining a numerical threshold value for δ NL m , above which the halo is considered to be formed. However, this implementation can introduce a small error in the determination of δ c . From equations (34) and (35), we can see that, for small x, the linear and nonlinear values of matter contrast differ by O x 4 . Thus, their initial values are slightly different. In the analytical solution for δ c , Eq. (38), this difference is naturally taken into account. If one neglects this difference in the initial conditions, the resulting δ c (z) presents a small spurious increase with redshift, moving away from the EdS value at high-z. This problem was noted in Herrera et al. [118] and further discussed in Pace et al. [119].
The approach presented in Pace et al. [119] is twofold: initiate the numerical integration at very high-z (z ∼ 10 5 ) and make a change of variable, δ m → 1/ f . The first measure diminishes the difference between δ NL m and δ L m at the beginning and, consequently, the spurious increase of δ c . The change of variable minimizes the numerical error in the determination of the moment of collapse.
Although this implementation give good results, the proper approach to accurately determine δ c numerically is to use the analytical solutions (34) and (35) to determine the linear and nonlinear initial conditions Batista and Marra [120]. These analytical solutions assume that the peculiar velocity is zero initially, but more general expressions can be found in Padmanabhan [77]. This method, however, might not be appliable to cosmologies in which Ω m is not very close to 1 at the beginning of the integration of the equations, like in Early DE models.
As we saw, the impact of homogeneous and CDE on δ c is small, see also figure 2. This suggests that the impact of DE fluctuations on structure formation occurs mainly via modifications on the matter growth function. This is somewhat unexpected because, as we will show later, at virialization, DE fluctuations can account for up to 10% of the total halo mass, whereas the change in δ c is below 1%. For possible contributions of DE fluctuations to the halo mass, see Creminelli et al. [38], Batista and Pace [49], Basse et al. [121].
This insensitivity of the collapse threshold on DE fluctuations can be understood as follows. First, δ c is defined at the collapse time, which numerically is implemented as a high value for δ NL m . Although δ NL de also grows, given the nature of the nonlinear evolution, the higher value of δ NL m will grow exponentially faster, thus making the DE contribution much less important at the collapse time. Second, the collapse threshold is given solely by the linear matter perturbation, which, although is impacted by DE perturbations, does not consider the direct contribution of DE linear perturbations.

Virialization threshold, δ v
It is important to note that, implicitly in the determination of δ c just described lies the assumption that DE fluctuations are not directly included in the quantities that define the collapse time and density threshold, namely δ NL m and δ L m . In the ΛCDM model, baryon and dark matter fluctuations are the only relevant ones, besides a small contribution of massive neutrinos 1 . Thus, the gravitational potential, which, for instance, will deflect light rays of background galaxies or set up the potential well that traps the hot intracluster gas, is entirely determined by the fluctuations in dark matter and baryon components. In the presence of DE fluctuations, the gravitational potential also depends on this new type of inhomogeneity. Therefore, it would be natural to redefine the threshold density, virial overdensity and growth function to take the contribution of DE fluctuations into account properly. For instance, equations (68) already suggests that the effective mass inside a shell of radius R includes DE fluctuations.
Then let us define effective quantities that include DE fluctuations. The total mass inside a shell of radius R is given by where M de = 4π 3 R 3ρ de δ de . Usually the background density of DE is not included in this mass definition, e.g, Creminelli et al. [48], Batista and Marra [120], Basse et al. [121]. In the case of CDE, its local EoS parameter gets less negative during the collapse. Consequently, locally, DE becomes more similar to pressureless matter, Mota and van de Bruck [36], Abramo et al. [124]. Such variations of the EoS may also be associated with soft-matter properties Saridakis [125]. The fraction of DE mass in the halo at the virialization time is Some authors have computed the values of with slightly different approaches to determine the virialization time, Creminelli et al. [38], Basse et al. [39], Batista and Pace [49], Batista and Marra [120], Heneka et al. [126]. In particular, using the approach summarized in equation (78), which assumes that CDE fluctuations behave as nonrelativistic particles, Batista and Marra [120] showed that | | can be up to 0.1, depending on w. In the case of phantom DE, this contribution is negative and positive for non-phantom. These values of raise an interesting question: if CDE can contribute up to 10% to the total halo mass, why its impact on the critical threshold is below 1%? One can also account for DE energy fluctuations in the density contrast, defined by This quantity was also used in nonlinear perturbation theory studies of CDE by Sefusatti and Vernizzi [127]. When defining the growth function as the change between clustering and homogeneous DE models with the same background can be about 3 − 7%, while the usual definition with δ m only differs about 1% Batista and Marra [120]. Moreover, this work showed that the D tot and D m (1 + ), where D m includes only the matter perturbation, are very similar, indicating consistency between the linear and nonlinear impact of CDE when using total quantities. It is important to note that, due to the DE contribution, the total mass inside the shell is not conserved during the evolution. After virialization, however, it has been argued that this contribution should be stable Creminelli et al. [38]. In Basse et al. [121], this effect was taken into account in the virial theorem for non-relativistic particles, yielding the following equation for virialization 1 2M tot Once the equations of the SCM are solved and δ NL m (z) and δ NL de (z) are known, the virialization time is determined when (78) is satisfied. For further discussion about virialization in the presence of DE, see Maor and Lahav [128], and for relativistic corrections in the virialization, Meyer et al. [129].
Having found the virialization redshift, one can compute the virialization threshold, including the contribution of DE fluctuations As can be seen in figure 3, the impact of CDE is larger than in δ c , reaching up to 4%. It is interesting to note that the magnitude of these differences is more consistent with the amount of DE fluctuations in halos described by ε than those observed in δ c . It is also important to highlight the behaviour of DE fluctuations. In figure 4, we show δ NL de (z v ). As can be seen, DE fluctuation can become nonlinear and have a mild decay at low-z, when the EoS of our examples approach −1. Clearly, the larger is |w a |, more significant are the DE fluctuations. We can also see in this plot the pathological case of δ NL de < −1 for phantom DE.
In the presence of CDE, the virialization overdensity, defined in (44), is given by In figure 5 we show the relative difference of ∆ v with respect to the ΛCDM one. Note that CDE makes ∆ v more similar to the ΛCDM values.

Halo mass functions
Now let us see how to use the findings of the SCM to estimate the observational impact of DE fluctuations on structure formation. The computation of the abundance of galaxy clusters relies on the Halo Mass Function (HMF), which gives the number density of halos per comoving volume. There are several developments for the HFM, including an analytic model based on the spherical collapse by Press and Schechter [130] and a semi-analytic approach based on ellipsoidal collapse Sheth and Tormen [131]. More recently, N-body numerical simulations were used to determine fitting functions for the HMF, e.g., Warren et al. [132], Tinker et al. [133]. Besides determining the abundance of galaxy clusters, the HMF can be used to compute the nonlinear power spectrum, Cooray and Sheth [134].
The proper approach to study the impact of CDE on the abundance of galaxy clusters would be to include its dynamical effects in the numerical simulations of structure formation. Codes with this capability have been developed only quite recently, but without determining the associated HMF. Let us first discuss developments based on analytic and semi-analytic HMF and then some relevant results from simulations.  We anticipate that there is no consensus of how to implement the impact of CDE in HMF. As a consequence, the predictions also vary, and the dection of DE fluctuations also depend on the particular implementation used to include the effects of CDE on halo abundances.
The Press-Schechter (PS) HMF Press and Schechter [130] assumes that the matter density field, δ, smoothed on some scale, has a Gaussian probability distribution The variance of the smoothed field, is given by where P(k) = |δ k | 2 is the linear matter power spectrum, is a top-hat window function with mass-scale relation given by (84) andρ m0 the background matter density now.
Assuming that regions with δ > δ c form bound structures, where δ c is a certain threshold value to be defined, the fraction of such objects with mass greater than M is given by The comoving number density of halos per mass interval is then given by where is the PS multiplicity function. The factor 2 was first introduced in an ad hoc manner, to guarantee the normalization of the HMF Press and Schechter [130], but it can be formally determined Peacock and Heavens [135], Bond et al. [136]. Assuming that only σ is mass-dependent, we have the usual form of the PS multiplicity function: In the first studies with the PS-HMF, the values of δ c assumed varied in the range (1, 10), e.g., Press and Schechter [130], Efstathiou et al. [137], Colafrancesco et al. [138], Gelb and Bertschinger [139]. Later, it became usual to use the constant EdS collapse threshold value δ c 1.69, e.g., Peacock and Heavens [135], Narayan and White [140]. In the presence of Λ (for open, flat and closed models) Lilje [141] reported 1.64 < δ c < 1.73 at low-z. As we saw in subsection (4.2), other works also determined fitting functions for these quantities for homogeneous DE models Kitayama and Suto [85], Weinberg and Kamionkowski [87], also finding a small deviation from the EdS value, of less than 1%.
The matter power spectrum can be numerically determined by codes like CAMB Lewis et al. [142] or CLASS Lesgourgues [143]. It can also be given by fitting functions, which, in most cases, can be separated in the following form where P p (k) is the primordial power spectrum associated to matter perturbations given by the inflationary model, T(k) is the transfer function (for instance, given by Eisenstein and Hu [144]) and the linear growth function of matter perturbations. Since the amount of DE is very small around the recombination time, the transfer function is essentially unaffected by DE. The main impact of DE, either smooth or inhomogeneous, is on the growth function, which strongly depends on the cosmological evolution at low-z.
In this context, for homogeneous DE models, the only relevant modification on the HMF occurs via the modifications on the growth function caused by different evolutions of w. Then it is expected that either analytic or numerical HMF can directly incorporate the impact of DE via the modifications of the growth function. Several works have studied homogeneous DE in such scenario, either with analytic approaches Percival [145], Le Delliou [146], Horellou and Berge [147], Liberato and Rosenfeld [148], Bartelmann et al. [149], Pace et al. [150,151] or numerical simulations Linder and Jenkins [152], Grossi and Springel [153].
For DE models with arbitrary c s , both δ c and D m become mass-dependent because DE fluctuations are enhanced above the sound horizon scale, and impact matter perturbations in a scale-dependent manner. This case is studied in Basse et al. [39,121,154]. Given that δ c (z) → δ c (M, z), the PS multiplicity function acquires an extra mass-dependent term Although this new mass dependence is a feature of CDE if c s is associated with some mass scale that can be probed by the observed abundance of galaxy clusters, it is much smaller than the mass dependence on σ.
In the limit of CDE, c s → 0, the growth of δ de has the same scale dependence of matter perturbations. Then both δ c and D m remain scale independent. In this case, the scale dependent feature due to CDE in the HMF vanishes, but its impact can be larger. Several papers have studied this scenario: Abramo et al. [37], Creminelli et al. [38], Batista and Pace [49], Nunes and Mota [89], Manera and Mota [92], Pace et al. [106], Batista and Marra [120], Heneka et al. [126], Abramo et al. [155]. Now let us consider the Sheth-Tormen HMF Sheth and Tormen [131] and how it can be modified to include the effects of DE fluctuations. Since it provides a better description of cluster abundances given by numerical simulations, we expect it to provide better estimates of the impact of DE fluctuations on clusters abundances. The original ST-HMF is given by  Despali et al. [81] have used several overdensities criteria to detect halos in simulations and fit the three parameters of ST-HMF. It was found that when halos are identified with the viral overdensity ∆ cc , defined in (43), the ST-HMF is nearly universal with respect to redshift and Ω m0 variations. These parameters can also be constrained by future observations, see Castro et al. [156].
Following these ideas, Batista and Marra [120] proposed that the proper threshold density for halo formation in the presence of CDE should be modified to include the contribution of DE fluctuations at virialization time, and used virial threshold, (45). But, in order to make a more conservative estimate, the usual parameter a was rescaled by a →ã = aδ 2 c /δ 2 v . With this choice, the ST-HMF is unchanged in the EdS limit. In figure 3, the evolution of δ v is shown for some homogenous and CDE models.
In practice, for a given mass scale, the abundance of massive halos strongly depends on the quantities In figure 6 we show the ratio of these functions with respect to the corresponding values in the ΛCDM model. The impact of CDE is larger in µ v and is nearly constant with redshiftdependent. Having in mind that, although δ de (z) decays at lower z, when w → −1, the quantity δ de (z)Ω de (z)/Ω m (z), which impacts both δ v and D tot , is nearly constant. Thus, the behaviour of ν v seems more consistent with the effective contribution of DE fluctuation in the dynamics of the collapse shown in figure (4). Another important point regarding the impact of CDE on HMFs is that they are not calibrated by numerical simulations. A possible approach to this issue, also used in the context of non-gaussianities LoVerde et al. [157] and baryonic feedback Velliscig et al. [158], is to multiply the more accurate numerically calibrated HMF by a factor that encodes the relative impact of DE fluctuations with respect to the usual homogeneous model Creminelli et al. [38], Batista and Marra [120], Heneka et al. [126] (dn/dM) c s =0 (dn/dM) c s =1 . (95) Moreover, one has to take into account how CDE affects the cluster mass. It is expected that the contribution of DE to mass shifts HMF according to M → M(1 − ) Creminelli et al. [38], Batista and Pace [49], Batista and Marra [120]. Batista and Marra [120] also analyzed modifications on the mass-scale relation, equation (84), and normalization condition. Whereas the former can be safely neglected, the latter is of the order of , possibly modifying HMF about a few per cent on all mass scales. These studies have shown that, depending on the evolution of w, the abundance of massive halos can change by 10 − 30% with respect to homogeneous DE models. But this change can be even larger very massive halos, M > 10 15 M , whose abundance is very sensitive to modifications on the exponential tail of HMF.

Numerical simulations
Recently, numerical simulations capable of treating CDE were developed. The first approach proposed was to include DE linear perturbations given by Einstein-Boltzmann solvers as a new source of the gravitational potential in N-Body simulations Dakin et al. [33]. Of course, this implementation can not deal with nonlinear DE fluctuations, but was crucial in confirming that linear and mildly nonlinear DE perturbations have a non-negligible impact on structure formation.
Hassani et al. [34] developed a code capable of describing nonlinear DE fluctuations. They showed that models with w = −0.9 and c s = 10 −3.5 do impact matter fluctuations and the gravitational potential on small scales. They also found that DE and matter fluctuations are correlated, as indicates the solution (14). However, HMF was not computed. Hassani et al. [35] also studied the imprint of CDE on observables associated with the gravitational potential. They confirmed the tendency of CDE to compensate for the changes in the background and found modifications of 2 − 5% on the observables studied. It is important to note that, due to the choice w = −0.9, the imprints of DE fluctuations found in these works are not as large as in models in which w is smaller at intermediate redshifts.
Although these numerical studies have confirmed several results from perturbation theory and the SCM, the distinct proposals to implement the impact of CDE on HMF were not yet tested by simulations. Unfortunately the results vary between these implementations, and there is no definitive prediction about the actual impact of CDE on the abundance of galaxy clusters. See also the discussion in Basse et al. [154]

Cosmological observables
As we saw, DE fluctuations can potentially impact the linear and nonlinear evolution of matter fluctuations and the gravitational potential. Thus, many observables like CMB anisotropies, linear and nonlinear matter power spectrum and growth rate can change due to the presence of a new clustering component. We have already discussed some of those effects, especially regarding HMF. Next, overview of some observational strategies discussed in the literature that can possibly detect CDE. The grouping of observables shown below is somewhat arbitrary because most of the works presented discuss and analyze combinations of them.

CMB and Large Scale Structure
Weller and Lewis [159] attempted to constrain the value of c s using the first WMAP data release and found that c s is unconstrained. In a similar analysis, Bean and Dore [42] reported a 1σ constraint c s < 0.2. Later, Hannestad [160] included Large Scale Structure data in the analysis, showing that c s remains essentially unconstrained. de Putter et al. [161] also reached the same conclusion, but showing that chances of detection of CDE are larger in Early DE models.
Takada [162] showed that a 2000deg 2 galaxy redshift survey at z 1, together with CMB information from Planck, can distinguish between smooth (c s = 1) and CDE with c s < 0.02 and w = −0.95.
Analyzing Early DE models and using CMB and LSS data, Bhattacharyya and Pal [163] reported c s = 0.37 and Ω de (a rec ) = 0.02. It was shown that, when c s is allowed to vary, the contribution of DE in the early Universe can be more significant.
It has also been shown that the cross-correlation between galaxy survey and the Integrated Sachs-Wolf effect is a promising technique to detect c s Hu and Scranton [164]. This idea was further explored in Ballesteros and Lesgourgues [45], Corasaniti et al. [165], Pietrobon et al. [166], Li and Xia [167].
The most recent analysis of CDE by the Planck team Ade et al. [54] indicates that c s is unconstrained. However, w was assumed constant. Since at late times we must have w −1, DE fluctuations are very small in this scenario, and the impact of c s on observables is negligible. Only models with time-varying w can present relevant DE fluctuations, like in the case o Early DE.

Higher order perturbation theory
CDE has also been studied in the framework of higher-order perturbation theory by Sefusatti and Vernizzi [127], D'Amico and Sefusatti [168], Anselmi et al. [169,170]. These works found an impact of a few per cent on nonlinear corrections to the linear power spectrum. These authors also understand that, in CDE models, the total perturbation is the relevant one to be considered, not only matter perturbations.

Weak lensing
Sapone et al. [171] investigated how tomographic weak lensing and galaxy redshift surveys can constrain DE sound speed. Considering w = −0.8 they found that, if c s < 0.01, the sound speed can be constrained. More studies in this direction include Ayaita et al. [172], Majerotto et al. [173]

Cluster abundances
Abramo et al. [155] have forecasted the constraining power of galaxy clusters surveys on c s . Although the implementation used is not entirely consistent with SCM in the sense of how c s is varied, what artificially enhances the dependence of cluster abundances on this parameter, some interesting conclusions were drawn. It was shown that future experiments like Euclid can play a decisive role in detecting DE fluctuations and that they impact the constraints on w 0 and w a at 10 − 30% level.
Basse et al. [154] has also forecasted how c s can be constrained by data from cluster abundances, CMB, cosmic shear and galaxy clustering correlations. It was found that, considering w 0 = −0.83 and w a = 0, future Euclid data can distinguish between c s = 1 and c s = 0. However, this result strongly depends on how the impact of CDE is implemented on HMF. In particular, it was found that, when considering the DE contribution to the cluster mass, the sensitivity on c s significantly degrades.
Appleby et al. [174] showed that, using the Euclid satellite cluster survey together with Planck data, Early DE models with c s < 0.1 and Ω de (a rec ) = 0.009 can be distinguished from models with c s = 1. In such models, DE is not negligible at high-z and the EoS parameter can be larger than −1 at high-z. Hence, the impact of DE perturbations is strongly enhanced.
The authors noted that c s is mainly constrained by CMB data. Although this work used the Tinker HMF Tinker et al. [133], which does not take into account the nonlinear effects of DE fluctuations, the reported insensitiveness of cluster abundance on clustering Early DE model is in accordance with Batista and Pace [49]. The reason for this is the following: whereas Early DE has a strong impact on the background evolution, decreasing the matter growth, if it clusters, the perturbations partially compensate for this change. Note, however, that this result was obtained using δ c and D m . When using δ v and D tot in HMF, it is expected that more pronounced changes due to CDE will be present.
Heneka et al. [126] used cluster data from several experiments plus CMB, Barion Acoustic Oscillations (BAO) and Supernova Ia (SNIa) observations to constraint cosmological parameters with c s = 1, c s = 0 and constant w. It was shown that the allowed regions in the σ 8 − Ω m0 plane change due to CDE. It was found that w < −1 is preferred by data and that σ 8 is reduced when c s = 0. This impact of CDE can alleviate the tension in cluster data found by Planck Ade et al. [175].
It is important to note that, a possible issue regarding constraints of CDE models with cluster data is associated with the observable-mass scaling relations, see Mantz et al. [176], Kravtsov and Borgani [177]. So far, no analytical or numerical study has been conducted to explore the impact of DE fluctuation on these relations. As discussed before, DE fluctuations certainly impact the gravitational potential, thus lensing signals, X-ray can SZ luminosities relations with the cluster mass can be affected.

Internal structure of galaxy clusters
The impact of CDE on the internal structure of halos was studied in Mota [178], Basilakos et al. [179]. As it would be expected from the results that R ta /R v < 0.5 in the presence of DE and its possible perturbations, the concentration parameter of galaxies clusters increases. Basilakos et al. [179] have analyzed the mass-concentration relation of four massive clusters, showing that the data is better described by CDE models.

S 8 tension and growth rate
Clustering DE models can be used to reduce the value of σ 8 , Kunz et al. [72]. The parameter S 8 = σ 8 √ Ω m0 /0.3 inferred by CMB observations using the ΛCDM model, Aghanim et al. [3], is about 3σ in tension with the values determined by weak lensing observations Abbott et al. [4], Asgari et al. [5], see also Di Valentino et al. [12] and references therein. To explore the impact of CDE -here we will focus on the variation of σ 8 only -we fix σ 8mod (z) = σ 8mod D mod (z), for a given model "mod", to have the same value as in the ΛCDM model at z rec where σ 8Λ = 0.8111 for Ω m0 = 0.3153 Aghanim et al. [3]. As can be seen in the left panel of figure (7), in the case where only matter perturbations contribute to the growth function, CDE make the evolution of σ 8 (z) more similar to ΛCDM. In the right panel, when DE perturbations are included in the growth function, we see that the impact of CDE is much larger. In this case, for w a > 0 (w a < 0), DE perturbations can make σ 8 (z) larger (smaller) than in ΛCDM.
Homogeneous phantom DE increase σ 8 because Ω de (z) grows rapidly at low-z, delaying the suppression of matter perturbations when compared to ΛCDM. When DE perturbations are allowed, given that δ de < 0, matter growth is suppressed. Considering the total growth function, this effect can potentially alleviate the S 8 tension.
On the other hand, homogeneous non-phantom DE decrease σ 8 because Ω de (z) starts to grow earlier than in ΛCDM model. Clustering DE now partially compensates for this decay in the growth of matter perturbations, enhancing σ 8 . For the case of the total growth function, CDE can potentially worsen the S 8 tension. An analysis of CDE using cluster, CMB, BAO and SNIa data was done by Heneka et al. [126]. It was found the σ 8 is reduced for c s = 0 with slight increase in Ω m0 . Considering the best fit values, the overall effect is a reduction of S 8 in CDE models. The results also show that, assuming constant w, phantom values are prefered by data.
Clustering DE also affects the growth rate of matter Batista [180], Mehrabi et al. [181,182], given by In particular, Mehrabi et al. [183] reported that CDE models can fit f (z)σ 8 (z) data better than ΛCDM.

Discussion
The variety of models that try to explain the cosmic acceleration is enormous. The most studied and tested ones either do not present DE fluctuations (Λ) or have negligible fluctuations on scales well below the horizon (quintessence). However, many other models can present a low c s value and relevant DE fluctuations on small scales. However, the actual amount of DE fluctuations also depends on the evolution of EoS. If w −1 throughout the cosmic evolution, DE fluctuations are very small, regardless of the value of c s .
Hence, the prospects to detect DE fluctuations on small scales, which would rule out Λ and quintessence as possible drivers of cosmic acceleration, strongly depend on how far from −1 the EoS can be in the past. Although Λ is in excellent agreement with almost all cosmological data available, there are no strong constraints on how much w can deviate from −1 at intermediate redshifts. For instance, the allowed parameter space for w 0 and w a for the case of c s = 1 is still very large Ade et al. [54]. As we reviewed, for constant EoS, DE models with w −0.9 and c s < 10 −2 leave a detectable impact on several cosmological observables. Being more conservative with the value of w now, we showed that, even with w 0 = −1, DE fluctuations become nonlinear when assuming |w a | > 0.2 and c s = 0, and influence structure formation.
The linear evolution and the corresponding impact of DE perturbations is well understood and implemented in numerical codes like CAMB and CLASS. The nonlinear evolution in the presence of CDE is, however, more complicated and much less developed. First efforts in this direction used the SCM and only recently numerical codes for the nonlinear evolution were developed. The early findings about the phenomenology of CDE are in agreement with those found by recent numerical results. These studies have shown that DE fluctuations can become nonlinear and do impact matter fluctuations, the formation of halos and the gravitational potential.
The abundance of massive galaxy clusters is certainly one of the most affected observables. However, up to now, no numerical simulations modelled the impact of CDE in HMF. The analytical motivated proposals for HMF in the presence of CDM are divided into two main groups: one in which the impact of CDE enters only via the matter quantities (δ c , D m , σ 8m ) and another based on the total weighted fluctuations (δ v , D tot , σ 8tot ). Naturally, the impact of CDE is enhanced in the latter type of modifications. Moreover, these different recipes can be combined with mass rescalings due to DE contribution, M → M(1 + ). However, there is no consensus on the literature on which is the accurate implementation to compute the effect of CDE on halo abundances.
An important theoretical issue with CDE models is associated with phantom models. For c s → 0, one can have δ de < −1 in the nonlinear regime. At first glance, this indicates that these models are inconsistent. A possible realization of phantom clustering models may include imperfect fluids Sawicki et al. [184].
Although CDE is the natural minimal generalization to quintessence, it is still not well explored in the literature. Studies which consider parametrizations of w are very common, but rarely consider c s < 1. As discussed, CDE can play a role in alleviating the S 8 tension. From another perspective, there is no indication that DE can not present relevant fluctuations on small scales. Therefore, these kinds of models deserve more attention from the cosmology community.
Besides these difficulties, preliminary forecasts studies indicate that future experiments like Euclid will be able to discriminate between homogeneous and CDE. The challenges to achieve this goal are big. Most of the constraining power of Euclid will come from nonlinear scales, and a good understanding of the nonlinear effects of CDE will be mandatory.