# Dynamical evolution of quintessence dark energy in collapsing dark matter halos

###### Abstract

In this paper, we analyze the dynamical evolution of quintessence dark energy induced by the collapse of dark matter halos. Different from other previous studies, we develop a numerical strategy which allows us to calculate the dark energy evolution for the entire history of the spherical collapse of dark matter halos, without the need of separate treatments for linear, quasi-linear and nonlinear stages of the halo formation. It is found that the dark energy perturbations evolve with redshifts, and their specific behaviors depend on the quintessence potential as well as the collapsing process. The overall energy density perturbation is at the level of for cluster-sized halos. The perturbation amplitude decreases with the decrease of the halo mass. At a given redshift, the dark energy perturbation changes with the radius to the halo center, and can be either positive or negative depending on the contrast of , and with respect to the background, where is the quintessence field. For shells where the contrast of is dominant, the dark energy perturbation is positive and can be as high as about .

###### pacs:

98.80.-k, 95.36.+x## I Introduction

The fast advance of different cosmological observations permits us to weigh the universe quantitatively Astier et al. (2006); Riess et al. (2007). It is known that the universe is currently dominated by the dark energy component, which accounts for about of the total energy density. Probing the nature of dark energy is thus one of the most challenging tasks in cosmological studies. The simplest dark energy candidate is the cosmological constant with the equation of state . While being consistent with observations, it runs into severe difficulties theoretically, namely the fine tuning problem and the coincidence problem. Various dynamical dark energy models have been proposed aiming to avoid such problems Copeland et al. (2006); Frieman et al. (2008); Linder (2008). The allowed parameter space for dynamical dark energy models derived from current observational data is still relatively large Komatsu et al. (2009); Li et al. (2008a). Tighter constraints are highly expected from future observational data with much improved quality and quantity. To realize the goal, however, it is crucial to understand in detail the observational impacts from the dark energy component. Different dark energy models give rise to different expansion histories of the universe, which in turn influence differently on the formation and evolution of large-scale structures Klypin et al. (2003); Macciò (2005); LeDelliou (2006); Mainini and Bonometto (2007). Besides, clustering can occur for dynamical dark energies themselves, and the clustering characteristics are highly model dependent. Therefore the observable imprints from dark energy perturbations can potentially affect the constraints on the properties of dark energy.

Dark energy perturbations associated with the linear evolution of the matter inhomogeneities have been studied extensively Unnikrishnan et al. (2008); Amendola (2004). Considerable effects can occur. For instance, it has been found that with the proper inclusion of dark energy perturbations, the constraints on the time-evolving equation of state of dark energy are weakened significantly due to their influence on the late ISW effects of CMB anisotropy Zhao et al. (2005, 2007); Li et al. (2008b). On the other hand, the dark energy clustering related to the nonlinear structure formation is less addressed. This is partially due to the numerical difficulties in solving the coupled perturbation equations for matter and dark energy in the nonlinear regime. However, many cosmological observations, such as the abundance of clusters of galaxies and weak gravitational lensing effects on small scales, are directly related to the nonlinear structure formation LaVacca and Colombo (2008). It has been proposed that the recently reported cold spots seen on the CMB map may arise from the Rees-Sciama effect, the nonlinear version of the ISW effects Granett et al. (2008); Sakai and Inoue (2008); Masina and Notari (2009). Therefore there is a need in understanding thoroughly the effects of dark energy on the formation of nonlinear structures. It is also of interest to investigate the corresponding clustering behavior of dark energy.

