The impact of massive neutrinos on the abundance of massive clusters
Abstract
We study the spherical, tophat 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 upperbounds 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 tophat 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 freestreaming 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 lineartheory 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 .
pacs:
98.80.Es,14.60.Pq,98.65.DxI Introduction
Galaxy clusters are the most massive, gravitationally bound objects in
the universe, therefore, their abundance and its redshiftevolution are
very sensitive to cosmology including the nature of dark energy
Vikhlinin et al. (2009); Rozo et al. (2010) as well as the primordial
nonGaussianity Dalal et al. (2008); Oguri (2009). There are various ongoing
and planned surveys capable of finding many clusters under a homogeneous
and wellcontrolled/calibrated selection; e.g., optical surveys such as
the Sloan Digital Sky Survey (SDSS) Koester et al. (2007), Subaru Hyper
SuprimeCam Survey
Miyazaki et al. (2006)
The neutrino oscillation experiments have revealed that the standard threeflavor neutrinos have nonzero masses, implying that the BigBang relic neutrinos contribute to the presentday 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, largescale 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 Nbody 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 massscale neutrinos have a fast freestreaming motion, larger than the gravityinduced bulk motion, at relevant redshift and it is difficult to represent the (perturbed) FermiDirac distribution with a finite number of Nbody particles at every spatial position (see Brandbyge and Hannestad, 2009, for an attempt on the gridbased 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 tophat spherical collapse model for a MDM model fully taking into account the effects of multicomponent 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 tophat 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 radiationdominated regime to the present time or from the sufficiently linear regime to the nonlinear regime, where the initial density perturbations are on superhorizon scales Naoz and Barkana (2007). Hence, we properly take into account the transition of perturbations from the super to subhorizon scales. Thirdly, to study the neutrino perturbations around the tophat 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), timevarying 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 lineartheory extrapolated critical density as indeed motivated by the tophat 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 7year result Komatsu et al. (2010); the presentday density parameters of matter and baryon are and , respectively; the dimensionless 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.
Ii Methodology
In this section, we develop a method for solving a spherical tophat collapse of CDM overdensity in a multicomponent system, which consists of radiation (R), baryon (b), cold dark matter (c) and massive neutrinos ().
ii.1 Evolutionary equation of tophat overdensity on superhorizon scale
To achieve a sufficient accuracy of the spherical tophat 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 clusterscale 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 superhorizonscale perturbations, we employ the following approach.
Since we are interested in a spherical tophat overdensity region, time evolution of such a tophat region can be most readily described by a perturbed Friedman universe (e.g., Padmanabhan (1993)). Here we meant by “perturbed” that the tophat region is described by a Friedman universe with a small positive curvature that corresponds to the tophat density perturbation, where the tophat region is embedded in the background Friedman universe that has a flat geometry. Thus the equation of motion for the boundary radius of tophat overdensity region is given as
(1) 
where is the radius of the tophat 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 nonrelativistic, respectively). Also note that dark energy contribution is negligible in a radiation dominated regime. The scale factor for the background universe obeys
(2) 
We use the linear theory predictions to determine the initial conditions of perturbations. The spherical tophat collapse model is equivalent to the case that the perturbations are solved under the synchronous gauge condition. In this case, the superhorizonscale 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 :
(3) 
where and . By inserting Eq. (3) into Eq. (1) at the initial time we can express the effective curvature parameter in terms of the initial overdensity as
(4) 
where we have used the adiabatic initial condition to reexpress radiation perturbation in terms of matter perturbation as . Hence Eq. (1) can be rewritten as
(5) 
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 tophat overdensity until the tophat 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 :
(6) 
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 subhorizon 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 perturbation
CDM plays a major role in the spherical collapse model. When the tophat 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:
(7) 
where is the radius of the tophat 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 tophat overdensity region discussed in the preceding section. The quantity is the mass fluctuation within the CDMoverdensity 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 tophat region holds:
(8) 
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 collapseredshift 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 tophat region, because baryon is dragged out of the CDM potential
well due to the
tight
coupling with photons before the decoupling epoch,
and neutrinos are freestreaming 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:
(9) 
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).
Baryon perturbation
Next let us consider evolution of baryon perturbation in a CDM tophat overdensity region.
After the horizoncrossing of the spherical tophat overdensity region until the decoupling epoch (), baryon is tightly coupled to photon and the baryon density perturbation cannot grow. More exactly, the baryonphoton 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 tophat overdensity profile as in CDM, the baryon perturbation becomes increasingly spatiallyextended than the CDM tophat 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 tophat region is during epochs from the horizoncrossing 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 tophat region had the horizoncrossing.
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 tophat region just like the CDM case (Eq. [7]):
(10) 
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 baryonphoton 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.
Neutrino perturbation
The neutrino perturbation on subhorizon scales cannot be captured by the CDM potential well due to the large freestreaming velocity. As a result, the neutrino perturbation around the CDM overdensity extends out to radius comparable with the freestreaming scale that is given as in the units of comoving scale; here is the velocity dispersion of the FermiDirac 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 FermiDirac 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):
(11)  
where and are the metric perturbations in the Newtonian gauge, and the perturbed distribution function is expanded in terms of the Legendre polynomials as
(12) 
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
(13) 
where denotes the Fouriertransformed coefficients of the th component (c, b, ). We insert the nonlinear tophat 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 tophat region as the coordinate center, which makes the potential dependent only on the length of wavevector ; . The corresponding potential for the lineartheory extrapolated density perturbations is computed from . Since we have assumed that CDM and baryon overdensities take spherical tophat profiles, whose radii are and , respectively, the Fouriertransformed counterparts of CDM and baryon perturbations are given analytically as
(14) 
where c or b, and are the mean overdensities of CDM and baryon within their respective tophat regions. We can compute the neutrino perturbation in Fourier space from the zeroth moment of the perturbed distribution function:
(15) 
The corresponding linear perturbation is given as , a standard output of the CAMB code.
We also need to compute the realspace 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 tophat regions are
(16)  
(17) 
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 tophat CDM overdensity
region
(18) 
where is the Fourier transform of the CDM tophat overdensity (Eq. [14]), 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 zeroth 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 tophat overdensity, where is the linear transfer of the zeroth moment of the perturbed distribution function. Similarly, we compute the higherorder 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 tophat 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 tophat collapse model in a multicomponent system of CDM, baryon and neutrinos.

