The small-scale Universe

Power spectrum for the small-scale Universe


The first objects to arise in a cold dark matter universe present a daunting challenge for models of structure formation. In the ultra small-scale limit, CDM structures form nearly simultaneously across a wide range of scales. Hierarchical clustering no longer provides a guiding principle for theoretical analyses and the computation time required to carry out credible simulations becomes prohibitively high. To gain insight into this problem, we perform high-resolution () simulations of an Einstein-de Sitter cosmology where the initial power spectrum is with . Self-similar scaling is established for and more convincingly than in previous, lower-resolution simulations and for the first time, self-similar scaling is established for an simulation. However, finite box-size effects induce departures from self-similar scaling in our simulation. We compare our results with the predictions for the power spectrum from (one-loop) perturbation theory and demonstrate that the renormalization group approach suggested by McDonald (2007) improves perturbation theory’s ability to predict the power spectrum in the quasilinear regime. In the nonlinear regime, our power spectra differ significantly from the widely used fitting formulae of Peacock & Dodds (1996) and Smith et al. (2003) and a new fitting formula is presented. Implications of our results for the stable clustering hypothesis vs. halo model debate are discussed. Our power spectra are inconsistent with predictions of the stable clustering hypothesis in the high- limit and lend credence to the halo model. Nevertheless, the fitting formula advocated in this paper is purely empirical and not derived from a specific formulation of the halo model.

Methods: N-body simulations — cosmology: dark matter

1 Introduction

In the standard cosmological paradigm, present-day structures arise from small-amplitude density perturbations which have their origin in the very early Universe. These primordial perturbations are presumed to form a Gaussian random field whose statistical properties are described entirely by the power spectrum, . While higher order statistics are required to describe the density field once nonlinearities develop, the power spectrum remains central to our understanding of structure formation.

During the radiation-dominated phase of the standard cold dark matter (CDM) scenario, the power spectrum evolves from its primordial, approximately power-law form, to one in which its logarithmic slope, , decreases from on large scales to on small scales6. At the start of the matter-dominated phase, which signals the beginning of structure formation, the dimensionless power spectrum, , decreases monotonically with scale. The implication is that structure forms from the bottom up. Hierarchical clustering, as this process has come to be known, is the central idea in our understanding of structure formation. Hierarchical clustering also explains why cosmological N-body simulations are able to provide a reasonable facsimile of true cosmological evolution; with enough dynamic range, simulations are able to follow the development of virialized, highly nonlinear systems on small scales while properly modeling the large-scale tidal fields that shape them. However, the dynamic-range requirement becomes increasingly difficult to achieve as since, in this limit, becomes independent of and structures collapse nearly simultaneously across a wide range of scales. Put another way, as , the infrared divergence of the power spectrum becomes increasingly problematic for numerical (as well as theoretical) studies.

Interest in the small-scale limit of the CDM hierarchy was prompted by the realization that dark matter halos have a wealth of substructure (Moore, et al., 1999a; Klypin, 1999) and that this substructure may have important implications for both direct and indirect dark matter detection experiments (See, for example, Stiff, Widrow, & Frieman (2001); Diemand, Kuhlen & Madau (2007); Kuhlen, Diemand, & Madau (2008) and Kamionkowski & Koushiappas (2008)). High-resolution simulations suggest that the subhalo mass function extends down to the dark matter free-streaming scale with approximately constant mass in substructure per logarithmic mass interval. These simulations probe structures which form from an initial power spectrum with . For example, in the simulation of the first CDM objects by Diemand, Kuhlen & Madau (2007), . With such extreme spectra, care must be taken in order to insure that the results are not corrupted by finite-volume effects. This issue, as it relates to the halo and subhalo mass functions, is discussed in Power & Knebe (2006) and Bagla & Prasad (2006) as well as in the companion to this paper, Elahi et al. (2008).

In this paper, we provide insight into the small-scale limit of CDM by focusing on scale-free cosmologies, that is Einstein-de Sitter cosmologies where the initial power spectrum is a power-law function of , . The guiding principle for understanding structure formation in these models is self-similar scaling which implies that the functional form of the dimensionless power spectrum is time-independent, up to a rescaling of the wavenumber (see below). Self-similar scaling provides a diagnostic test of whether a simulation has sufficient dynamic range (Jain & Bertschinger, 1998). Contact with the standard CDM cosmology is made by treating parameter as a proxy for scale: corresponds to cluster scales and to galactic scales. The limit corresponds to the bottom of the CDM hierarchy.

We perform N-body simulations with and , and compare our results directly with theoretical models. The computation costs to conduct credible simulations with are prohibitively high (see below) and even our results must be considered suspect because of finite box-size effects. We compare our results for the power spectra in the quasilinear regime with predictions from one-loop perturbation theory and demonstrate that for and , the agreement is vastly improved if one implements the renormalization-group approach suggested by McDonald (2007).