There have been analyses on the formation of nonlinear structures, such as the virial radius and the average mass density of halos, taking into account the existence of the background dark energy Wang and Steinhardt (1998); Maor and Lahav (2005); Mota and van de Bruck (2004); Nunes and Mota (2006). Parameterized dark energy perturbations and their correlations with the structure formation are investigated in some of the studies Langlois and Vernizzi (2007); Abramo et al. (2007); Maor and Lahav (2005); Mota (2008). Recently, under the framework of spherical collapse of halos, the perturbations of quintessence fields induced by the structure formation are analyzed Dutta and Maor (2007); Mota et al. (2008). Different from the parameterization approach, in which the equation of state and the perturbation of dark energy are regarded as being independent with each other, such analyses can present a consistent picture on the dark energy perturbation that often correlates with the evolution of the equation of state strongly. Dutta and Maor Dutta and Maor (2007) mainly study the dark energy perturbation in the linear regime of the halo formation. Mota, Shaw Silk Mota et al. (2008) consider the dark energy clustering in the process of nonlinear structure formation employing a matched asymptotic expansion analysis method, in which, they adopt different parameter expansions for the linear, quasi -linear and nonlinear stages, respectively. In this paper, we present a numerical method within the spherical collapse model that allows us to calculate the co-evolution of the halo formation and the dark energy clustering during the whole process, from linear to nonlinear epoches. In this approach, we do not need to treat the linear, quasi-linear and nonlinear evolutions separately. We analyze the quintessence dark energy perturbation paying particular attention to its correlation with the background field evolution, i.e., the evolution of the equation of state of dark energy.

The paper is organized as follows. In section 2, we present the model and the corresponding formulae. Section 3 explains in detail our numerical approach. Our results are shown in section 4. Section 5 contains discussions.

## Ii The Model and Formulae

### ii.1 Two-components spherical collapse model

In the framework of the spherical collapse of dark matter halos, we consider the induced perturbation of dark energy. We start with the line element of a spherically symmetric system, which can be written as

(1) |

where and R are functions of both time and radius . Note the corresponding in the Friedmann-Robertson-Walker (FRW) metric, where is the scale factor of the universe. The coordinate system is chosen such that it is comoving with the dark matter component. Thus the component of the energy-momentum tensor for dark matter is equal to the matter density, and the other components are zero. For the quintessence model, the energy-momentum tensor is

(2) |

and

(3) |

where V() is the potential of the quintessence. The total energy-momentum can therefore be expressed as

(4a) | |||

(4b) | |||

(4c) | |||

(4d) | |||

(4e) |

where the dot denotes the time derivative and the prime denotes the spatial derivative .

From the Einstein field equations, we have

(5a) | |||||

(5b) | |||||

(5c) |

The evolution of the quintessence scalar field, which is not coupled to the dark matter component except gravitationally, is governed by

(6) |

where .

### ii.2 Tolman approximation

The perturbation of non-coupled quintessence dark energy is expected to be small. Thus to a good approximation, we can ignore the time-spatial term in the right hand side of the equation (5c). By time integrating its left hand side, we obtain a relation between and R, which follows

(7) |

where is a constant of integration and a function of only, which reflects the local spatial curvature. Then we have

(8) |

This is the well-known Tolman metric that is often utilized to compute the evolution of spherical objects. For a homogeneous universe, the Tolman metric reduces to the FRW metric and . With the Tolman metric, the Einstein equations become

(9a) | |||

and | |||

(9b) |

### ii.3 Spatial Curvature

In the Tolman metric approximation, the term does not depend on time. Therefore we can calculate this quantity at any chosen time. In our consideration, we calculate at the initial redshift chosen to be . Because the dark energy component is negligible at such a high redshift, the equations (9a) and (9b) can be further simplified. By combining the two equations, we obtain

(11) |

where denotes the mass enclosed within the comoving radius r, and

(12) |

Therefore depends on the initial dark matter density profile. In this paper, we adopt the following density distribution

(13) |

where is the initial background dark matter density, is the initial central density contrast in the halo region, and is related to the characteristic scale of the halo and will be discussed further later. For /, such a perturbation profile has an overdense/underdense zone in the central region and an underdense/overdense transitional middle region which compensates for the central region so that the overall perturbation is zero.

### ii.4 The double exponential quintessence model

In this paper, we consider the quintessence model with a commonly adopted double exponential potential

(14) |

