The impact of massive neutrinos on the abundance of massive clusters
We study the spherical, top-hat collapse model for a mixed dark matter model including cold dark matter (CDM) and massive neutrinos of mass scales ranging from to a few eV, the range of lower- and upper-bounds implied from the neutrino oscillation experiments and the cosmological constraints. To develop this model, we properly take into account relative differences between the density perturbation amplitudes of different components (radiation, baryon, CDM and neutrinos) around the top-hat CDM overdensity region assuming the adiabatic initial conditions. Furthermore, we solve the linearized Boltzmann hierarchy equations to obtain time evolution of the lineariezed neutrino perturbations, yet including the effect of nonlinear gravitational potential due to the nonlinear CDM and baryon overdensities in the late stage. We find that the presence of massive neutrinos slows down the collapse of CDM (plus baryon) overdensity, however, that the neutrinos cannot fully catch up with the the nonlinear CDM perturbation due to its large free-streaming velocity for the ranges of neutrino masses and halo masses we consider. We find that, just like CDM models, the collapse time of CDM overdensity is well monitored by the linear-theory extrapolated overdensity of CDM plus baryon perturbation, smoothed with a given halo mass scale, if taking into account the suppression effect of the massive neutrinos on the linear growth rate. Using these findings, we argue that the presence of massive neutrinos of mass scales or 0.1 eV may cause a significant decrease in the abundance of massive halos compared to the model without the massive neutrinos; e.g., by 25% or factor 2, respectively, for halos with and at .
Galaxy clusters are the most massive, gravitationally bound objects in
the universe, therefore, their abundance and its redshift-evolution are
very sensitive to cosmology including the nature of dark energy
Vikhlinin et al. (2009); Rozo et al. (2010) as well as the primordial
non-Gaussianity Dalal et al. (2008); Oguri (2009). There are various ongoing
and planned surveys capable of finding many clusters under a homogeneous
and well-controlled/calibrated selection; e.g., optical surveys such as
the Sloan Digital Sky Survey (SDSS) Koester et al. (2007), Subaru Hyper
Miyazaki et al. (2006)
The neutrino oscillation experiments have revealed that the standard three-flavor neutrinos have non-zero masses, implying that the Big-Bang relic neutrinos contribute to the present-day mean matter density by at least about 0.4 and 0.8% if the neutrinos follow the normal mass and inverted mass hierarchies corresponding to the lower bounds on the sum of neutrino masses, and eV, respectively (e.g., Takada et al., 2006). While the absolute neutrino mass scale is still unknown, large-scale structure probes provide a powerful means of constraining the neutrino mass Hu et al. (1998); Abazajian and Dodelson (2003); Wang et al. (2005). In fact the cosmological probes have put currently the most stringent upper bound on the sum of neutrino masses; – 0.8 eV (95% C.L.), the different bounds for different probes that employ the different level of assumptions on nonlinear structure formation Ichiki et al. (2009); Vikhlinin et al. (2009); Saito et al. (2011). The ongoing and upcoming cosmological surveys promise to further tighten the neutrino mass constraint and have the potential to detect the absolute mass scale, rather than the upper bound, if systematic errors are well under control. See Abazajian et al. (2011) for the current status of the cosmological neutrino constraints and the future prospect.
However, the previous cluster cosmology experiments rest on the use of halo mass function calibrated based on N-body simulations that ignore the effect of massive neutrinos (e.g., Tinker et al. (2008); Bhattacharya et al. (2011)). Although there are several attempts to simulate nonlinear structure formation in a mixed dark matter (CDM plus massive neutrinos) model (see Klypin et al. (1993, 1997) for the pioneer work and Brandbyge et al. (2008); Brandbyge and Hannestad (2009); Brandbyge et al. (2010); Viel et al. (2010) for the recently revisited attempts), it is still very difficult to accurately simulate the structure formation especially for the neutrino mass scales of a few 0.1 eV or lighter, because such light mass-scale neutrinos have a fast free-streaming motion, larger than the gravity-induced bulk motion, at relevant redshift and it is difficult to represent the (perturbed) Fermi-Dirac distribution with a finite number of N-body particles at every spatial position (see Brandbyge and Hannestad, 2009, for an attempt on the grid-based simulation of neutrinos). There are also several attempts Saito et al. (2008); Wong (2008); Saito et al. (2009); Lesgourgues et al. (2009) aimed at developing the perturbation theory based approach to analytically model the nonlinear structure formation for a MDM model by extending the linear perturbation theory Ma and Bertschinger (1995). However, the perturbation theory based model is only valid up to the quasi nonlinear regime and breaks down for the nonlinear regime relevant for halo formation.
Therefore, the purpose of this paper is to develop a top-hat spherical collapse model for a MDM model fully taking into account the effects of multi-component system (radiation, baryon, CDM and neutrinos) (Tomita, 1969; Gunn and Gott, 1972, for the pioneer work on the spherical collapse model for a CDM model) (also see Peebles, 1980; Padmanabhan, 1993). There are several key ingredients to include in developing this model. Firstly, we carefully account for differences in the density perturbation amplitudes of each component around the top-hat CDM overdensity region assuming the adiabatic initial conditions as predicted from standard inflation scenario. Secondly, to achieve the desired accuracy, we solve time evolution of the density perturbations from the deeply radiation-dominated regime to the present time or from the sufficiently linear regime to the non-linear regime, where the initial density perturbations are on super-horizon scales Naoz and Barkana (2007). Hence, we properly take into account the transition of perturbations from the super- to sub-horizon scales. Thirdly, to study the neutrino perturbations around the top-hat CDM overdensity, which cannot be treated as a fluid, we properly solve the linearized Boltzmann hierarchy equations Ma and Bertschinger (1995) taking into account the nonlinear gravitational potential well due to the nonlinear CDM and baryon perturbations.
The spherical collapse model is an approximated method for studying the nonlinear dynamics due to the unrealistic symmetry assumed. Nevertheless, this gives a very useful tool for studying various effects on nonlinear structure formation; the spatial curvature or the cosmological constant Lahav et al. (1991); Lilje (1992), time-varying and/or clustered dark energy Wang and Steinhardt (1998); Bartelmann et al. (2006); Abramo et al. (2007); Francis et al. (2009); Eggers Bjaelde and Wong (2010); Creminelli et al. (2010), the modified dark matter scenario Oguri et al. (2003), the modified gravity scenario Schäfer and Koyama (2008); Schmidt et al. (2010); Borisov et al. (2011), and the effect of baryon perturbation Naoz and Barkana (2007). As for the effect of massive neutrinos, we know that the effect is small given the current upper bounds on the sum of neutrino masses, a few 0.1 eV (at most a few % contribution to the matter density). Hence we expect that, once the nonlinear dynamics is realized for a MDM model, we can perturbatively include the effect of massive neutrino on the CDM simulation based predictions by slightly modifying the model ingredients. For example, the halo mass function is given as a function of the peak height , where is the rms linear mass fluctuations of halo mass scale at redshift and is the linear-theory extrapolated critical density as indeed motivated by the top-hat spherical collapse model Press and Schechter (1974). Given these facts, we may be able to infer the effect of massive neutrinos on the halo mass function once the effects on and are realized. Thus, along this approach, we will also discuss the impact of massive neutrinos on the abundance of massive halos.
Throughout this paper, we employ, as our fiducial cosmological model, a flat -dominated CDM model that is consistent with the WMAP 7-year result Komatsu et al. (2010); the present-day density parameters of matter and baryon are and , respectively; the dimension-less Hubble parameter is ; the spectral tile and normalization parameter of the primordial power spectrum are , and , respectively. In most parts of our paper, when adding massive neutrinos to the fiducial cosmological model, we vary the CDM density parameter by fixing the total dark matter density to , where the energy density parameter of massive neutrino is specified by the sum of neutrino masses as Takada et al. (2006). We assume the three neutrino species, and assume one species among them is massive for simplicity.
In this section, we develop a method for solving a spherical top-hat collapse of CDM overdensity in a multi-component system, which consists of radiation (R), baryon (b), cold dark matter (c) and massive neutrinos ().
ii.1 Evolutionary equation of top-hat overdensity on superhorizon scale
To achieve a sufficient accuracy of the spherical top-hat collapse, we set up the initial conditions of perturbations in a deeply radiation dominated era, in our case. The primary reason of this high initial redshift is to reproduce the results in the earlier work Naoz and Barkana (2007) at the limit of massless neutrino case (), where Naoz and Barkana (2007) studied the effect of baryon perturbation (without massive neutrinos) on halos relevant for first stars at redshift . Note that, for cluster-scale halo formation, we can start from the later initial redshift such as -, but here we use the initial redshift in order to retain a broader coverage of our model validity. In such a radiation dominated era, the perturbations of interest are on superhorizon scales. To avoid gauge ambiguities that may arise in dealing with such superhorizon-scale perturbations, we employ the following approach.
Since we are interested in a spherical top-hat overdensity region, time evolution of such a top-hat region can be most readily described by a perturbed Friedman universe (e.g., Padmanabhan (1993)). Here we meant by “perturbed” that the top-hat region is described by a Friedman universe with a small positive curvature that corresponds to the top-hat density perturbation, where the top-hat region is embedded in the background Friedman universe that has a flat geometry. Thus the equation of motion for the boundary radius of top-hat overdensity region is given as
where is the radius of the top-hat overdensity region of matter (CDM and baryon) and radiation, and are the mean energy densities of matter (CDM plus baryon) and radiation, and and are their overdensities (e.g., defined as ), respectively. The dot notation denotes the time derivative, and the constant is the effective curvature parameter which is given in terms of the initial overdensity (see below). In this regime, massive neutrinos with mass scales we are interested in are still relativistic and contribute the radiation. Note because of on superhorizon scales for adiabatic initial conditions, where is the mean density of total matter ( or when the massive neutrino is relativistic or non-relativistic, respectively). Also note that dark energy contribution is negligible in a radiation dominated regime. The scale factor for the background universe obeys
We use the linear theory predictions to determine the initial conditions of perturbations. The spherical top-hat collapse model is equivalent to the case that the perturbations are solved under the synchronous gauge condition. In this case, the superhorizon-scale perturbations grow as Naoz and Barkana (2007); Ma and Bertschinger (1995). Assuming this growing mode and using the mass conservation yield the initial condition for the velocity :
where we have used the adiabatic initial condition to re-express radiation perturbation in terms of matter perturbation as . Hence Eq. (1) can be re-written as
The relation between and follows from the mass conservation; . Hence, given the initial conditions on (or ) and , Eq. (5) can be solved numerically to obtain time evolution of the top-hat overdensity until the top-hat region enters into the horizon.
Eq. (5) is an exact equation that can be applied even if the perturbation amplitude is large (unrealistic though). In the radiation dominated regime, the linear theory gives a good approximation. By linearizing Eq. (5) we can derive the differential equation which governs time evolution of the linear perturbation :
where we have used the fact at the initial redshift. Again note to a good approximation in this regime.
ii.2 Evolutionary equations of perturbations on sub-horizon scales
When the overdensity region enters into the horizon, perturbations of different components evolve in different ways; the CDM overdensity continues to grow, and baryon and neutrinos cannot grow together with CDM. We below describe our treatments of each component’s dynamics on subhorizon scales after the horizon crossing.
CDM plays a major role in the spherical collapse model. When the top-hat overdensity region enters into the horizon, we use the following equation, obtained in the Newtonian gauge, in order to solve the dynamics up to the nonlinear collapse of spherical CDM overdensity region:
where is the radius of the top-hat region of CDM overdensity, and and denote the mean energy and pressure densities, which determine the cosmic expansion history over the range of radiation, matter and dark energy dominated eras. Note at the horizon crossing, where is the radius of the initial top-hat overdensity region discussed in the preceding section. The quantity is the mass fluctuation within the CDM-overdensity sphere, and includes contributions from CDM, baryon and neutrino perturbations: . Note that we ignore perturbations of dark energy in this paper (see Abramo et al., 2007; Eggers Bjaelde and Wong, 2010, for the spherical collapse model with dark energy perturbations, but without massive neutrinos).
For CDM perturbation, the mass conservation within the top-hat region holds:
where is the overdensity at time , and the quantities with subscript “” denote their quantities at the horizon crossing. For the perturbations of interest, the horizon crossing is earlier than the decoupling epoch, where the perturbations are well in the linear regime. In the equation above we have used at the initial redshift, so used the fact . We determine , and by matching the values to those from Eq. (5) at the horizon crossing. We again stress that, by matching these initial conditions from superhorizon to subhorizon scales, we can achieve a sufficiently accurate setup of the initial conditions needed for the nonlinear spherical collapse dynamics, even for a high collapse-redshift such as .
To solve Eq. (7), we need to specify the mass
fluctuations of baryon and massive neutrinos, and
, within the spherical region of radius at
each time step. However, unlike CDM, the mass fluctuations are not
conserved within the CDM top-hat region, because baryon is dragged out of the CDM potential
well due to the
coupling with photons before the decoupling epoch,
and neutrinos are free-streaming out of the CDM potential well due
to their large thermal velocities
After the decoupling epoch , baryon becomes cold, and follows the mass conservation. We will below describe our treatments of baryon and neutrino perturbations in subsequent subsections.
Linearizing Eq. (7) yields the following equation to describe the evolution of linear CDM perturbation:
where is the linear CDM overdensity, and and are the linear density perturbations averaged within the sphere of radius . Here the comoving radius of is set to the same in the linear theory: . The initial conditions for and are set by matching to Eq. (6) at the horizon crossing. Before the decoupling epoch and until sufficiently higher redshift before halo formation, the equation above gives a good approximation to Eq. (7).
Next let us consider evolution of baryon perturbation in a CDM top-hat overdensity region.
After the horizon-crossing of the spherical top-hat overdensity region until the decoupling epoch (), baryon is tightly coupled to photon and the baryon density perturbation cannot grow. More exactly, the baryon-photon fluid oscillates according to the acoustic sound wave in the CDM potential well – the baryon acoustic oscillations (BAO) Hu and Sugiyama (1995). The characteristic scale of this clustering is about 150 Mpc for our fiducial cosmological model. Therefore, even if the baryon perturbation initially had a top-hat overdensity profile as in CDM, the baryon perturbation becomes increasingly spatially-extended than the CDM top-hat region as time goes by until . The linear perturbation theory predicts that the baryon density perturbation becomes much smaller in the amplitude than the CDM density perturbation at the decoupling epoch for length scales of interest; at for scales of interest.
Therefore we simply assume that the baryon density perturbation averaged within the CDM top-hat region is during epochs from the horizon-crossing to the decoupling epoch; , where is specified once a background cosmological model is given, e.g., from the CAMB code Lewis et al. (2000). We have checked that the spherical collapse of CDM overdensity is not changed even if we instead use the different assumption; , where is the baryon overdensity when the CDM top-hat region had the horizon-crossing.
After the decoupling epoch baryon becomes a cold component, and can cluster together with CDM perturbation. We can thus use the spherical collapse equation to solve nonlinear evolution of baryon perturbation in the CDM top-hat region just like the CDM case (Eq. ):
Here is the radius of baryon overdensity region, which is chosen from the radius satisfying at . After the decoupling epoch, the mass conservation holds: . The initial conditions are set to at , which corresponds to , motivated by the fact that baryon-photon coupling prevents baryon perturbation from growing before the decoupling epoch. Thus the baryon velocity perturbation is different from that of CDM perturbation, given as , so the time evolution of baryon radius differs from the CDM radius . As times goes by, the baryon perturbation eventually catches up with the CDM perturbation as we will explicitly show below.
The neutrino perturbation on sub-horizon scales cannot be captured by the CDM potential well due to the large free-streaming velocity. As a result, the neutrino perturbation around the CDM overdensity extends out to radius comparable with the free-streaming scale that is given as in the units of comoving scale; here is the velocity dispersion of the Fermi-Dirac distributed neutrinos (see Appendix A in Takada et al. (2006) for the definition). For example, for neutrinos of mass scale eV, Mpc at . To model these physical processes we use the modified CAMB code to solve time evolution of neutrino clustering, where we properly take into account the nonlinear gravitational potential due to nonlinear CDM and baryon perturbations in the late stage (see Saito et al. (2008) for the similar approach to solving the evolution of mildly nonlinear perturbations).
To be more precise, we solve the linearized Boltzmann equation that governs time evolution of the neutrino distribution function . Here is the conformal time, and denote the comoving momentum and its direction, is the proper energy times scale factor , is the background distribution function (the Fermi-Dirac distribution) and is the perturbed distribution function. The Boltzmann equation in an expanding universe can be reduced to the following hierarchical equations in Fourier space Ma and Bertschinger (1995):
where and are the metric perturbations in the Newtonian gauge, and the perturbed distribution function is expanded in terms of the Legendre polynomials as
Our main interest is to study the impact of massive neutrinos on the spherical collapse of CDM overdensity in the matter or dark energy dominated era. At redshifts after the decoupling epoch, the difference between metric perturbations and , which arises from anisotropic stress, is negligible: . In this case, the potential perturbation is given by the Poisson equation as
where denotes the Fourier-transformed coefficients of the -th component (c, b, ). We insert the nonlinear top-hat overdensities of CDM and baryon into the Poisson equation above to compute the potential including the nonlinear contribution. Note that we take the center of the top-hat region as the coordinate center, which makes the potential dependent only on the length of wavevector ; . The corresponding potential for the linear-theory extrapolated density perturbations is computed from . Since we have assumed that CDM and baryon overdensities take spherical top-hat profiles, whose radii are and , respectively, the Fourier-transformed counterparts of CDM and baryon perturbations are given analytically as
where c or b, and are the mean overdensities of CDM and baryon within their respective top-hat regions. We can compute the neutrino perturbation in Fourier space from the zero-th moment of the perturbed distribution function:
The corresponding linear perturbation is given as , a standard output of the CAMB code.
We also need to compute the real-space overdensity of neutrino perturbations at each time step, which is needed for the spherical collapse model of CDM and baryon overdensities. The radial profile of neutrino density perturbation and its mass fluctuation within the CDM or baryon top-hat regions are
We use the publicly available code FFTLog Hamilton (2000) to compute the density profile at each time step, because the code allows for a fast computation of the integration involving the Bessel function kernel. Similarly, the density profile and the mass fluctuation for the linear neutrino perturbation can be obtained by using the linear density perturbation instead of .
We employ the decoupling epoch to set up the initial
conditions of neutrino perturbations as in the baryon perturbations.
Since the neutrino perturbations at are well in the
linear regime, we can determine the initial conditions by matching to
the linear perturbation theory predictions. To be more precise, assuming
the adiabatic initial conditions, we can determine the neutrino density
perturbation at around the top-hat CDM overdensity
where is the Fourier transform of the CDM top-hat overdensity (Eq. ), and the functions and are the transfer functions of massive neutrinos and CDM at , respectively. We use the CAMB outputs to obtain the transfer functions. The ratio takes into account the relative amplitude difference of neutrino and CDM perturbations at under the adiabatic initial conditions. We compute the zero-th moment of the perturbed distribution at the initial time by multiplying the CAMB output with so that Eq. (15) gives the neutrino density perturbation around the CDM top-hat overdensity, where is the linear transfer of the zero-th moment of the perturbed distribution function. Similarly, we compute the higher-order function () from the CAMB outputs multiplied by . We can obtain the subsequent evolution of neutrino perturbations by solving the Boltzmann equation hierarchies Eq.(11) given these initial conditions at .
Our approach above is still an approximation; we used the linearized Boltzmann equations where each term in the hierarchies depends linearly on perturbation quantities. In other words, even though we include the effect of nonlinear gravitational potential, we ignore nonlinear terms in the Boltzmann equations, e.g., the term proportional to , which can be important if the neutrino perturbation itself becomes nonlinear. We will come back to this issue later.
ii.3 Summary: our recipe for solving the spherical top-hat collapse model
|superhorizon||subhorizon before the decoupling ()||subhorizon after|
|CDM||perturbed FRW (Eq. 5)||sph. collapse model (Eq. 7)||sph. collapse model (Eq. 7)|
|baryon||perturbed FRW (Eq. 5)||sph. collapse model (Eq. 10)|
|massive-||perturbed FRW (Eq.5)||linearized Boltzmann eqs. (Sec. II.2.3)|
|radiation||perturbed FRW (Eq.5)||–||–|
Here is a quick summary of the procedures we take for solving the spherical top-hat collapse model in a multi-component system of CDM, baryon and neutrinos.
Choose a target mass scale of halo, , to determine the comoving scale of the spherical top-hat CDM overdensity region via the relation .
Solve the linear CDM perturbation of the comoving scale from the initial time to the present time based on the linear perturbation theory, assuming the adiabatic initial conditions. In these calculations we properly take into account the fact that the density perturbations are on super-horizon scales in early epochs, enter into the horizon, and evolve on sub-horizon scales.
Choose a target collapse redshift that corresponds to the halo formation. Then, as a first guess, normalize the initial top-hat CDM overdensity, , for a given cosmological model in such a way that the linear-theory extrapolated overdensity satisfies the condition , a prediction for the collapse redshift for Einstein de-Sitter model, a CDM model without baryon and massive neutrino contributions.
Solve the nonlinear evolution of CDM overdensity until . Iteratively solve the spherical collapse model by changing the initial overdensity amplitude until the CDM overdensity becomes to collapse at the target collapse redshift, . Also obtain the linear-theory extrapolated overdensity for CDM or CDM plus baryon perturbations, or .
Table 1 also gives a quick summary of these treatments, clarifying which equations we use for solving the spherical top-hat collapse model.
It would also be useful to explicitly list the assumptions we employ for the spherical collapse model:
We used the linearized Boltzmann hierarchy equations to solve the time evolution of linearized neutrino perturbations. However, we include the effect of nonlinear gravitational potential due to the nonlinear CDM and baryon density perturbations.
We assumed the top-hat profiles of CDM and baryon perturbations.
We set the baryon and neutrino density perturbations to zero, i.e. , before the decoupling epoch, because we found it gives a good approximation compared to a more rigorous calculation under the adiabatic initial conditions.
For the second assumption above, rigorously speaking, even if we consider a spherically symmetric top-hat overdensity region at the initial time, the radial profile become changed by the presence of radiation, baryon and neutrino perturbations, which go out of the CDM overdensity region after the horizon crossing. As a result the overdensity region no longer obeys a top-hat profile. However, a spherical top-hat collapse model is anyway an approximated method for studying the nonlinear dynamics of the initial overdensity regions, preferentially representing density peaks in the primordial perturbations. Hence the top-hat overdensity region can be interpreted as the average density contrast around such density peaks after making a top-hat filtering with the smoothing scale of target halo mass scale (e.g., Padmanabhan, 1993, for a similar discussion). For these reasons, our treatment of assuming a spherically symmetric, top-hat overdensity for CDM perturbation is adequate enough for our purpose. Our main goal is to study the impact of massive neutrinos on the spherical collapse by comparing the results with and without massive neutrino contribution. Also note that our method includes the limit of spherical collapse for a pure CDM model without baryon and massive neutrino contributions, by imposing .
iii.1 Spherical collapse in a mixed dark matter model
To compute the spherical collapse model, we assume a flat-geometry cold dark matter, dominated cosmological model (CDM) that is consistent with the WMAP results Komatsu et al. (2010). We further need to specify neutrino parameters. In this paper, we assume standard three flavors of neutrinos. Since structure formation is sensitive only to the sum of the three-species neutrino masses, we assume, for simplicity, that only one species of neutrinos are massive and other two species are massless. In this case, the neutrino free-streaming scale is shortest for a fixed total neutrino mass (or a fixed ), and therefore the neutrino has the largest ability to cluster on small scales compared to a case that the total neutrino mass is split into different species. We then study how the spherical collapse is affected by massive neutrino assuming the mass scale ranging from 0.05 to a few 0.1 eV. This range of mass scales covers the lower and upper bounds on the neutrino mass implied from the neutrino oscillation experiments and the cosmological constraints (Ichiki et al., 2009; Saito et al., 2011; Abazajian et al., 2011). In most parts of this paper, when we add the massive neutrinos, we keep the energy density of “dark matter” () and other parameters fixed, and vary .
Fig. 1 shows how the CDM top-hat overdensity grows as a function of cosmic time. We assumed for neutrino mass, a mass scale close to the lower bound implied from the normal mass hierarchy, and halo mass scale corresponding to the comoving radius, Mpc, for our fiducial CDM model. By considering massive halos, we can estimate the largest effect of neutrinos on the spherical top-hat collapse, as such a massive halo has the ability to trap neutrinos around it due to the deepest gravitational potential well. We, as a working example, set the initial CDM overdensity so that the top-hat region collapses at redshift for a model without massive neutrino.
The plot also shows the baryon density contrast within the CDM top-hat region as well as the radial profile of neutrino perturbations out to radii outside the top-hat region. Note that, for the baryon perturbation, we show only its density contrast within the CDM top-hat region (more exactly speaking, the baryon overdensity region is slightly more spatially-extended than the CDM as we described above). At sufficiently high redshifts such as , the baryon and neutrino perturbations are smaller in the amplitudes than the CDM density contrast. Then at lower redshifts, the baryon perturbation eventually catches up with the CDM perturbation. For this particular case, at lower redshifts , the comoving radius of CDM top-hat region starts to shrink, entering into the nonlinear regime (or equivalently deviating from the linear evolution). The figure explicitly demonstrates that the CDM and baryon perturbations, i.e. cold components, can collapse together having at the collapse redshift.
On the other hand, the neutrinos of this mass scale cannot catch up with the CDM perturbation due to the large free-streaming velocity. To be more precise, the present-day free-streaming scale in the comoving scale unit is (see Appendix A in Takada et al. (2006)), which is much larger than a few Mpc, a scale of the virial radius of massive halos. The plot shows that the neutrinos are indeed clustering around the CDM top-hat region, and become to have the radial profile. The neutrino perturbation peaks at the center of the CDM top-hat region, but the density contrast is still smaller than unity, so not yet in the highly nonlinear regime. More precisely, the neutrino perturbation averaged within the CDM top-hat region at the collapse redshift. We have also checked that, when neutrino mass is in the range smaller than a few 0.1 eV, the neutrino density contrast grows only up to the weakly nonlinear regime a few even at the collapse time (for the case of eV, at the collapse redshift). Note that, for the halo of and the collapse redshift , and for and eV, respectively). Therefore our approach using the linearized Boltzmann equations for neutrino perturbations is approximately validated. Once the halo is formed via virialization of the kinetic and gravitational bound energies, the neutrino would become stably clustered around the halo region as studied in Singh and Ma (2003); Abazajian et al. (2005). Nevertheless the resulting neutrino overdensity is much smaller than the CDM and baryon perturbations in a halo region, therefore we ignore the neutrino mass contribution to the halo mass in the following analysis for simplicity.
Fig. 2 shows the time evolution of density contrasts within the top-hat region for each component. Note that we computed the neutrino density contrast by averaging the density profile within the top-hat region. The CDM and baryon collapse at having . The neutrinos are affected by the nonlinear clustering of CDM perturbations, but do not enter into the highly nonlinear regime.
In Fig. 3, we compare the spherical collapses of CDM perturbation for models with and without baryon perturbation and/or massive neutrino contributions. We used the same initial conditions of CDM perturbation except for the result without baryon perturbation (labelled as “w/o and ”). For the model without baryon perturbation (thin solid curve), we set the CDM perturbation amplitude to match the CDM amplitude at for our fiducial model (without massive neutrino), but set . In this case, if we set a sufficiently early collapse redshift, this model gives the collapse redshift given by , the case for an Einstein de-Sitter model. At later collapse redshift, the cosmological constant becomes dominant in the cosmic expansion, and the collapse redshift differs from the Einstein de-Sitter prediction. This result is compared with other curves. First, the dashed curve shows the collapse of CDM perturbation when the baryon perturbation is added. The presence of baryon perturbation, which cannot grow during epochs from the horizon enter until the decoupling epoch, delays the spherical collapse. The delay is rather significant, from to , and therefore the result suggests that we need to carefully take into account the effect of baryon perturbation on the nonlinear structure formation in N-body simulations, especially when we set up the initial conditions (we will discuss this issue later). The bold solid and dotted curves show the spherical collapse when massive neutrino is further added for mass scales of and eV, respectively. The collapse redshifts are further delayed as and 0.45, respectively. These mass scales correspond to the lower bounds on total neutrino mass for the normal and inverted mass hierarchies. Thus, adding the smoother, massive components into the CDM perturbation progressively delays the spherical collapse.
We then study the linear-theory extrapolated overdensity at the collapse redshift (we will often call it the critical overdensity hereafter). In Fig. 4 we show the critical overdensity of CDM plus baryon perturbation, , smoothed with length scales and corresponding to halo mass scales and , respectively (the subscript “” stands for CDM plus baryon). The overdensity can serve as a clock to infer the collapse redshift, because it can be easily computed once the initial power spectra of CDM and baryon perturbations, halo mass scale and cosmological model are specified. In an Einstein de-Sitter universe (), which includes CDM alone, the critical overdensity can be derived analytically Tomita (1969); Gunn and Gott (1972) and is found to be , independently of halo mass and collapse redshift. Even for a model with curvature or dark energy contribution, the critical density to a good approximation, independently of halo mass Lahav et al. (1991); Lilje (1992); Wang and Steinhardt (1998). The top solid curve shows the critical overdensity when baryon is included, but massive neutrino is ignored. The critical overdensity differs from 1.686, and the change of is due to the presence of baryon perturbation, which has a smaller amplitude than the CDM perturbation at higher redshifts (e.g., see Fig. 1). Our result is consistent with the result in Naoz and Barkana (2007), where they studied the effect of baryon perturbation on the spherical collapse at high redshift for much smaller halos that are relevant for first stars. Although the presence of smoother baryon perturbation delays the spherical collapse, it leads to the smaller than 1.686. However, note that the linearly-extrapolated overdensity for CDM perturbation alone, , is indeed greater than 1.686. That is, the initial CDM top-hat overdensity with greater amplitude than expected from is needed so that it collapses at a given collapse redshift . The curve peaks around () having . At the lower redshifts than , especially at , the critical overdensity becomes smaller than the peak value due to the effect of the cosmological constant. When the cosmological constant or more generally dark energy becomes to dominate the cosmic energy density, the accelerating cosmic expansion slows down the growth of CDM plus baryon perturbation, and delays the spherical collapse. This yields the smaller critical overdensity. Both the linear and nonlinear growths of CDM plus baryon perturbation are delayed by the cosmic acceleration, and the linear growth is more suppressed than the nonlinear growth, because the spherical collapse eventually separates from the cosmic expansion in the nonlinear stage, and becomes more affected by the self-gravity of nonlinear CDM plus baryon overdensity. For these reasons, when the growth of density perturbations is suppressed by the faster cosmic expansion than the Einstein de-Sitter model, it generally leads to the smaller critical overdensity than 1.686.
The other curves in Fig. 4 show the results for when including the massive neutrino for a fixed total matter density . The presence of massive neutrino further delays the spherical collapse (see Fig. 3), and in turn leads to the smaller critical overdensity , the same trend for the effects of baryon perturbation and the cosmological constant. However, the massive neutrino only decreases by less than 0.1% compared to the solid curve, for these halo mass and neutrino mass scales and over a range of redshifts we have studied. This small change in can be contrasted with the effect on the linear growth rate; the growth rate is suppressed by the amount of at relevant redshift compared to the growth rate without the massive neutrino Hu et al. (1998); Takada et al. (2006), corresponding to 1.6 and 3.2% suppression for the neutrino mass scales of 0.05 and 0.1 eV, respectively. The results imply that the neutrino effect on the spherical collapse is well captured by the linear growth rate of CDM plus baryon perturbation.
The different curves show that the change in is not monotonic with changing neutrino masses, when keeping the present-day dark matter density fixed. This non-trivial dependence can be understood as follows. The neutrino effect on the spherical collapse arises from its effect on the cosmic expansion history and the gravitational collapse of CDM perturbation. First, the presence of massive neutrino leads to a faster cosmic expansion during the neutrino was relativistic, which slows down the growth of CDM perturbation. Note that the neutrino becomes non-relativistic when . Secondly, the neutrino perturbation does contribute to the nonlinear gravitational collapse of CDM perturbation, and therefore accelerates the spherical collapse to some extent. The net effect arises from these competing effects. We can study which of these two effects is more important as follows. Fig. 5 shows the results when we ignore the neutrino perturbation () in both the spherical collapse calculation as well as the linear growth calculation for CDM plus baryon perturbation. The figure shows that, with increasing the neutrino mass scale, the spherical collapse more delays and the critical overdensity becomes monotonically smaller. Hence, comparing Figs. 4 and 5 manifests that the neutrino perturbation does contribute to the spherical collapse, which differs from the effect of smooth dark energy model Wang and Steinhardt (1998); Pace et al. (2010).
In Fig. 6, we summarize dependences of the critical overdensity on halo mass and neutrino mass, assuming the collapse redshift . It can be found that, for a fixed halo mass scale, the critical overdensity first decreases with increasing the neutrino mass from 0.05 eV, but then starts to increase at greater neutrino mass scales from some mass scale. The turnover neutrino mass scale slightly changes with halo mass scale. This non-trivial dependence arises depending on which of the two competing effects discussed above dominates. If we ignore the neutrino perturbation, the critical overdensity decreases with increasing neutrino mass independently of halo mass.
iii.2 The impact of massive neutrinos on halo mass function
In this subsection, we estimate the impact of massive neutrinos on the halo mass function that is one of the most important observables for cluster surveys.
As we have shown, the effect of the massive neutrino on the nonlinear gravitational collapse of CDM plus baryon perturbation is well captured by the linear growth rate. In other words, the nonlinear neutrino clustering around the CDM overdensity does not largely change the nonlinear dynamics, and therefore is very unlikely to change structural properties of mass distribution within a halo. We here assume that the halo mass function for a MDM model can be obtained from a mapping of the mass function in CDM models without massive neutrino. That is, we assume that the mapping of halo mass function can be obtained by assuming that (1) only the cold component (CDM plus baryon) can collapse to form halos, and (2) the halo mass function for a MDM model can be obtained just by replacing the linear-theory mass fluctuations appearing in the mass function for a CDM model with the corresponding mass fluctuation of CDM plus baryon perturbation for a MDM model:
where the function is the fitting formula that is obtained based on a suit of N-body simulations for CDM models. The previous works have shown that the fitting formula is well characterized in terms of the peak height, Sheth and Tormen (1999); Tinker et al. (2008), where is the linear-theory extrapolated critical overdensity for halo formation at a given redshift and is the linear rms mass fluctuation smoothed with the halo mass scale and at redshift . In Eq. (19), we assumed that we can obtain the halo mass function for a MDM model simply by using the peak height for CDM plus baryon perturbation as well as by using the prefactor , the mean mass density of CDM plus baryon, because the cold component is the collapsing component to form halos. To be more precise, the rms mass fluctuation of halo mass for CDM plus baryon perturbation is defined as
where is the linear power spectrum of CDM plus baryon perturbation at target redshift , and the Fourier-transformed top-hat filter: . The filtering scale and halo mass are related via ( is the present-day mean mass density of CDM plus baryon).
As for the fitting formula, we use the formula in Bhattacharya et al. (2011) that is obtained from N-body simulations for a range of CDM models varying around the fiducial cosmological model consistent with the WMAP data:
where , , and . For the peak height , Bhattacharya et al. (2011) simply used the fixed critical overdensity , the value of Einstein de-Sitter model, and then found the best-fit parameters and so on by fitting the functional form above with the mass function measured from simulations for variant CDM models. More exactly speaking, when we have the effects of baryon perturbation and cosmic acceleration, the critical overdensity is changed from the Einstein de-Sitter value . However, the change is very small, less than a percent level (see Fig. 4), and therefore it was assumed that the change of is absorbed by tuning the fitting model parameters. If the change of is properly taken into account, the fitting will yield slightly different best-fit model parameters of and so on. Furthermore, although the presence of baryon perturbation changes the collapse of CDM perturbation (see Fig. 3), we here assume that the simulations in Bhattacharya et al. (2011) properly take into account the effect of baryon when setting up the initial conditions of N-body simulations (see below for a further discussion). To compute in Eq. (19), we use the CAMB code Lewis et al. (2000) to compute the transfer functions. The CAMB outputs include the effect of massive neutrinos or baryon perturbations on the growth of CDM perturbation.
The upper panel of Fig. 7 compares the halo mass functions at for different models with and without massive neutrino contribution. The lower panel explicitly shows the ratio of the mass functions with and without massive neutrino. Here we assumed and eV for the neutrino mass scale, which are close to the lower bounds of the normal mass hierarchy (NH) and the inverted hierarchy (IH) that are implied from the terrestrial experiments. Hence either of these results would inevitably exist in our universe. The presence of massive neutrino decreases the abundance of halos, more significantly for more massive halos that reside in the exponential tail of mass function. The decrease in the halo abundance is up to a factor 2 around . This change can be compared to the effect on the linear mass fluctuation such as ; the neutrino of these mass scales decreases only by a few percent for neutrinos of these mass scales. Again the higher sensitivity of halo mass function to neutrino mass is through the exponential tail of mass function at massive halo ends. The thin solid curve in the lower panel (although almost overlapped with the dotted curve) shows the ratio when further taking into account the change in the critical density in the mass function (Eq. 19); more explicitly we decreases the critical density by 0.03%, a maximum change implied from Fig. 6 for the case of eV. It is clear that the change in the critical density due to the massive neutrinos causes a negligible effect on the halo mass function.
Fig. 8 shows the similar results, but for higher redshifts and , respectively. The decrease in the abundance of cluster scale-halos is more significant at higher redshifts.
One may think whether or not the effect of massive neutrino on the halo mass function is mostly described by the change of , the normalization parameter of power spectrum amplitudes often used in the literature. Fig. 9 compares the halo mass functions for a MDM model with eV and for a CDM model where at is lowered so as to match the value for the MDM model; more precisely, is changed to 0.486 from 0.500. Note that both the models have the same . The CDM model with the lowered roughly reproduces the decrease in the halo abundance. However, the two curves do not exactly agree because of the difference in the linear power spectra of CDM and baryon perturbations. Nevertheless, this gives a justification of the neutrino mass constraint derived in Vikhlinin et al. (2009), where the neutrino mass constraint is obtained from the allowed range of values that are derived by comparing the observed abundance of X-ray luminous clusters with the model halo mass functions varying within CDM models without massive neutrino contribution.
iii.3 Discussion: Cosmological parameter degeneracies
When adding the massive neutrinos for different masses, we have so far kept the total dark matter density, , fixed. For a more practical perspective, the CMB information give precise constraints on the CDM and baryon densities, and , as well as the curvature parameter or equivalently the total energy density, . Massive neutrinos with small mass scales of a few eV were relativistic before the decoupling epoch, and do not affect the CMB observables. Therefore the CMB observables cannot well constrain the neutrino mass of the small mass scales, leaving degeneracies in cosmological parameters. Given these facts one might think that, when adding the massive neutrinos, we should keep these CMB-constrained parameters fixed. If we assume a flat geometry, this is equivalent to varying either the Hubble parameter or the energy density of the cosmological constant with fixing the CMB parameters above. For example, when the neutrino mass eV is added, this leads to or from the fiducial values or , respectively.
Fig. 10 shows the critical density for halos of for a MDM model with various neutrino mass scales, where we varied either or by the amount determined by the neutrino mass scale, but fixing . The results are similar to Fig. 4; the effect of massive neutrinos on the critical density is very small for the range of cosmological models. Fig. 11 shows how the MDM models alter the halo mass function compared to the case without massive neutrinos, which can be compared with our fiducial case where the massive neutrinos of eV are added by varying the CDM density parameter . The parameter change of or also leads to the smaller abundance of massive halos as in the case changing , but the decrease is slightly smaller than the case when changing .
In this paper, we have developed a method to solve the nonlinear dynamics of top-hat CDM overdensity region including the effects of baryon perturbation and massive neutrinos. In developing the spherical collapse model, we properly set up the initial conditions of each components (baryon, CDM and neutrinos), which have different amplitudes and profiles, assuming the adiabatic initial conditions (see Fig. 1). In fact we found that the nonlinear dynamics is very sensitive to detailed setup of the initial conditions of top-hat CDM perturbation, more precisely and the velocity of top-hat radius . For example, we cannot employ the linear-theory prediction for an Einstein de-Sitter model, , to set up the initial conditions, e.g., even at an epoch in the sufficiently linear regime such as the decoupling epoch , because this solution ignores that the CDM perturbation is affected by the presence of baryon and massive neutrino.
Since we cannot treat the neutrinos as a perfect fluid, we properly solved the linearized Boltzmann hierarchy equations to compute time evolution of linearized neutrino perturbations, where we include the effect of nonlinear gravitational potential due to the nonlinear CDM and baryon perturbation in the late stage. For neutrino mass scales lighter than a few eV, the range inferred from the neutrino oscillation experiments and the cosmological constraints, the neutrino perturbation stays in the quasi nonlinear regime, (see Fig. 2). This gives a justification of our treatment where we used the linearized Boltzmann equations. As for an improved modeling, one can further include the nonlinear terms such as the coupling term between the nonlinear gravitational potential and the perturbed phase-space density of neutrinos in order to solve the time evolution of neutrino perturbations in a perturbation theory manner.
By solving the spherical collapse model for cosmological models around a CDM model that is consistent with the WMAP data, we found that both the neutrino and baryon perturbations delay the collapse of CDM overdensity compared to a model with CDM alone (Fig. 3). However, interestingly we found that the collapse redshift can be well monitored by the linear-theory extrapolated overdensity of CDM (plus baryon) perturbation(s) for the ranges of neutrino masses ( a few 0.1 eV) and halo mass scales we have considered. This result is promising because the linear-theory extrapolated overdensity (the critical density) can be accurately computed using the linear perturbation theory, once cosmological model and neutrino mass are specified. In other words, we found that the massive neutrinos with the rang of mass scales lead to only a small change in the critical density by % compared to the model without massive neutrino, but with the same (Figs. 4, 5, and 6).
Given the results of the spherical collapse model, we gave Eq. (19) to estimate the halo mass function including the effect of massive neutrinos, where the effect of massive neutrinos are properly taken into account in the linear mass fluctuations of CDM and baryon perturbations at a given redshift, smoothed with a given halo mass scale; (also see Brandbyge et al., 2010, for the similar discussion). Using the equation, we found that the presence of massive neutrinos with and 0.1 eV, the lower-bound mass scales of normal and inverted mass hierarchies, respectively, may cause a significant decrease in the abundance of massive halos; more specifically, up to a factor of 2 for halos with and at (see Figs. 7 and 8). Thus our results imply that massive neutrinos, which should exist in our universe, relax to some extent a possible tension that the cutting-edge SZ experiments could not find as many massive clusters as what was originally expected Marriage et al. (2011); Williamson et al. (2011). This needs to be further studied more carefully. Since it is still challenging to accurately simulate nonlinear structure formation in a MDM model, especially for such light neutrino mass scales of a few eV (see Brandbyge and Hannestad, 2009; Brandbyge et al., 2010; Viel et al., 2010, for the attempts), the analytical model developed in this paper will give a useful tool or at least useful guidance for interpreting ongoing and upcoming wide-area surveys of massive clusters.
Our findings also propose several applications. First, as we stressed above, a careful setup of the initial conditions is very important in order to have an accurate nonlinear dynamics, for a multi-component system with CDM, baryon and neutrinos. This implies that it is very important to set up the accurate initial conditions for cosmological simulations including the effect of baryon such as smoothed particle hydrodynamical (SPH) simulations (see Yoshida et al., 2003; Tseliakhovich and Hirata, 2010; Naoz et al., 2011, for the similar discussion). Since the spherical collapse model gives an exact solution of the nonlinear dynamics, albeit an unrealistic symmetry assumed, we can explore how to set up the initial conditions by combining the spherical collapse model with the linear and/or perturbation theory predictions. For example, it was shown that using the second-order Lagrangian perturbation theory allows one to set up more accurate initial conditions of N-body simulations that are simulations for a model with CDM-alone or single cold component Scoccimarro and Sheth (2002); Jenkins (2010). We can extend this analysis to a multi-component system; we can apply the second-order Lagrangian perturbation theory to CDM and baryon perturbations separately by taking into account the different growth rates, and then can study how the improved initial conditions can reproduce the exact solution of spherical collapse for CDM and baryon perturbations starting from a given initial redshift. Such a study will give a useful guidance for exploring how to set up the initial conditions for CDM and baryon particles in a SPH-type simulation. This can be further extended to a case further including the neutrino particles. These are our future study, and will be presented elsewhere.
Secondly, several studies recently claimed that detected massive clusters at high redshifts beyond may give a tension of CDM structure formation model (Jee et al., 2009; Holz and Perlmutter, 2010, also see references therein). It is indeed interesting to explore whether or not these particular catalogs of clusters, which are found by different observations/surveys under different selection functions, can falsify the CDM predictions as explored in Mortonson et al. (2011); Williamson et al. (2011). However, the effect of massive neutrinos has been ignored in the previous studies. Again, the method developed in this paper can be used to address how the presence of the high- massive clusters may falsify a more realistic cosmological model that includes massive neutrino contribution. This study will be presented elsewhere.
Acknowledgement: We thank Scott Dodelson, Salman Habib, Wayne Hu, Chung-Pei Ma, Daisuke Nagai, Smadar Naoz, Shun Saito, Roman Scoccimarro, Tristan Smith, and David Spergel for useful discussion. This work is supported in part by the Grant-in-Aid for the Scientific Research Fund (Nos. 21740177 and 23340061), by JSPS Core-to-Core Program “International Research Network for Dark Energy”, by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and by the FIRST program “Subaru Measurements of Images and Redshifts (SuMIRe)”, CSTP, Japan.
- E-mail address: email@example.com
- E-mail address: firstname.lastname@example.org
- See http://cmb.as.arizona.edu/ẽisenste/acousticpeak/acoustic_physics.html for a pedagogical illustration of density evolutions of different components around the CDM density perturbation peak in the linear regime.
- Again also see http//cmb.as.arizona.edu/~eisenste/acousticpeak/acoustic_physics.html for the similar method for the linear perturbations
- A. Vikhlinin, A. V. Kravtsov, R. A. Burenin, H. Ebeling, W. R. Forman, A. Hornstrup, C. Jones, S. S. Murray, D. Nagai, H. Quintana, et al., Astrophys. J. 692, 1060 (2009), eprint 0812.2720.
- E. Rozo, R. H. Wechsler, E. S. Rykoff, J. T. Annis, M. R. Becker, A. E. Evrard, J. A. Frieman, S. M. Hansen, J. Hao, D. E. Johnston, et al., Astrophys. J. 708, 645 (2010), eprint 0902.3702.
- N. Dalal, O. Doré, D. Huterer, and A. Shirokov, Phys. Rev. D 77, 123514 (2008), eprint 0710.4560.
- M. Oguri, Physical Review Letters 102, 211301 (2009), eprint 0905.0920.
- B. P. Koester, T. A. McKay, J. Annis, R. H. Wechsler, A. E. Evrard, E. Rozo, L. Bleem, E. S. Sheldon, and D. Johnston, Astrophys. J. 660, 221 (2007), eprint arXiv:astro-ph/0701268.
- S. Miyazaki, Y. Komiyama, H. Nakaya, Y. Doi, H. Furusawa, P. Gillingham, Y. Kamata, K. Takeshi, and K. Nariai, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (2006), vol. 6269 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series.
- T. A. Marriage, V. Acquaviva, P. A. R. Ade, P. Aguirre, M. Amiri, J. W. Appel, L. F. Barrientos, E. S. Battistelli, J. R. Bond, B. Brown, et al., Astrophys. J. 737, 61 (2011), eprint 1010.1065.
- R. Williamson, B. A. Benson, F. W. High, K. Vanderlinde, P. A. R. Ade, K. A. Aird, K. Andersson, R. Armstrong, M. L. N. Ashby, M. Bautz, et al., ArXiv e-prints (2011), eprint 1101.1290.
- S. Wang, Z. Haiman, W. Hu, J. Khoury, and M. May, Physical Review Letters 95, 011302 (2005), eprint arXiv:astro-ph/0505390.
- M. Oguri and M. Takada, Phys. Rev. D 83, 023008 (2011), eprint 1010.0744.
- M. Takada, E. Komatsu, and T. Futamase, Phys. Rev. D 73, 083520 (2006), eprint arXiv:astro-ph/0512374.
- W. Hu, D. J. Eisenstein, and M. Tegmark, Physical Review Letters 80, 5255 (1998), eprint arXiv:astro-ph/9712057.
- K. Abazajian and S. Dodelson, Physical Review Letters 91, 041301 (2003), eprint arXiv:astro-ph/0212216.
- K. Ichiki, M. Takada, and T. Takahashi, Phys. Rev. D 79, 023520 (2009), eprint 0810.4921.
- S. Saito, M. Takada, and A. Taruya, Phys. Rev. D 83, 043529 (2011), eprint 1006.4845.
- K. N. Abazajian, E. Calabrese, A. Cooray, F. De Bernardis, S. Dodelson, A. Friedland, G. M. Fuller, S. Hannestad, B. G. Keating, E. V. Linder, et al., ArXiv e-prints (2011), eprint 1103.5083.
- J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, G. Yepes, S. Gottlöber, and D. E. Holz, Astrophys. J. 688, 709 (2008), eprint 0803.2706.
- S. Bhattacharya, K. Heitmann, M. White, Z. Lukić, C. Wagner, and S. Habib, Astrophys. J. 732, 122 (2011), eprint 1005.2239.
- A. Klypin, J. Holtzman, J. Primack, and E. Regos, Astrophys. J. 416, 1 (1993), eprint arXiv:astro-ph/9305011.
- A. Klypin, R. Nolthenius, and J. Primack, Astrophys. J. 474, 533 (1997).
- J. Brandbyge, S. Hannestad, T. Haugbølle, and B. Thomsen, J. Cosmo. Astropart. Phys. 8, 20 (2008), eprint 0802.3700.
- J. Brandbyge and S. Hannestad, J. Cosmo. Astropart. Phys. 5, 2 (2009), eprint 0812.3149.
- J. Brandbyge, S. Hannestad, T. Haugboelle, and Y. Y. Y. Wong, ArXiv e-prints (2010), eprint 1004.4105.
- M. Viel, M. G. Haehnelt, and V. Springel, J. Cosmo. Astropart. Phys. 6, 15 (2010), eprint 1003.2422.
- S. Saito, M. Takada, and A. Taruya, Physical Review Letters 100, 191301 (2008), eprint 0801.0607.
- Y. Y. Y. Wong, J. Cosmo. Astropart. Phys. 10, 35 (2008), eprint 0809.0693.
- S. Saito, M. Takada, and A. Taruya, Phys. Rev. D 80, 083528 (2009), eprint 0907.2922.
- J. Lesgourgues, S. Matarrese, M. Pietroni, and A. Riotto, J. Cosmo. Astropart. Phys. 6, 17 (2009), eprint 0901.4550.
- C. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), eprint arXiv:astro-ph/9506072.
- K. Tomita, Progress of Theoretical Physics 42, 9 (1969).
- J. E. Gunn and J. R. Gott, III, Astrophys. J. 176, 1 (1972).
- P. J. E. Peebles, The large-scale structure of the universe (1980).
- T. Padmanabhan, Structure Formation in the Universe (1993).
- S. Naoz and R. Barkana, Mon. Not. Roy. Astron. Soc. 377, 667 (2007), eprint arXiv:astro-ph/0612004.
- O. Lahav, P. B. Lilje, J. R. Primack, and M. J. Rees, Mon. Not. Roy. Astron. Soc. 251, 128 (1991).
- P. B. Lilje, Astrophys. J. Letters 386, L33 (1992).
- L. Wang and P. J. Steinhardt, Astrophys. J. 508, 483 (1998), eprint arXiv:astro-ph/9804015.
- M. Bartelmann, M. Doran, and C. Wetterich, J. Cosmo. Astropart. Phys. 454, 27 (2006), eprint arXiv:astro-ph/0507257.
- L. R. Abramo, R. C. Batista, L. Liberato, and R. Rosenfeld, J. Cosmo. Astropart. Phys. 11, 12 (2007), eprint 0707.2882.
- M. J. Francis, G. F. Lewis, and E. V. Linder, Mon. Not. Roy. Astron. Soc. 393, L31 (2009), eprint 0810.0039.
- O. Eggers Bjaelde and Y. Y. Y. Wong, ArXiv e-prints (2010), eprint 1009.0010.
- P. Creminelli, G. D’Amico, J. Noreña, L. Senatore, and F. Vernizzi, J. Cosmo. Astropart. Phys. 3, 27 (2010), eprint 0911.2701.
- M. Oguri, K. Takahashi, H. Ohno, and K. Kotake, Astrophys. J. 597, 645 (2003), eprint arXiv:astro-ph/0306020.
- B. M. Schäfer and K. Koyama, Mon. Not. Roy. Astron. Soc. 385, 411 (2008), eprint 0711.3129.
- F. Schmidt, W. Hu, and M. Lima, Phys. Rev. D 81, 063005 (2010), eprint 0911.5178.
- A. Borisov, B. Jain, and P. Zhang, ArXiv e-prints (2011), eprint 1102.4839.
- W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).
- E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G.&nb