To predict the full non-linear power spectrum, one must resort to semi-analytic models such as the ones described in Hamilton et al. (1991) and Peacock & Dodds (1996). These models are based on the stable-clustering hypothesis (Peebles, 1974; Davis & Peebles, 1977) which holds that gravitationally-bound systems decouple from the rest of the Universe once they collapse. An alternative, known as the ‘halo model’ (Ma & Fry, 2000; Peacock & Smith, 2000; Seljak, 2000), allows for the continual accretion of mass onto existing haloes. The density field is treated as a distribution of mass concentrations, each characterized by a density profile. The power spectrum then involves the convolution of this density profile with the halo mass function. Simulations by Smith et al. (2003) of structure formation in a number of scale-free cosmologies demonstrate a clear departure from the stable clustering hypothesis and appear to support the halo model, at least qualitatively.

Peacock & Dodds (1996) and Smith et al. (2003) provide fitting formulae for nonlinear power spectra. Formally, these formulae apply to all initial power spectra with though they are calibrated using simulations with . One goal of this paper is to provide an alternative fitting formula which applies when .

The overall layout of this paper is as follows: In §2, we present background material including a discussion of self-similar scaling, perturbation theory, and semi-analytic models. We describe our simulations in §3 and our results in §4. We also provide a new and improved fitting formula for the nonlinear power spectra. In §5, we discuss the implications of our results for the halo model and stable clustering hypothesis. We conclude, in §6, with a summary and a discussion of directions for future investigations.

2 Preliminaries

2.1 Statistics of the Density field

In keeping with standard definitions (e.g. Peacock 1999), we express real-space density perturbations as deviations from the mean background density, , and then construct the -space representation as follows:


The power spectrum, , (strictly speaking a spectral density) is defined by the relation,


where denotes an ensemble average, is the Dirac delta function, and where the -dependence is dropped for notational simplicity. Note that both and have units of volume. So long as the density field is statistically isotropic, encapsulates the wavenumber dependence. A useful quantity is the dimensionless power spectrum, , which measures the power per logarithmic wavenumber bin. Given these definitions, the wavenumber at which corresponds to the dividing line between the linear and nonlinear scales.

2.2 Initial Conditions

Initial conditions for our N-body simulations are specified at an early epoch when the density field is accurately described by linear perturbation theory. We assume that prior to this epoch, the power spectrum is a power-law function of between appropriately chosen high and low wavenumber cutoffs. That is,




is the scale factor, and is a normalization constant. (Here and throughout, we assume an Einstein-de Sitter cosmology. For other cosmological models, is replaced by the linear growth factor, .) The functions and truncate the power spectrum above and below , respectively. We assume and where is the size of the simulation “box” and is the Nyquist wavenumber of the initial particle distribution. We also follow Kudlicki, Plewa, & Rozyczka (1996) in using the truncation functions


where is the Heaviside function. A similar initial set-up is used in our renormalization group calculations (see section 4.1). Our assumptions imply that at this early epoch , so that all modes within the simulation or computation box are initially in the linear regime.

Utilizing these definitions, we describe the nonlinear scale by the condition , i.e., . Since is the only preferred length scale, the power spectrum should evolve according to the self-similar scaling ansatz


(see Jain & Bertschinger (1996, 1998) for a discussion and review of the literature). Note that since , we can also write as a function of .

2.3 Perturbation Theory

The evolution of the power spectrum from the linear regime to the mildly nonlinear or quasilinear regime can be estimated via perturbation theory (PT, see, for example, Bernardeau et al. (2002)). The starting point is an expansion for the density perturbation field of the form:


where is again the scale factor, characterizes linear density fluctuations, and denote terms of order . So long as the density field is statistically isotropic, the perturbative expansion can be written in terms of the power spectrum:


where . In practice, the series is rarely carried beyond second order in . The calculations are aided by a diagrammatic scheme in which the higher-order terms are described as “loop corrections” to the “tree-level” term, (Scoccimarro & Frieman, 1996). The one-loop correction, , comprises two distinct terms or diagrams, and . These terms involve integrals over the linear power spectrum, . Explicit expressions can be found in Scoccimarro & Frieman (1996), McDonald (2007) and elsewhere. For scale-free models, they combine to give


where , which can be calculated analytically, is positive for and negative for (Scoccimarro & Frieman, 1996; Bernardeau et al., 2002). Hence, represents a “critical index” where nonlinear corrections are vanishingly small (Scoccimarro, 1997).

2.4 Renormalization Group Approach