where and are parameters representing the profiles of the two exponentials, respectively, and reflects the height of the potential. In our analyses, we take and as our fiducial model. We also vary with to see the effects of the potential on dark energy perturbations. The starting redshift is taken to be . The initial velocity of the field is taken to be . The parameter and the initial value of are adjusted so that we have and at . Specifically, we first choose a pair of reasonable values of and the initial , and evolve the background field to the present to calculate and at . By comparing the obtained and with and , we adjust the values of and the initial to compute the evolution of the field again. Such a process is repeated until and are reached. It is known that the double exponential quintessence model has a late-time attractor solution, in which the final and are mostly determined by the potential, i.e., the values of , and , but not very sensitive to the initial conditions of Barreiro et al. (2000). Indeed, such a tracking behavior has been seen in our analyses. On the other hand, we start from a relatively low redshift (), and choose initially. Thus albeit not an extremely fine-tuning process, a certain level of adjustment on the initial value of is needed in our analyses.

Including the perturbations, we have

(15a) | |||

(15b) | |||

and | |||

(15c) |

In Fig. 1, we show the evolution of the background equation of state with for , respectively. The corresponding evolution of the time derivative of the background quintessence field, , are shown in Fig. 2. It is seen that with the increase of , the potential gets steeper and the term is larger. As a result, the equation of state increases. For and , in the redshift ranges and , respectively. It will be seen later that the evolutionary behavior of the background field affects the dark energy perturbations significantly, thus in general, the two cannot be treated independently.

## Iii Numerical method

### iii.1 Iterative algorithm

We first compute the local metric in the region of the spherical dark matter halo. We incorporate the local dark matter density and the dynamical evolution of the background quintessence field into our calculation. The dark energy perturbations are neglected in the metric calculations. Thus equation(9a) can be written as follows

(16) |

where stands for background. Equation (9b) has been used in deriving the expression of given in part C of the previous section.

Equation (16) is an ordinary differential equation. We then solve it at different shells equally spaced in . The initial conditions for are chosen to be and at , where is the scale factor of the universe at .

With the obtained metric, we thus can analyze the spatial perturbation of dark energy. We rewrite the equation(II.2) in the following form

(17) | |||||

Equation (17) can be regarded as a differential equation of with respect to with the right hand side as the source of the spatial perturbations of . We solve this equation iteratively at each time . First, we insert the background and into the right hand side of equation (17). With this source term, we solve at each . Then for each obtained , we find its corresponding . This is done by looking through the evolution of the background field , and locating its value at the place where . Then at this time, different value is obtained at different . In the next iteration, this set of are used in calculating the source term, which lead to a new set of . We repeat such an iteration process until the solutions on converge. For each iteration, the shooting method from the center outward is applied to solve this spatial differential equation Press et al. (2002). The boundary conditions are taken to be and . As shown in Eq. (13), we consider compensated matter density perturbations. Thus is defined to be the shell inside which the average matter density perturbation is close to zero. The evolution of the shell at follows that of the background universe.

In comparison with the method solving the partial differential equation directly, our numerical algorithm not only saves time but also avoids propagating numerical errors step by step in the time advancement.

### iii.2 Initial conditions

In this study, we consider flat cosmologies with the parameters being taken from WMAP5 results Dunkley et al. (2009). Specifically, , and , where is the dimensionless matter density of the universe, is the Hubble constant in units of , and is the linearly extrapolated matter density fluctuation on top-hat scale of . The initial redshift is taken to be . We consider overdense regions with typical scales of clusters and galaxies, respectively. For galaxies, we take , and for clusters of galaxies, we take . The corresponding linear comoving scales are calculated through , where is the critical density of the universe at present. This gives rise to and , respectively. The initial perturbations over the considered scales are respectively set to be and , estimated by employing the linear power spectrum of Eisenstein & Hu Eisenstein and Hu (1998).

We now derive the relation between and in equation (13) and and . From the definition of , we have

where

(18) |

Then we obtain