Choose a target mass scale of halo, , to determine the comoving scale of the spherical tophat 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 superhorizon scales in early epochs, enter into the horizon, and evolve on subhorizon scales.

Choose a target collapse redshift that corresponds to the halo formation. Then, as a first guess, normalize the initial tophat CDM overdensity, , for a given cosmological model in such a way that the lineartheory extrapolated overdensity satisfies the condition , a prediction for the collapse redshift for Einstein deSitter 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 lineartheory 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 tophat 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 tophat 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 tophat 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 tophat profile. However, a spherical tophat 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 tophat overdensity region can be interpreted as the average density contrast around such density peaks after making a tophat 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, tophat 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 Results
iii.1 Spherical collapse in a mixed dark matter model
To compute the spherical collapse model, we assume a flatgeometry 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 threespecies 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 freestreaming 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 tophat 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 tophat 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 tophat region collapses at redshift for a model without massive neutrino.
The plot also shows the baryon density contrast within the CDM tophat region as well as the radial profile of neutrino perturbations out to radii outside the tophat region. Note that, for the baryon perturbation, we show only its density contrast within the CDM tophat region (more exactly speaking, the baryon overdensity region is slightly more spatiallyextended 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 tophat 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 freestreaming velocity. To be more precise, the presentday freestreaming 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 tophat region, and become to have the radial profile. The neutrino perturbation peaks at the center of the CDM tophat 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 tophat 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 tophat region for each component. Note that we computed the neutrino density contrast by averaging the density profile within the tophat 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 deSitter model. At later collapse redshift, the cosmological constant becomes dominant in the cosmic expansion, and the collapse redshift differs from the Einstein deSitter 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 Nbody 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 lineartheory 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 deSitter 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 linearlyextrapolated overdensity for CDM perturbation alone, , is indeed greater than 1.686. That is, the initial CDM tophat 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 selfgravity 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 deSitter 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 presentday dark matter density fixed. This nontrivial 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 nonrelativistic 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 nontrivial 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 lineartheory 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:
(19) 
where the function is the fitting formula that is obtained based on a suit of Nbody 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 lineartheory 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
(20) 
where is the linear power spectrum of CDM plus baryon perturbation at target redshift , and the Fouriertransformed tophat filter: . The filtering scale and halo mass are related via ( is the presentday 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 Nbody simulations for a range of CDM models varying around the fiducial cosmological model consistent with the WMAP data:
(21) 
where , , and . For the peak height , Bhattacharya et al. (2011) simply used the fixed critical overdensity , the value of Einstein deSitter model, and then found the bestfit 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 deSitter 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 bestfit 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 Nbody 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 scalehalos 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 Xray 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 CMBconstrained 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 .
Iv Summary
In this paper, we have developed a method to solve the nonlinear dynamics of tophat 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 tophat CDM perturbation, more precisely and the velocity of tophat radius . For example, we cannot employ the lineartheory prediction for an Einstein deSitter 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 phasespace 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 lineartheory 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 lineartheory 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 lowerbound 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 cuttingedge 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 widearea 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 multicomponent 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 secondorder Lagrangian perturbation theory allows one to set up more accurate initial conditions of Nbody simulations that are simulations for a model with CDMalone or single cold component Scoccimarro and Sheth (2002); Jenkins (2010). We can extend this analysis to a multicomponent system; we can apply the secondorder 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 SPHtype 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, ChungPei 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 GrantinAid for the Scientific Research Fund (Nos. 21740177 and 23340061), by JSPS CoretoCore 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.
Footnotes
 Email address: ichiki@a.phys.nagoyau.ac.jp
 Email address: masahiro.takada@ipmu.jp
 http://sumire.ipmu.jp/en/
 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
References
 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:astroph/0701268.
 S. Miyazaki, Y. Komiyama, H. Nakaya, Y. Doi, H. Furusawa, P. Gillingham, Y. Kamata, K. Takeshi, and K. Nariai, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series (2006), vol. 6269 of Society of PhotoOptical 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 eprints (2011), eprint 1101.1290.
 S. Wang, Z. Haiman, W. Hu, J. Khoury, and M. May, Physical Review Letters 95, 011302 (2005), eprint arXiv:astroph/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:astroph/0512374.
 W. Hu, D. J. Eisenstein, and M. Tegmark, Physical Review Letters 80, 5255 (1998), eprint arXiv:astroph/9712057.
 K. Abazajian and S. Dodelson, Physical Review Letters 91, 041301 (2003), eprint arXiv:astroph/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 eprints (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:astroph/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 eprints (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:astroph/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 largescale 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:astroph/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:astroph/9804015.
 M. Bartelmann, M. Doran, and C. Wetterich, J. Cosmo. Astropart. Phys. 454, 27 (2006), eprint arXiv:astroph/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 eprints (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:astroph/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 eprints (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