PT breaks down when the loop corrections, which are formally divergent, become comparable to the tree-level terms. The renormalization group (RG) scheme proposed by McDonald (2007) alleviates this problem essentially by updating the perturbative expansion as the system evolves. RG removes the divergences in the PT expansion, leaving behind a well-behaved, renormalized power spectrum. Operationally, one begins with an initial power spectrum and takes a small step forward in time using one-loop perturbation theory. The new power spectrum is used as the initial condition for the next timestep. This procedure is accomplished by solving the integro-differential equation


with the initial conditions, . Here, and with the proviso that rather than is used in evaluating and (see McDonald (2007)). Details on our own scheme for solving this equation are given in Section 4.1.

As discussed in McDonald (2007), the method has a number of limitations. In particular, Eq. 11 ignores higher-order terms (two-loop and beyond) in . Moreover, decaying mode solutions are not included. Thus, our RG results should be interpreted with a degree of caution.

The approach described here is an example of how RG techniques can be used to remove secular divergences in differential equations. Crocce & Scoccimarro (2006) and Crocce & Scoccimarro (2008) outline an alternative method to study the nonlinear evolution of large-scale structure which also employs RG techniques. Their formalism is conveniently represented in terms of Feynmann diagrams and is closer in spirit to RG applications in high energy and statistical physics.

2.5 Nonlinear Regime — Stable Clustering vs. Halo Model

N-body simulations provide the most direct means to determine the nonlinear power spectrum. The so-called HKLM procedure provides an avenue by which one can begin to understand the simulation results from a theoretical basis (Hamilton et al., 1991; Peacock & Dodds, 1996). The approach yields fitting formulae for the nonlinear power spectrum (or two-point correlation function, ) in different cosmological models. The key assumption is the existence of a one-to-one relation between the nonlinear power spectrum at wavenumber and the linear power spectrum at an earlier epoch and smaller wavenumber, . The starting point is the mapping


which is the analogue of the real-space relation,


that results from associating the volume averaged correlation function, , with the local over-density. The nonlinear power spectrum may then be written as a function of the linear power spectrum at the earlier epoch:


where is an appropriately chosen fitting formula and, by assumption, . Note that for scale-free models, Eq. 12 together with Eq. 7 implies Eq. 14. In the linear regime must be the identity function, i.e., , and the power spectrum grows with the scale factor as . In the nonlinear regime, Hamilton et al. (1991) and Peacock & Dodds (1996) appeal to the “stable clustering hypothesis”, which holds that highly nonlinear structures decouple from the expansion. Under this assumption, and the asymptotic behaviour of must be given by .

To estimate , Peacock & Dodds (1996) advocate the fitting formula


where the parameters , , and are determined by fitting Eq. 15 to power spectra measured in simulations. The parameters are understood as follows: determines a second-order departure from linear growth, and control the behaviour in the quasilinear regime, controls the amplitude of the asymptote and shapes the transition between the two regimes. Best-fit parameters are expressed as functions of , for example, . The expressions for the other parameters similarly diverge as though one must bear in mind that they are based on results from simulations with .

The assumption of stable clustering has been challenged by various groups (Ma & Fry, 2000; Peacock & Smith, 2000; Seljak, 2000; Smith et al., 2003) on the basis that dark matter haloes continually accrete matter and never fully decouple from the rest of the Universe. An alternative approach is provided by the halo model in which the density field is given as a distribution of mass concentrations (haloes) which evolve and have their own internal structure. The two-point correlation function comprises a one-halo term, which is associated with the correlation of mass within a single halo, and a two-halo term, which is associated with the correlation between different haloes (Ma & Fry, 2000; Peacock & Smith, 2000; Seljak, 2000; Smith et al., 2003). Since the power spectrum is the Fourier transform of the two-point correlation function, the associated components of the power spectrum, and , involve integrals over the halo mass function and Fourier-transformed halo density profile. We return to this point in Section 5.

Motivated in part by the separation of components used in the halo model, Smith et al. (2003) construct a fitting formula in the form


for the power spectra measured in their simulations. By construction, dominates the power spectrum in the quasilinear regime and is meant to account for halo-halo correlations while dominates the power spectrum in the nonlinear regime and is meant to account for single halo correlations. However, the model is purely empirical and not calculated directly from the halo model. It is also worth noting that their formula has eight free parameters, three more than that of Peacock & Dodds (1996).

3 Simulations

N-body simulations are carried out with scale-free initial conditions and , , , using the parallel tree-PM code GADGET-2 (Springel, 2005). Initial conditions are generated on a regular grid using a second order-Lagrangian Perturbation Theory (2LPT) code (Crocce, Pueblas, & Scoccimarro, 2006; Thacker & Couchman, 2006). The primary benefit of 2LPT is to reduce the impact of spurious transient modes which arise from the truncation of the perturbative expansion. Since these modes decay, their impact is to delay the time in the simulation at which credible statistics can be calculated (Scoccimarro, 1998; Crocce, Pueblas, & Scoccimarro, 2006). The simulations are run with a softening length of the initial interparticle spacing. GADGET-2 parameters such as the opening angle used in constructing the particle tree and maximum time step criterion were set to their default value.