In our analyses, we take , therefore we have . The largest radius in our calculation is set to be for the galaxy scale and for the cluster scale. Beyond the largest radius, all the relevant quantities merge into the background ones.

For the curvature term , as we discussed previously, it can be calculated at . Thus we have

(19) |

We normalize with the largest scales in our calculations.

### iii.3 Virialization

In the general spherical collapse model, the virialization process is not naturally included. It is known that the violent realization and phase mixing play important roles in the final stage of the formation of dark matter halos. However, detailed understanding and modelling of the virialization process are still lacking. In Shaw and Mota (2008); Mota et al. (2008), they fit the simulation results to describe the virialization stage.

Here we adopt the following description for the virialization. The standard spherical model is used to calculate the metric evolution of for each shell before its turn-around, denoted as . After the turn-around, the time derivative of is modified such that it is approaching zero at the virialization radius . Specifically, we use

(20) |

where the stands for the standard spherical collapse model. It is noted that in our modelling, the overdense region has a r-dependent density profile. Consequently, different shells reach turn-around at different time. Thus Eq. (20) models the virialization for regions inside individual shells. In other words, the virialization starts from the inner region and gradually extends outward. To certain extents, we believe that our modelling grabs the basic picture of the virialization.

Fig. 3 shows for some shells, where we take instead of 2.0 from simulation results used in Mota et al. (2008); Hamilton et al. (1991). We also calculate the mean overdensity , where is the mean mass density within and the mass contained inside is computed by Eq. (12), and is the background density at the virialization time. Here the virialization redshift for a shell is defined to be the redshift when the shell collapses to the center in the case without considering the virialization treatment. For different shells, the results are shown in Table 1. It is seen that for regions where the virialization happens at high redshift, , the well known value for . At later times, as the contribution from the dark energy component becomes significant, decreases. At , we have . Our results are consistent with those from the simple consideration of energy partitions Mota (2008), indicating the meaningfulness of our treatment for the virialization.

In Fig. 4, we show the time evolution of the relative density profile, where the horizontal axis is with being the value at and the verticle axis is the relative density . As expected, the relative density gets higher and more centrally concentrated as time evolves. On the other hand, it is noted that our prescription of the virialization Eq. (20) only accounts for the fact that for a virialized shell. It does not include considerations about the detailed virialization processes, such as shell crossing, phase mixing and violent relaxation. Thus in our treatment, the evolution of is basically an initial value problem, and the final relative density profile is closely related to its initial one. Since the NFW profile found from numerical simulations Navarro et al. (1996)Navarro et al. (1997) must be at least partially related to the virialization processes, we do not expect that the NFW profile can be naturally obtained in our simple analyses presented here without adjusting the initial relative density profile carefully.

## Iv Result

In this section, we present the results on the calculated dark energy perturbations.

Fig. 5 shows the time evolution of for four spatial shells with , , and , respectively, for the cluster-scale halo, where is the dark energy density perturbation, and is in unit of the largest radius ( in this case) in our calculation. The corresponding spatial variations of at different redshifts are shown in Fig. 6. For the purpose of easy comparison with the results of Mota et al. (2008), we follow them to show instead of . Under certain approximations, it is shown in Mota et al. (2008) that the quantity can be directly related to the matter density perturbations. This motivates them to focus on . In our following presentations, besides , we also show at various places. We can see that the amplitude of is except for the central region where can reach as high as about . At , the dark energy perturbations are positive everywhere in the calculated region, showing an initial concentration of dark energy in the region. With the evolution of the halo formation, the outer region appears as a dark energy void with negative , while the inner region (except for ) has a positive dark energy perturbation. As time proceeds, the void region expands till very late time in the future when the dark energy perturbation turns to be positive again everywhere in the region.

In the following, we will analyze different contributions to the dark energy density perturbations in detail. Generally, there are three terms to determine the total energy density of a scalar field, namely, the kinetic energy term , the potential energy term and the term related to the spatial variation of the field . For the very central shell with , we choose a boundary condition with , and thus for this shell, only the kinetic and the potential terms exist. In Fig. 7, the upper panel shows the perturbations of the field (dotted line) and its time derivative (dashed line) with respect to the background quantities. Because of its slower expansion rate than the background, the dragging term is smaller in equation (II.2). Therefore the evolution of the field is faster, which leads to smaller value of and the corresponding change of comparing to those of the background quantities. The middle panel shows the corresponding energy perturbations. The lower panel presents , and in logarithm scale, where the thick and thin lines correspond to positive and negative quantities, respectively. It is seen that the potential energy perturbation is always negative, while the kinetic energy perturbation can be positive or negative depending on the shape of the potential. At early times with , the kinetic energy perturbation is dominant, and thus we have . At later times, the potential energy perturbation becomes the dominant term, which gives rise to the negative dark energy perturbation, i.e., the dark energy void. It is noted that as the shell turns over, the term changes sign, and thus it is no longer the friction term, but a driving force for the local evolution of the field. It causes the deviation of from the background, and eventually leads to the domination of the kinetic energy perturbation and thus positive dark energy perturbations again at very late time with . It can be seen that the behavior of the dark energy perturbations is closely related to the form of the potential of the field, which in turn determines the evolution of equation of state of the background dark energy. Fig. 8 shows the time evolution of the dark energy perturbations of the central shell for different potentials with (dotted line), (solid line) and (dashed line), respectively. The dependence of on is clearly seen. Therefore the dark energy perturbations and its background E.o.S are in general highly correlated, and cannot be treated independently.

For shells between the center and the outer boundary, the spatial variation of the field contributes significantly to the dark energy perturbations. Because of the homogeneity of the background, the spatial perturbation term always contributes positively to dark energy fluctuations. In Fig. 9, we show different parts of dark energy perturbations for the shell with . The upper panel shows the perturbations of (dotted line) and (dashed line), respectively. In comparison with those of (the upper panel of Fig. 7), these perturbations are smaller. The middle panel shows the potential (dotted line), kinetic (dashed line), spatial (dash-dotted) terms, and the total (solid line) of dark energy perturbations. For this shell, the spatial term is the major perturbation term over a large redshift range with . The lower panel shows the corresponding quantities over in log-scale. It is seen that for , the spatial term is much larger than the other two terms, and can be as large as about . It is also noted that and the perturbation of are correlated with each other.

From Fig. 6, we see that the peak perturbation occurs at at high redshifts (solid and dashed lines). It gradually moves outward, and eventually reaches . We find that the peak position is in good accordance with the shell within which the virialization happens. Such a shell corresponds to the position of maximum .

Fig. 10 shows the average dark energy density perturbations within (solid line), (dashed line), (dotted line) and (dash-dotted line), respectively. Firstly, we see the notable high peak at for the solid curve, which is contributed by the spatial derivative term . As we average over larger regions outward, the dark energy void balances its concentration in the inner region and results smaller amplitude of dark energy perturbations. The inserted region of the plot shows the zoom-in pattern of the corresponding perturbations in very late times. The spatial variation of the dark energy perturbations can be relevant in considering the ISW effect in a collapsing overdense region, for which we will investigate further in our future studies.

To see the dependence of our results on the chosen initial density profile, we also perform the analyses with a compensated top-hat initial density distribution. The mass contained in the overdense region and the central density perturbation are set to be the same as those of the compensated Gaussian profile of Eq. (18). Thus in the top-hat profile, the overdense region is narrower than that of the Gaussian one. Fig. 11 shows the spatial distribution of at with the solid and the dashed lines for the top-hat and the Gaussian profiles, respectively. It is seen that the overall behaviors of the two are similar, with positive dark energy perturbations in the inner regions and negative ones in the outer regions. The narrower range of the positive perturbations in the top-hat case is closely related to its narrower matter overdense region than that of the Gaussian case. Since the profile of the initial matter density perturbations affects the detailed collapsing process, the specifics of dark energy perturbations are expected to depend on the initial matter density profile. However, the overall level of the dark energy perturbations is anticipated to be most sensitive to the total mass contained in the overdense region. This can be seen from the following discussions. On the other hand, it is worthwhile to investigate in detail the dependence of the dark energy perturbations on the initial mass density profile, another direction that we can pursue in the future.