Table 3 summarizes key features of the simulations used in this study such as the number of simulation particles and the epochs at which the power spectra are measured. The latter are expressed in terms of the ratio where is defined as the scale factor at the epoch when the mode on the scale of the box is equal to the nonlinear scale, that is, . With this definition, .

As approaches , the absence of modes beyond the box scale induces an error in the nonlinear power spectrum. Smith et al. (2003) adopt the ad hoc criterion that the missing variance, , associated with these modes satisfies the condition


where, to a good approximation, with . Figure 1 shows for the outputs of our six simulations. We see that the final two outputs in our high-resolution and simulations and all but the first few outputs in our simulation fail to meet the Smith et al. (2003) criterion. We return to this point below.

Power spectra are calculated using a cubic mesh with side length (Table 1). We set so long as is a power of 2. Here, is the side length of the initial grid of simulation particles. When is not a power of 2 we set equal to the first power of 2 larger than . The power spectra are calculated using piecewise quadratic spline interpolation (Hockney & Eastwood, 1981) and adjusted to account for the strong filtering of this mass-assignment scheme. No correction is made for shot noise.

Table 1: Summary of Simulations
initial scale factor output scale factors
-1 0.0014 0.026, 0.11, 0.19, 0.21, 0.24, 0.30
-2 0.028 0.052, 0.13, 0.23, 0.31, 0.42, 0.49
0.010 0.031, 0.054, 0.17, 0.29, 0.51, 0.67
0.006 0.024, 0.043, 0.063, 0.17, 0.30, 0.36
-2.25 0.009 0.018, 0.021, 0.039, 0.084, 0.221,0.394
-2.5 0.015 0.033, 0.074, 0.15, 0.23, 0.26, 0.29, 0.33

The ratio of the dimensionless power spectrum at the Nyquist frequency, , to that at the box scale, , provides a measure of a simulation’s dynamic range. For a scale-free power spectrum


Thus, if particles are required to achieve a scaling solution over a reasonable range in when , (Jain & Bertschinger, 1998), particles are required at , particles are required for , and particles are required for . We set the particle number for the simulation on the basis of these arguments.

We fully anticipated that self-similar scaling would not be achieved in our simulation. Our run, carried out with illustrates the difficulties that arise as one attempts to simulate highly negative spectral indices. In practice is the largest simulation we can perform within word-addressing limits. It is also clear from this discussion that running an simulation at will yield little improvement over our simulation since, apparently, one requires . Further discussions of the difficulties in simulating scale-free models with can be found in Smith et al. (2003) and Elahi et al. (2008).

4 Results

In Figure 2 we plot the power spectrum from our simulation at the six epochs listed in Table 1. The power at a given wavenumber increases with time while (in this figure, the wavenumber where the power spectrum deviates from the linear form, ), decreases. These results are in close agreement with previous simulations. The cumulative halo distribution is also in good agreement with the expected results (Elahi et al., 2008). In Figure 3 we test the self-similar scaling ansatz by plotting the dimensionless power spectra as a function of for each of our four high-resolution simulations. Taken together, the spectra from different epochs yield a composite power spectrum. Consider, first, the case (upper left panel). The power spectrum covers 7 orders of magnitude in or, equivalently, 3-4 orders of magnitude in . The fact that the composite power spectrum is very nearly a single-valued function of indicates that self-similar scaling is essentially achieved. Note also that for most values of . This result is consistent with PT (See Eq. 10 and note that ).

As with the run, the composite spectra for and show excellent consistency with the scaling hypothesis. However, a departure from self-similar scaling is observed at large in the simulation highlighting the difficulty of simulating spectral indices. Note that the high- feature in some of the early timesteps of our simulation is a remnant of the grid used in setting up the initial conditions. The feature is subdominant to the physical small-scale power at later times in the simulation.