We now compare dark energy perturbations induced by the formation of different sized halos. Fig. 12 shows the dark energy perturbations at of the galaxy-sized halo with (solid line) and the cluster-sized halo with (dashed line). It is seen that the first transition from positive to negative dark energy perturbations occur at for both cases. As we discussed previously, this transition reflects the transition of the field from the fast evolutionary stage with a steep potential to the slow evolutionary stage with a relatively flat potential. Thus it is mainly determined by the behavior of the field potential itself, which explains the similarity seen in the two cases. On the other hand, the second transition from negative to positive perturbations in late times is closely related to the formation of the halos. It happens earlier for galaxies than that of clusters. This is because of the higher density perturbations initially set for the galaxy-sized halo, which leads to its earlier formation than that of the cluster. For the overall amplitude of dark energy perturbations, it is about an order of magnitude larger for the cluster case than that of the galaxy case. Thus it is mainly the mass of the system rather than the initial amplitude of the dark matter density perturbation that determines the overall amplitude of the quintessence dark energy perturbation. This is consistent with the expectation that quintessence dark energy perturbations show their effects mainly on very large scales.

## V Discussion and Conclusion

In this paper, we analyze the quintessence dark energy perturbations induced by the spherical collapse of dark matter halos. We develop an iterative algorithm which allows us to solve the field perturbations consistently from linear to quasi-linear and to nonlinear stages of the halo formation. We find that the dark energy perturbations depend sensitively on the field potential, and thus on the background equation of state. In the linear regime of the halo formation, the slower expansion near the overdense region leads to the local evolution of the field to be ahead of the background field. When the quintessence field is in its fast evolutionary stage, is accelerating quickly. The perturbation of the kinetic dark energy term is positive, and dominates over the negative perturbation of the potential energy term. Adding in the positive perturbation from the spatial term results a dark energy enhancement everywhere around an overdense region at the very early stage of the evolution. As the field enters the region where the potential gets flatter, the term is smaller, and the potential dark energy perturbation takes over. This gives rise to a negative perturbation, namely, dark energy void, in central and outer regions where the spatial variation term is small in accordance with our chosen boundary conditions. On the other hand, for regions in between, the term is relatively large and there are strong dark energy concentrations. In the late stage when the considered region collapses and virializes, the friction term changes sign to become a driving source for the field evolution. Thus the time evolution of is accelerated although the potential is nearly flat. This leads to the domination of the positive kinetic dark energy perturbation again in very late stages. The overall dark energy perturbation is at the level for cluster-sized halos and for galactic halos.

The work by Dutta and Maor (2007) is the earliest one to analyze the dynamical dark energy evolution in the framework of spherical collapse of dark matter halos. They employ a linearized numerical approach, and thus only consider the linear stage of the halo formation. They find a scalar field void associate with the overdense region, and further argue that it should be deepened in the nonlinear stage. While our results from the linear stage of the halo formation are in qualitative agreement with theirs, some differences are worth mentioning. First, in their calculation of the dark energy perturbations, the term associated with the spatial variation of the field is not included. It is seen from our analyses that in some spatial regions, this term can contribute significantly to positive dark energy perturbations to form a dark energy concentration rather than a void. Secondly, our results (e.g., Fig. 7) show that the field void is indeed deepened initially as expected, but then gradually becomes shallower. This behavior can be understood by noting that in the regions where the potential of the field is relatively flat, the velocity of the field decreases. As a result, the background field gradually catches up the field around the overdense region, and decreases. Therefore the depth that the field void can go depends sensitively on the dark energy potential, and thus the background equation of state.

Our study is mostly related to the study in Mota et al. (2008). In their analyses, they adopt the method of matched asymptotic expansions , in which different expansion parameters are used at different regions, mainly in the interior region and in the exterior region. The two are required to match with each other in the intermediate region. With such approximations, they derive a relation between the dark energy perturbations and the peculiar velocity and the average matter density contrast . Different treatments for the linear, quasi-linear and nonlinear evolution of the halo formation are then applied to obtain the information of and , which finally give rise to the dark energy perturbations. On the other hand, we employ a full numerical treatment. The metric evolution is calculated in the Tolman approximation with the effect of dark energy perturbations neglected. We introduce an analytical description to mimic the virialization of dark matter halos, which allows to carry out the metric calculation all the way from linear to nonlinear stages of the halo formation consistently. With the obtained metric, an iterative algorithm is used to compute the dark energy perturbations. Given different approaches, our results are in good accordance with those of Mota et. al. Mota et al. (2008). Both have seen the transition for the dark energy perturbation from positive to negative and to positive again. However, some quantitative differences exist. Considering the central region, the second transition from negative to positive dark energy perturbations happens later from our results. We notice that in Mota et. al. Mota et al. (2008) the metric calculation in nonlinear stages is done by ignoring the dark energy component completely. In our analyses, the dark energy perturbations are neglected, but the background dark energy component is kept in computing the metric evolution.

For a test, we have done a calculation on metric by dropping the dark energy component. The result on dark energy perturbations is shown in Fig. 13 for the central region (dashed line) together with the result with background dark energy component included in metric calculations (solid line). It is seen that the second transition indeed occurs earlier for the dashed line. This basically reflects that the existence of background dark energy affects the epoch of halo formation, which in turn affects the behavior of dark energy perturbations. We also find that while both catch the contribution of the term to dark energy perturbations, we obtain a larger amplitude of the term than that of Mota et. al. Mota et al. (2008). This is closely related to the different initial density profiles as well as the different descriptions of the virialization process in the two studies. This, to certain extents, indicates that dark energy perturbations are sensitive to the details of their perturbing matter distributions.

Our algorithm can be readily applied to study dark energy perturbations of different models. Further investigations on their cosmological implications and seeking for possible observational effects of dark energy perturbations, e.g., the ISW effects around superclusters and clusters of galaxies, are highly desired and will be carried out as our future tasks.

###### Acknowledgements.

We are very grateful to the referee for the encouraging comments and suggestions. We thank Xinmin Zhang, Hong Li, Junqing Xia, and Mingzhe Li for their helpful discussions. This research is supported in part by the NSFC of China under grants 10373001, 10533010 and 10773001, and the 973 program No.2007CB815401.## References