The effect of resolution in achieving self-similar scaling is illustrated in Figure 4 where we compare the spectra from the three simulations. A departure from self-similar scaling is apparent in the simulation and quite severe for . It may be that these simulations are over-evolved, as suggested by Figure 1. In any case, based upon the scaling arguments presented in the previous section, an simulation should be comparable to an simulation. Hence, it is not surprising that the departures from self-similar scaling seen in the right-hand panels of Figure 4 are comparable to those seen in our simulation (lower-right panel of Figure 3. The departure manifests itself in a suppression of power at small or large scales, as expected since power is missing due to the finite size of the simulation volume.

4.1 Perturbative or Quasilinear Regime

We now focus on the quasilinear regime in order to illustrate the improvement renormalization group methods brings to perturbation theory. We implement the RG approach by solving the Eq. 11 assuming an initial scale-free power spectrum with or . The initial spectrum is “evolved” forward in time (or equivalently, scale factor ) using a 4th-order Runge Kutta scheme with an adaptive stepsize (see, for example, Press et al. (1992)). Each Runge Kutta step requires an evaluation of the and integrals. These integrals have the same functional form as those that appear in Scoccimarro & Frieman (1996) except that is replaced by which, in turn, is updated at each step. As with N-body simulations, we must truncate the initial power spectrum at both high and low wavenumbers. Otherwise, the integrals would diverge. Scoccimarro & Frieman (1996) use sharp cutoffs which are convenient for power-law spectra with integer since analytic expressions can be derived. Smooth cutoffs are more manageable for the RG analysis where and must be evaluated numerically (Scoccimarro, private communication).

To make contact with our discussion in Section 2.2 we refer to the IR and UV cutoffs as and , respectively and assume an initial power spectrum of the form




The ratio , which corresponds to the dynamic range of the calculation, is set to while , which determines the sharpness of the k-space cutoffs, is set equal to . For , and have terms of order which cancel, leaving behind a residual term of order . The challenge, numerically, is to determine the surviving terms which can be much smaller than the terms that cancel, especially for large and small . Here, we use the Romberg integration routine from Press et al. (1992).

The solution to Eq. 11 yields an evolutionary sequence for the power spectrum, . Departures from the linear power spectrum increase with beginning at high wavenumber. We evolve until the is roughly equal to the geometric mean of and . Since our dynamic range is a full three orders of magnitude greater than is found in our N-body simulations, finite box effects are much less a concern here. Moreover, our results are insensitive to the form of the cutoff functions, and , since they are derived in a region well inside the computation box.

In Figure 5 we show the measured in the mildly nonlinear regime together with predictions from PT, RG-improved PT, and the Zel’dovich approximation. The latter is discussed in Taylor & Hamilton (1996). Also shown are the fitting formulae of Peacock & Dodds (1996) and Smith et al. (2003). Note that in the case, slowly decreases with increasing for and it is difficult to discern the true self-similar evolution of the power spectrum given that the actual initial power spectrum has a large- cutoff. On the other hand, for and , it is clear that RG does the best job of tracking the power spectrum in the quasilinear regime the RG power spectra does not quite capture the rapid evolution of the measured power spectra. The situation is less clear for where the validity of the simulation is in doubt. The question remains as to whether agreement between RG and the simulations might be improved by refinements in the RG analysis, such as those suggested by McDonald (2007).

4.2 Nonlinear Regime

In Figure 6, we plot the complete power spectrum for our high-resolution, simulation together with the predictions of Peacock & Dodds (1996) and Smith et al. (2003). Also shown is our own fitting formula given by




The form of this formula is motivated by that of Peacock & Dodds (1996) but has one additional parameter to allow for a more general behaviour in the limit. The parameters, derived by performing a nonlinear least-squares fit (see, for example, Press et al. (1992)), are given in Table 2. The lower panel in Figure 6 shows the logarithmic slope . Note that monotonically decreases with near the Nyquist wavenumber.

In Figures 7 and 8, we show and for our high-resolution, and simulations. Again, decreases monotonically in the large- limit. Nevertheless, by design, virtually all fitting formulae (including our own) have a power-law form in the high- limit, that is, as . The lesson is that fitting formula should not be extrapolated to scales below the smallest scales probed by the simulation used in their construction.

Neither the Peacock & Dodds (1996) nor Smith et al. (2003) fitting formulae do a particularly good job of fitting the power spectra from our simulations. For , the Smith et al. (2003) formula provides a reasonable fit up to but decreases too rapidly beyond this point. Conversely, Peacock & Dodds (1996) predict that is constant in the large- limit whereas the measured power spectrum shows a clear decline. The discrepancies between predicted and measured power spectra for are equally severe. By contrast, Eq. 21 provides an excellent fit to the nonlinear power spectra from our high-resolution simulations. And while it has one more parameter than the fitting formula of Peacock & Dodds (1996), it has two fewer than that of Smith et al. (2003).

The case, shown in Figure 9, is difficult to analyse because of the departure from self-similar scaling. In this plot, we truncate the power spectra from different timesteps at large , where the effects of aliasing is apparent, and at small , where the effects of missing large-scale power is apparent. We contend that this procedure yields a roughly continuous curve, which provides a reasonable facsimile of the true power spectrum. The plausibility of this procedure is illustrated in the left-hand panels of Figure 4 where one can imagine carrying out a similar procedure with our , results to yield an approximate form for the power spectrum from our high-resolution simulation.

-1 -0.158 0.181 0.0729 1.571 1.845 2.664
-2 -0.0312 0.690 0.478 1.243 1.266 8.647
-2.25 3.471 4.038 0.348 1.413 1.372 0.659
-2.5 174.0 110.2 4.532 1.492 1.231 0.879
Table 2: Parameters for fitting formula, Eq. 21

5 Halo Model Revisited

In this section, we explore the halo model vis-à-vis our simulation results in more detail. Our focus here is on the high- limit of the nonlinear power spectrum where the one-halo term dominates. The term can be expressed as an integral over the halo mass function, , and the Fourier-transformed density profile of a single halo, . Following Seljak (2000), we use the peak height as the integration variable where is the critical overdensity for spherical collapse ( in an Einstein-de Sitter universe) and is the rms mass overdensity for a spherical region of radius . For scale-free models, . One finds




is a dimensionless form for the halo mass function.

A commonly used fitting formulae for halo density profiles take the form


where and (Navarro, Frenk, & White, 1996; Moore, et al., 1999b; Kravtsov et al., 1998). The characteristic halo scale length, , depends on the halo mass through the relation where is the virial radius and is the halo concentration parameter. Cosmological simulations indicate that the concentration parameter decreases with mass, roughly as a power-law (see, for example, Navarro, Frenk, & White (1996) and Bullock et al. (2001)). With this in mind, we follow Seljak (2000) and Ma & Fry (2000) and write so that .

Our focus is on the power spectrum in the high- limit where the Fourier transform of may be approximated by a step function,


(At wavenumbers above , decreases as but the Heaviside function provides a suitable form for our discussion.) Press & Schechter (1974) and Sheth & Tormen (1999) provide analytic expressions for . In the small- limit (i.e., small or large limit), one finds where () for Press & Schechter (1974) (Sheth & Tormen (1999)). Putting all this together, we find that in the large- limit where


Ma & Fry (2000). Note, however, that is independent of the cusp parameter indicating that the power spectrum in the strongly nonlinear regime is insensitive to the structure of the inner halo.

In Figure 10 we plot our results for the asymptotic slope of the nonlinear power spectrum, , together with those from the Smith et al. (2003) simulations. We make the point of displaying our results as upper bounds on since the slope of appears to be a decreasing function of as (see Figures 6-9). For and , these bounds agree with the quoted values from the Smith et al. (2003) simulations. Furthermore, the functional dependence of on from their fitting formula appears to be consistent with our and results.

Our results suggest that increases with increasing and that . In other words, as , the nonlinear dimensionless power spectrum becomes independent of (i.e., equal power per logarithmic wavenumber bin) just as with the linear power spectrum. Formulations of the halo model with constant, nonzero cannot reproduce this behaviour. To illustrate this point, we plot these predictions for assuming , as in Seljak (2000), and either the Press & Schechter (1974) or Sheth & Tormen (1999) values for . These predictions are inconsistent with our simulation results and those of Smith et al. (2003).

Clearly, the dependence of the concentration parameter on halo mass is central to the development of the halo model. Bullock et al. (2001) devised a toy model to explain the concentration-mass relation seen in simulations. In the case of scale-free cosmologies their model predicts which, when combined with Eq. 27, would seem to yield the desired behaviour for in the limit. The lower panel in Figure 10 shows this prediction. While it does better than constant- versions of the halo model, the predicted tends to lie above the values obtained in the simulation.

6 Summary and Conclusions

The fundamental tenet of the hierarchical clustering scenario is that small-scale objects form earlier than large-scale ones. A corollary of this statement is that individual structures can be identified with a specific wavenumber range of the primordial power spectrum according to their mass. In CDM cosmologies, the logarithmic slope or spectral index of the primordial power spectrum runs from 1 at large scales to -3 at small scales. Thus, as we increase the dynamic range in our simulations and push to smaller and smaller scales, we probe structures that form from density perturbations with a power spectrum approaching . However this limit represents a singular case where the dimensionless power spectrum is independent of scale and structures across a wide range in mass collapse nearly simultaneously. The nature of structure formation changes and the computing requirements for performing simulations increase dramatically.

This work and our companion paper, Elahi et al. (2008), provide inside into the underlying physics of CDM models by considering scale-free cosmologies. We focus here on the nonlinear power spectrum and in Elahi et al. (2008), on the distribution of subhaloes. The evolution of the power spectrum in scale-free cosmologies is remarkably simple — the dimensionless power spectrum, when written as a function of the ratio , is time-independent. Obviously in simulations, the finite computation volume breaks the scale-free nature of the problem and leads to departures from the scaling solution. The dimensionless power spectrum provides a simple test of whether finite volume effects have corrupted the simulation (Jain & Bertschinger, 1998).

Our high-resolution and simulations demonstrate the scaling solution across the simulation volume while showing clear differences with simulations performed at lower resolution. Moreover, our results differ markedly from the the fitting formula provided by Peacock & Dodds (1996) and Smith et al. (2003). A plausible power spectrum for was constructed by stitching together outputs from different timesteps. Though it shows significant, and entirely expected departures from the scaling solution, it represents our best estimate of the power spectrum for models with this small. We summarize our results for our four high-resolution simulations by means of a simple fitting formula. Future work will fill in the gaps in (e.g., and ) and enable use to develop a model for power spectra of arbitrary and therefore arbitrary shape.

The renormalization group improvements to perturbation theory developed in McDonald (2007) represent a promising avenue for studying scale-free models with and likewise, the low-mass limit of the CDM hierarchy. Not surprisingly, the calculations become more difficult as . Our analysis of the and cases confirms McDonald’s claim that RG does lead to an improvement in the predictions of perturbation theory when compared to simulations with the caveat that, as becomes more negative, the RG-predicted power spectrum fails to capture the rapid rise of the power spectrum seen in the simulations. Our analysis for the case is less conclusive but departures in the simulations from the scaling solution suggest that the problem may reside there rather than in the PT analysis. In principle, RG-improved PT can yield a handle on the form of the power spectrum in the mildly nonlinear regime even as .

The halo model provides a theoretical framework for understanding the two-point correlation function and nonlinear power spectrum. Our results, with respect to this model, are somewhat inconclusive. We agree with Smith et al. (2003) that the stable clustering hypothesis of Hamilton et al. (1991) and Peacock & Dodds (1996) fails. On the other hand, the halo model appears to have difficulty reproducing the relation between the asymptotic slope of the power spectrum and . We must therefore settle for an empirical fitting formula for the nonlinear power spectrum.

It is a pleasure to thank P. McDonald, R. Scocciamarro, and C. Orban for useful conversations. We thank R. Smith for carefully reading our manuscript and providing valuable suggestions. We also thank C. Orban for uncovering out an important typo in one of the equations. PJE acknowledges financial support from the Natural Science and Engineering Research Council of Canada (NSERC). RJT and LW acknowledge funding by respective Discovery Grants from NSERC. RJT is also supported by grants from the Canada Foundation for Innovation and the Canada Research Chairs Program. Simulations and analysis were performed on the computing facilities at the High Performance Computing Virtual Laboratory at Queen’s University, SHARCNET, Arizona State University Fulton High Performance Computing Initiative and the Computational Astrophysics Laboratory at Saint Mary’s University.
Figure 1: Missing variance, , as given in Eq.17, for the six simulations described in the text with line types/colours as follows: solid/blue — ; dotted/red — , ; short-dashed/red — , ; long-dashed/red — , ; short-dashed-dot/green — ; long-dashed-dot/cyan — . The output number, corresponds to the values listed in Table 1. For , we show for the last six outputs. The horizontal black curve corresponds to the criterion adopted by Smith et al. (2003).
Figure 2: Power spectrum, as a function of wavenumber for the simulation. Different colours and symbols correspond to different output times as given in Table 3. The sequence, from early to late times is magenta-blue-cyan-green-brown-red, or, alternatively, open square-filled square-open triangle-filled triangle-open circle-filled circle.
Figure 3: Dimensionless power spectrum, as a function of wavenumber . Colours and symbols are the same as in Figure 2. Bottom panel in each quadrant shows the ratio . Upper left — ; Upper right — ; Lower left — ; Lower right — .
Figure 4: Dimensionless power spectrum, , and the ratio for from simulations with different numbers of particles. Black points are power spectra at different timesteps measured in our highest resolution simulation. Superimposed in colour are measurements from the (left) and (right) simulations. As in Figure 2, the sequence, from early to late times is magenta-blue-cyan-green-brown-red or, alternatively, open square-filled square, open triangle-filled triangle-open circle-filled circle.
Figure 5: The ratio as a function of from our four high-resolution simulations as labelled in each panel. Black points are from the simulation with the different symbols representing results from different outputs (from early to late outputs: solid squares, open triangles, solid triangles, open circles, solid circles). Line colours/types are: blue/solid — RG; green/dotted — one-loop; cyan/dot-dashed — Zel’dovich approximation; red/short-dashed — Smith et al. (2003) fitting formula; magenta/long-dashed — Peacock & Dodds (1996) fitting formula.
Figure 6: The ratio as a function of for the , simulation. Shown is the full range in probed by the simulation is shown. Black points are from the simulation. Red curve is the Smith et al. (2003) fitting formula. Magenta curve is the Peacock and Dodds fitting formula. Blue curve is our own fitting formula, Eq. 21. Plotted in the lower panel is the logarithmic slope, of the power spectrum (see text).
Figure 7: Same as Figure 6 but for .
Figure 8: Same as Figure 6 but for .
Figure 9: Same as Figure 6 but for .
Figure 10: Asymptotic slope of the power spectrum, , as a function of the slope, , of the linear power spectrum. The black points are from our simulations while the green points are from the Smith et al. (2003) simulations. In the upper left panel we show the predictions of the halo model (Eq. 27) for and (Press & Schechter). The upper right panel shows the prediction of the halo model assuming and (Sheth & Tormen). The lower left panel shows the predictions from Peacock & Dodds (1996) (magenta curve) and Smith et al. (2003) (red curve). The lower right panel, shows the predictions of the halo model using the prescription from Bullock et al. (2001) for .


  1. affiliation:
  2. affiliation:
  3. affiliation:
  4. affiliation:
  5. affiliation:
  6. More precisely, the logarithmic slope decreases from to logarithmic corrections where is the spectral index of the primordial power spectrum. An analysis of the three-year WMAP data by Spergel (2007) indicates that . However, for the sake of argument, we will set since the difference is not relevant for our discussion.


  1. Bagla, J. S. & Prasad, J. 2006, MNRAS, 370, 993
  2. Bernardeau, F. et al. 2002, Physics Reports, 367, 1
  3. Bullock, J. et al. 2001, MNRAS, 321, 559
  4. Crocce, M., Pueblas, R., & Scoccimarro, R. 2006, MNRAS, 373, 369
  5. Crocce, M. & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533
  6. Crocce, M. & Scoccimarro, R. 2006, Phys. Rev. D, 73, 063519
  7. Davis, M. & Peebles, P. J. E., 1977, ApJS, 34, 425
  8. Diemand, J., Kuhlen, M. & Madau, P. 2007, ApJ, 657, 262
  9. Elahi, P. et al. 2008, arXiv:0811.0206; MNRAS 2009, in press
  10. Hamilton, A. J. S. et al. 1991, ApJ, 374, L1
  11. Hockney, R. W., & Eastwood, J. W., 1981, Computer Simulation Using Particles, New York: McGraw-Hill.
  12. Jain, B. & Bertschinger, E. 1996, ApJ, 456, 43
  13. Jain, B. & Bertschinger, E. 1998, ApJ, 509, 517
  14. Kamionkowski, M. & Koushiappas, S. M. 2008, arXiv:0801.3269
  15. Klypin, A. et al. 1999, ApJ, 522, 82
  16. Kravtsov, A. V. et al. 1998, ApJ, 502, 48
  17. Kudlicki, A., Plewa, & Rozyczka, M. 1996, Acta Astronomica, 46, 297
  18. Kuhlen, M., Diemand, J. & Madau, P. 2008, arXiv:0805:4416
  19. Ma, C. & Fry, J. N. 2000, ApJ, 543, 503
  20. McDonald, P. 2007, Phys. Rev. D, 75, 043517
  21. Moore, B. et al. 1999, ApJ, 524, L19
  22. Moore, B. et al. 1999, MNRAS, 310, 1147
  23. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  24. Peebles, P. J. E., ApJ, 189, L51
  25. Peacock, J. A., 1999, ”Cosmological Physics”, Cambridge University Press, Cambridge.
  26. Peacock, J. A. & Dodds, S. J. MNRAS, 280, L19
  27. Peacock, J. A. & Smith, R. E. 2000, MNRAS, 318, 1144
  28. Power, C. & Knebe, A. 2006, MNRAS, 370, 691
  29. Press, W. H. et al. 1986, Numerical Recipes, (Cambridge: Cambridge University Press)
  30. Press, W. H. & Schechter, P. 1974, ApJ, 187, 425
  31. Scocciamarro, R., 1997, ApJ, 487, 1
  32. Scocciamarro, R., 1998, MNRAS, 299, 1097
  33. Scoccimarro, R. & Frieman, J. A. 1996, ApJ, 473, 620
  34. Seljak, U. 2000, MNRAS, 318, 203
  35. Sheth, R. K. & Tormen, G. 1999, MNRAS, 308, 119
  36. Smith, R. E. et al. 2003, MNRAS, 341, 1311
  37. Spergel, D. N. et al. 2007, ApJS, 170, 377
  38. Springel, V. 2005, MNRAS, 364, 1105
  39. Stiff, D., Widrow, L. M., & Frieman, J. 2001, Phys. Rev. D, 64, 083516
  40. Taylor, A. N. & Hamilton, A. J. S. 1996, MNRAS, 282, 76
  41. Thacker, R. J. & Couchman, H. M. P., 2006, Int. J. High Perf. Comp. & Net., 4, 303
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description