- Astier et al. (2006) P. Astier, J. Guy, N. Regnault, R. Pain, E. Aubourg, D. Balam, S. Basa, R. G. Carlberg, S. Fabbro, D. Fouchez, et al., Astron. Astrophys. 447, 31 (2006).
- Riess et al. (2007) A. G. Riess, L.-G. Strolger, S. Casertano, H. C. Ferguson, B. Mobasher, B. Gold, P. J. Challis, A. V. Filippenko, S. Jha, W. Li, et al., Astrophys. J. 659, 98R (2007).
- Copeland et al. (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006).
- Frieman et al. (2008) J. Frieman, M. Turner, and D. Huterer, Annual Review of Astronomy and Astrophysics 46, 485 (2008).
- Linder (2008) E. V. Linder, Gen. Relativ. Gravit. 40, 329 (2008).
- Komatsu et al. (2009) E. Komatsu, J. Dunkley, M. R. Nolta, C. L. Bennett, B. Gold, G. Hinshaw, Jarosik, D. Larson, M. Limon, L. Page, et al., Astrophys. J. S. 180, 330 (2009).
- Li et al. (2008a) H. Li, J.-Q. Xia, Z. Fan, and X. Zhang, J. Cosmol. Astropart. Phys. 10, 046 (2008a).
- Klypin et al. (2003) A. Klypin, A. V. Macciò, R. Mainini, and A. Bonometto, Astrophys. J. 599, 31 (2003).
- Macciò (2005) A. V. Macciò, Mon. Not. Roy. Astron. Soc. 361, 1256 (2005).
- LeDelliou (2006) M. LeDelliou, J. Cosmol. Astropart. Phys. 01, 021 (2006).
- Mainini and Bonometto (2007) R. Mainini and S. Bonometto, J. Cosmol. Astropart. Phys. 09, 017 (2007).
- Unnikrishnan et al. (2008) S. Unnikrishnan, H. K. Jassal, and T. R. Seshadri, Phys. Rev. D. 78, 123504 (2008).
- Amendola (2004) L. Amendola, Phys. Rev. D. 69, 103524 (2004).
- Zhao et al. (2005) G.-B. Zhao, J.-Q. Xia, M. Li, B. Feng, and X. Zhang, Phys. Rev. D. 72, 123515 (2005).
- Zhao et al. (2007) G.-B. Zhao, J.-Q. Xia, B. Feng, and X. Zhang, Int. J. Mod. Phys. D 16, 1229 (2007).
- Li et al. (2008b) H. Li, J.-Q. Xia, G.-B. Zhao, Z.-H. Fan, and X. Zhan, Astrophys. J. 683, L1 (2008b).
- LaVacca and Colombo (2008) G. LaVacca and L. P. L. Colombo, J. Cosmol. Astropart. Phys. 04, 007 (2008).
- Granett et al. (2008) B. R. Granett, M. C. Neyrinck, and I. Szapudi, Astrophys. J. 683, 99 (2008).
- Sakai and Inoue (2008) N. Sakai and K. T. Inoue, Phys. Rev. D. 78, 063510 (2008).
- Masina and Notari (2009) I. Masina and A. Notari, J. Cosmol. Astropart. Phys. 02, 019 (2009).
- Wang and Steinhardt (1998) L. M. Wang and P. J. Steinhardt, Astrophys. J. 508, 483 (1998).
- Maor and Lahav (2005) I. Maor and O. Lahav, JCAP 0507, 003 (2005).
- Mota and van de Bruck (2004) D. F. Mota and C. van de Bruck, Astron. Astrophys. 421, 71 (2004).
- Nunes and Mota (2006) N. J. Nunes and D. F. Mota, Mon. Not. Roy. Astron. Soc. 368, 751 (2006).
- Langlois and Vernizzi (2007) D. Langlois and F. Vernizzi, J. Cosmol. Astropart. Phys. 02, 017 (2007).
- Abramo et al. (2007) L. R. Abramo, R. C. Batista, L. Liberato, and R. Rosenfeld, J. Cosmol. Astropart. Phys. 11, 012 (2007).
- Mota (2008) D. F. Mota, J. Cosmol. Astropart. Phys. 09, 006 (2008).
- Dutta and Maor (2007) S. Dutta and I. Maor, Phys. Rev. D 75, 063507 (2007).
- Mota et al. (2008) D. F. Mota, D. J. Shaw, and J. Silk, Astrophys. J. 675, 29M (2008).
- Barreiro et al. (2000) T. Barreiro, E. J. Copeland, and N.J.Nunes, Phys. Rev. D 61, 127301 (2000).
- Press et al. (2002) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C++ (Cambridge University Press, 2002).
- Dunkley et al. (2009) J. Dunkley, E. Komatsu, M. R. Nolta, D. N. Spergel, Larson, G. Hinshaw, L. Page, C. L. Bennett, B. Gold, N. Jarosik, et al., Astrophys J. S. 180, 306 (2009).
- Eisenstein and Hu (1998) D. J. Eisenstein and W. Hu, Astrophys. J. 496, 605 (1998).
- Shaw and Mota (2008) D. J. Shaw and D. F. Mota, Astrophys. J. S. 174, 277 (2008).
- Hamilton et al. (1991) A. J. S. Hamilton, P. Kumar, E. Lu, and AlexMatthews, Astrophys. J. 374, L1 (1991).
- Navarro et al. (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996).
- Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997).