Exponents of non-linear clustering in scale-free one dimensional cosmological simulations
aa r X i v : . [ a s t r o - ph . C O ] S e p Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 21 October 2018 (MN L A TEX style file v2.2)
Exponents of non-linear clustering in scale-free onedimensional cosmological simulations
David Benhaiem , Michael Joyce and Fran¸cois Sicard Laboratoire de Physique Nucl´eaire et Hautes ´Energies, Universit´e Pierre et Marie Curie - Paris 6, CNRS IN2P3 UMR 7585,4 Place Jussieu, 75752 Paris Cedex 05, France Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR 6303 CNRS-Universit´e de Bourgogne,9 Avenue A. Savary, BP 47 870, F-21078 Dijon Cedex, France
21 October 2018
ABSTRACT
One dimensional versions of dissipationless cosmological N-body simulations have beenshown to share many qualitative behaviours of the three dimensional problem. Theirinterest lies in the fact that they can resolve a much greater range of time and lengthscales, and admit exact numerical integration. We use such models here to study hownon-linear clustering depends on initial conditions and cosmology. More specifically,we consider a family of models which, like the three dimensional EdS model, lead forpower-law initial conditions to self-similar clustering characterized in the strongly non-linear regime by power-law behaviour of the two point correlation function. We studyhow the corresponding exponent γ depends on the initial conditions, characterized bythe exponent n of the power spectrum of initial fluctuations, and on a single parameter κ controlling the rate of expansion. The space of initial conditions/cosmology dividesvery clearly into two parts: (1) a region in which γ depends strongly on both n and κ and where it agrees very well with a simple generalisation of the so-called stableclustering hypothesis in three dimensions, and (2) a region in which γ is more or lessindependent of both the spectrum and the expansion of the universe. The boundaryin ( n, κ ) space dividing the “stable clustering” region from the “universal” region isvery well approximated by a “critical” value of the predicted stable clustering expo-nent itself. We explain how this division of the ( n, κ ) space can be understood as asimple physical criterion which might indeed be expected to control the validity ofthe stable clustering hypothesis. We compare and contrast our findings to results inthree dimensions, and discuss in particular the light they may throw on the questionof “universality” of non-linear clustering in this context. Key words:
Cosmological structure formation, gravitational clustering, N -body sim-ulation Cosmological N body simulations are the primary instru-ment used to make theoretical predictions for structure for-mation in current models of the universe. Analytical un-derstanding of the non-linear regime, crucial for many non-trivial tests of these models, remains very poor despite themany phenomenological results derived in simulations. Bet-ter understanding of the physics of non-linear clusteringcould help to better constrain and control cosmological sim-ulations, which despite their ever increasing power and so-phistication are still subject to questions about their ac-curacy and resolution. In this work we use simplified onedimensional (1D) “toy models” of the full three dimensional(3D) problem to study non-linear structure formation. Morespecifically we study, using a family of such models, how the properties of non-linear clustering depend on both initialconditions and cosmology. Our main result is that, for coldinitial conditions, we observe that the large phase space ofinitial conditions and cosmology we explore breaks up intotwo parts: on the one hand, a region in which there is anapparent “universality” of the non-linear clustering, i.e., inwhich it depends very weakly both on initial conditions andon the cosmological expansion; and on the other hand, aregion where the non-linear clustering depends strongly onboth the initial conditions and expansion history. In the lat-ter region, the exponent characterizing the two point cor-relations is in excellent agreement with a simple analyticalprediction which follows from the so-called stable clusteringhypothesis, appropriately generalized to this class of mod-els. Further the apparently well localized boundary dividing c (cid:13) D. Benhaiem, M. Joyce and F. Sicard the phase space into these two disjoint regions can in fact becharacterized using the stable clustering prediction, in termsof a critical value of the corresponding exponent. The non-linear clustering appears to be strictly universal in the limitin which there is no expansion, but also very well approx-imated for steeper initial spectra when there is expansion.We compare our results with those in three dimensions, anddiscuss how our results may provide a framework for under-standing the much discussed apparent “universality” of haloprofiles in this context.The exploration of 1D models as a tool for under-standing the physics of non-linear structure formationin a cosmological context goes back at least as far asthe work of Melott (1982, 1983) for the case of hotdark matter models. Different initial conditions and vari-ants of the model have been discussed by a number ofother authors (Yano & Gouda 1998; Miller & Rouet 2002;Aurell & Fanelli 2002; Miller & Rouet 2006; Valageas 2006;Miller et al. 2007; Gabrielli et al. 2009; Miller & Rouet2010a). There is also an extensive literature on the statisti-cal mechanics of finite one dimensional self-gravitating sys-tems (see e.g. Joyce & Worrakitpoonpon (2011) and refer-ences therein), as well as some studies which use this case toexplore issues in cosmological N body simulations: Binney(2004) uses it to probe discreteness effects, and Schulz et al.(2012) the issue of universality in cold collapse.The study reported here is a direct continuation of thatin Joyce & Sicard (2011). This latter paper showed, for thethree 1D models studied in the previous literature, the verystrong qualitative similarities in the temporal developmentof clustering to that in the 3D case. For cold cosmologicallike initial conditions clustering is hierachical and driven bythe linear amplification of fluctuations, leading, for the caseof power law initial conditions, to self similar evolution ofthe non-linear correlations. The exponent characterising thenon-linear clustering was, further, found to be in very goodagreement with that predicted by the so-called “stable clus-tering” hypothesis. In this paper we extend this study to aone parameter family of models which allows continuous in-terpolation between the static model and the two expandingmodels previously considered. This allows us to probe fullythe extent of validity of stable clustering, where it breaksdown and what happens in this case. The study reveals thatthere is a fairly abrupt switchover from the validity of stableclustering to an apparent universality in the clustering.In the next section we define the class of 1D models westudy, discussing also their relation to 3D cosmological mod-els. This presentation improves and simplifies that given inJoyce & Sicard (2011). For the family of scale-free models wefocus on, we derive, more rigorously than in Joyce & Sicard(2011), the prediction for the behaviour of an isolated struc-ture and use it to derive the exponents predicted for thenon-linear two point correlation when the stable clusteringhypothesis is made. In the following section we discuss ournumerical simulations, and in particular how we calibratethem using “exact” simulations. The next section gives ourresults, and in the last section we discuss their interpretationand possible relation to the three dimensional problem.
The class of 1D models we study are defined as follows. Wetake the equations of motion in comoving coordinates of 3Dcosmological N body simulations, and simply replace the3D Newtonian gravitational interaction with the 1D gravi-tational interaction. This gives 1D equations of motion d x i dt + 2 H dx i dt = − ga J X j = i sgn( x i − x j ) (1)where x i are the particle positions on the line. The super-script ‘J’ in the force term indicates that the sum, whichextends over an infinite distribution of masses with non-zeromean mass density, is regularized, just as in 3D, by subtrac-tion of the contribution of the mean mass density, i.e., theforce term is sourced only by fluctuations about the meandensity. As we will discuss further below this sum can bewritten explicitly in a simple manner.The coupling g is the 1D analogy of Newton’s constant(multiplied by the particle mass). To establish a more directrelation between the 1D and 3D models, one can considerthe “particles” in the 1D system to represent infinite parallelsheets of zero thickness embedded in 3D. In this case one hasthe identification g = 2 πG Σ where Σ is the surface massdensity of the sheets. The mean (comoving) mass density ρ of the corresponding 3D universe and the number density n of the 1D system are then related as ρ = Σ n , and therefore gn = 2 πGρ . (2)For any given 3D cosmological model which provides a scalefactor a ( t ) and associated Hubble constant H ( t ) we thenhave a unique corresponding 1D model. The case a ( t ) = 1(and H ( t ) = 0) on the other hand, defines a model analogousto a 3D universe without expansion (which is well definedbut does not correspond to a cosmological model derivedfrom general relativity).The one parameter family of models we will study isa continuous interpolation between the 1D model obtainedusing the 3D EdS cosmology, and the static model . To definethem precisely, and motivate their choice, it is convenientto change time variable in Eq. (1) defining τ = R dta / . Thisgives the equations of motion in the form d x i dτ + Γ( τ ) dx i dτ = − g J X j = i sgn( x i − x j ) , (3)i.e., in which all effects of the cosmology appear only as afluid type damping term, withΓ = 12 a / H = 12 a dadτ . (4)For a 3D EdS universe with comoving mass density ρ , i.e., a = ( t/t ) / where t = 1 / √ πGρ , we obtain a Γ which is independent of time and given byΓ = 13 t = p πGρ / p gn / c (cid:13) , 000–000 xponents of non-linear clustering in one dimension the dependence of non-linear clustering on this parameter,as well as on the initial conditions. Γ can simply be con-sidered as a control parameter for understanding the roleof expansion in the determination of the properties of thenon-linear clustering. While most previous studies have con-sidered either the EdS model (also known as the “quintic”model for reasons which will be recalled below), or the staticmodel, one other model in this family, corresponding to thecase Γ = √ gn , introduced in Miller & Rouet (2002), hasbeen studied in Miller & Rouet (2006); Miller et al. (2007);Miller & Rouet (2010a) as well as in Joyce & Sicard (2011).We note, using Eq. (4), that this class of models is ob-tained by taking the functional dependence of the 3D EdSexpansion law, but allowing a freedom in the normalisationof the expansion rate to the matter density, i.e., we can ob-tain these models taking H = κ πGρ / a where κ is apositive constant. Eq. (4) then still holds, and thus a = e τ ,but instead of (5) we have, using Eq. (2), thatΓ = κ p gn / . (6)In terms of 3D cosmological models the case κ >
1, i.e.,Γ > Γ EdS , can be considered to be a model in which thereis not just “ordinary” matter but an additional pressure-less component of the energy density which does not clus-ter, e.g., a homogeneous scalar field with appropriate equa-tion of state (see e.g. Ferreira & Joyce (1998) and referencestherein). Equivalently one can consider that the κ = 1 mod-els are obtained using the “correct” normalisation for theHubble constant, H = 8 πGρ /
3, but choosing the 1D cou-pling freely instead of imposing the relation Eq. (2). It isimportant to note that, given that 3D Hubble law is im-posed “by hand” on the 1D model, there is in principle nocorrect value for Γ. Our point of view here is simply to useΓ (or κ ) as a control parameter in trying to understand therole of expansion in non-linear structure formation. In otherwords we are seeking to improve our qualitative understand-ing only through the study of these toy models. In the context of the problem of structure formation thefamily of models we study has the interest of being, likeEdS models and static models in 3D, scale-free: the expan-sion introduces no additional time or length scales. Thishas as a consequence that, for cold particles with initialdensity fluctuations with a power law spectrum, one ex-pects to obtain, at sufficiently long times, temporal evo-lution which is “self-similar”, i.e., the evolution in time ofthe spatial correlation functions is equivalent to a simpletime dependent rescaling of the spatial variables. The rea-son why this occurs is that there is, in this case, only asingle physically relevant length scale in the problem: thescale at which fluctuations are of order one, of which theevolution is predicted by linear theory . Further if the clus-tering in the strongly non-linear regime does not dependon this scale one would expect to obtain a scale-free be-haviour of spatial correlations in this regime. Such power In the discrete problem there is at least one additional lengthscale — the initial interparticle spacing — but physical (contin-uum limit) results should not depend on it. law behaviour has been observed in both expanding (see,e.g., Efstathiou et al. (1988); Smith et al. (2003) and refer-ences therein) and static 3D simulations (Baertschiger et al.2007a,b, 2008) for a range of power law initial conditions. 3Dsimulations, however, can show such behaviour over a verylimited spatial range, making it difficult to establish if it isassociated with a true scale invariance of the non-linear clus-tering. The study of 1D models in the family we are consider-ing (Miller & Rouet 2006; Miller et al. 2007; Miller & Rouet2010a; Joyce & Sicard 2011), with much greater spatial res-olution and numerical accuracy, shows very convincing ev-idence that these models do indeed give rise to power-lawclustering indicative of truly scale-invariant clustering.A central problem of non-linear structure formation isthen that of understanding how the exponent (or possiblyexponents) characterizing the non-linear clustering is deter-mined. In 3D EdS cosmology Peebles (see Peebles (1980))proposed many years ago a simple explanation of how suchscale invariant clustering might arise: if highly non-linearstructures (of all different sizes) are supposed to evolve es-sentially independently after their formation, they will viri-alize and remain stable in physical coordinates. From thisone may derive an exponent which depends only on theexponent in the power spectrum of initial density fluctua-tions. We derive here now, more rigorously than the heuristicderivation given in Joyce & Sicard (2011), the prediction ofsuch a stable clustering hypothesis for the 1D models westudy. To do so we simply derive and then analyse the equa-tions of motion of a finite sub-system.
In the simulations which we consider here, we start, as in 3Dcosmological simulations, from a particle configuration gen-erated by applying small displacements to an infinite perfectlattice configuration. At any moment in time such a config-uration can be fully specified by giving the displacements u i of particles from the sites of a perfect lattice. The positionof the i -th particle with respect to some arbitrary origin, is x i = iλ + u i , where λ = n − . In a particle labelling with x i > x i − for all i , the displacements are uniquely defined,and such that no two particles cross when they are appliedto the lattice. In this case it has been shown (Gabrielli et al.2009) that the force on the particle i , given by the sum onthe right side of Eq. (3), with an appropriate specificationof the regularisation, has the simple exact expression − g J X j = i sgn( x i − x j ) = 2 gn ( u i − h u i ) (7)where h u i is the average particle displacement. This resulthas been shown for an infinite system with statisticallytranslationally invariant displacements in Gabrielli et al.(2009), with the sole assumption of decay of correlations ofthe displacements at asymptotically large scales. The sameresult has been shown to hold in Miller & Rouet (2010b) forthe case of an infinite periodic system, using an appropriateregularisation of an Ewald summation. In this case, as wewill now show explicitly, the expression (7) can be written c (cid:13) , 000–000 D. Benhaiem, M. Joyce and F. Sicard as − g J X j = i sgn( x i − x j ) = − g [ N L ( i ) − N R ( i )] + 2 gn [ x i − h x i p ](8)where N L ( i ) ( N R ( i )) is the number of particles in the pe-riodic cell to the left (right) of particle i , and h x i p is theaverage particle position, i.e., centre of mass of the particlesin the cell .Let us consider now any finite subsystem S of the infi-nite system defined simply by the N particles, with i = 1 ...N say, in some interval [ x , x + L ] at a given instant in time,Substituting u i = x i − iλ in (7), and noting that the numberof particles in the interval to the left (right) of the i th par-ticle is given by N L ( i ) = i − N R ( i ) = N − i ), we obtainthat the force on the particle i may be written F ( x i ) = − g [ N L ( i ) − N R ( i )]+2 gn [ x i −h x i S ] − gn [ h u i−h u i S ](9)where h x i S is simply the position of the centre of mass ofthe subsystem, and h u i S the average displacement of theparticles in the subsystem. If the subsystem chosen is theperiodic cell of an infinite periodic system, we have h u i = h u i S and thus obtain (8). For any generic subsystem S , onthe other hand, the last term on the right in (9) is the samefor all particles in the subsystem. Defining y i = x i −h x i S wethus obtain the equation of motion for any particle in thesubsystem, relative to the CM of the subsystem, as d y i dτ + Γ dy i dτ = − g X j = i,j ∈S sgn( y i − y j ) + 2 gn y i (10)where the sum extends over the particles in the subsystem.This remains valid for as long as no other particle outsidethe system crosses the particles at its extremities. Eq. (10)are thus exactly the equations of motion of a subsystem aslong as it is isolated in this sense. For Γ = 0 we can transform to the coordinates˜ t = 3Γ e Γ τ/ , r i = e τ/ y i (11)in which Eq. (10) becomes d r i d ˜ t = − g X j = i,j ∈S sgn( r i − r j ) + 2˜ t [1 + 9 gn Γ ] r i . (12)These coordinates have been chosen in order to remove thedamping term in the equations of motion, and all depen-dence on the expansion now appears only in a time depen-dent force term, of which the amplitude decays as 1 / ˜ t . Atsufficiently long times, therefore, an isolated subsystem ofthe infinite 1D “expanding” model has, to an arbitarily goodapproximation, the equations of motion of a strictly isolatednon-expanding 1D self-gravitating system. We will refer tothe coordinates in (11) as quasi-physical coordinates : theyare the closest analogy we have to what are called physicalcoordinates in the three dimensional problem (i.e. ~r = a ( t ) ~x , Contrary to what is stated in Miller & Rouet (2010b), the ex-pressions for the force in the two papers are in agreement. t ). In these coordinates the 3D cosmological equations ofmotion are equivalent to those for an infinite Newtonianself-gravitating system expanding about a centre in an in-finite (static) space. Correspondingly, any subsystem (e.g.galaxy) taken far from other mass in an expanding universewill evolve in these coordinates like an isolated Newtoniansystem (see e.g. Joyce & Sylos Labini (2012) for a deriva-tion). The additional term in 1 / ˜ t in the 1D model, whichwe cannot get rid of by a coordinate transformation, appearsbecause we have imposed in the derivation of the modelthe 3D expansion law rather than using the analogy of thislaw in 1D (which gives a collapse in a finite time, qualita-tively different to that of the 3D expansion which we wish tomodel). It is important to note that, re-expressed in terms ofthe three dimensional scale factor a ( t ), the transformationin (11) from the “comoving” coordinates y i to the “quasi-physical” coordinates in (11) is r i = e τ/ y i ≡ a / y i , and not r i = ay i as one might naively expect . The second term on the right hand side of Eq. (10) may berewritten simply as [ κ + 1] gn y i The first term is of or-der gN ∼ gn S L where n S is the average comoving densityof the subsystem. Thus if we consider any isolated subsys-tem which is already highly non-linear (with n S ≫ n ) weare necessarily in the regime where the second term maybe neglected. In this case the isolated subsystem will thusbehave, in quasi-physical coordinates, like a truly isolated1D self-gravitating system. Just as in 3D, such a system, forlarge N , in the collisionless limit, attains on a few dynamicaltime scales ( ∼ / √ gn S ) a stationary virialized state (see e.g.Joyce & Worrakitpoonpon (2011) and references therein),i.e., a system which becomes time independent in the quasi-physical coordinates r i . More generally, even if not strictlyvirialized, it rapidly attains a well defined average size aboutwhich virial oscillations persist. In order to derive the exponent in what follows we will alsouse the self-similar behaviour of the evolution, which, for thetwo point correlation function ξ ( x, τ ) means that ξ ( x, τ ) = ξ (cid:16) xR s ( τ ) (cid:17) (13)in the spatial range where the self-similarity is observed .The function R s ( τ ) is derived from linear theory, which de-scribes a self-similar evolution from power law initial con-ditions once the growing mode dominates. In the 1D model We note that Yano & Gouda (1998) give an incorrect general-ization of the stable clustering hypothesis to the 1D EdS modelassuming that stability will be attained in the coordinates ax . Asa result they arrive at the incorrect conclusion from their numer-ical study that stable clustering is not observed in the model. We follow the standard notation here in which ξ ( x ) de-notes the reduced two point correlation function, i.e., ξ ( x ) = n − h n (0) n ( x ) i − n ( x ) is the local particle density and h· · · i denotes the average over the statistically translationally in-variant point process. See e.g. Peebles (1980) or Gabrielli et al.(2005). c (cid:13) , 000–000 xponents of non-linear clustering in one dimension ξ ( x ) x/LX min X max ξ (x min ) ξ (x max ) Power law fit ξ (x,t) Figure 1.
Two point correlation function measured in an evolvedsimulation, with κ = 1 and n = 2. The values of x min and x max define the range in which there is (by eye) an apparent power-lawbehaviour. the linear theory growth of perturbations can be easily de-rived by using the expression for the force Eq. (7) in termsof the displacements u i . Setting, without loss of generality,the mean displacement h u i to zero, we have simply d u i dτ + Γ du i dτ = 2 gn u i . (14)which is characterized by the growing mode u i ∝ e α Γ τ = a α where α = 14 [ − r κ ] . (15)For κ = 1 we thus recover the same growth law as inthe standard 3D EdS model. Assuming a power law powerspectrum of density fluctuations P ( k ) ∝ k n (and thus ξ ∝ /x n +1 at large x ) the evolution is self-similar with R s = e α n Γ τ = a α n . (16) Shown in Fig. 1 is a typical two point correlation ξ ( x ) ob-tained from the evolution of power law initial conditions inthe class of models we study. Between a lower cut-off scale x min , and an upper cut-off scale x max a simple power lawbehaviour of the correlation function is observed ξ ∼ x − γ .The exponent γ may be written as γ = − ln ( ξ ( x max )) − ln ( ξ ( x min ))ln ( x max ) − ln ( x min ) . (17)The temporal evolution of the correlation function isfound to be self-similar, down to the scale x min , below whichit is broken. The scale x max is in the range in which self-similarity applies, and thus scales in proportion to R s , with ξ ( x max ) a fixed amplitude corresponding to the transitionbetween the linear and non-linear regime.If we make now the hypothesis that the structures giv-ing rise to the correlation measured at the lower cut-off are stable in quasi-physical coordinates , i.e., these structures be-have like isolated approximately virialized structures, wehave that x min ∼ a − / . Further, since for ξ ≫ ξ ( x )mesures simply the average density at distance x from anoccupied point normalized to the mean density, we have that ξ ( x min ) ∼ a / . It follows that the exponent tends asymp-totically (at sufficiently long times) to γ sc ( n, κ ) = 2 κ ( n + 1) κ (2 n −
1) + 3 √ κ + 24 (18)This result is identical to that derived in Joyce & Sicard(2011) using more heuristic arguments based on the scalingof the energy . It tells us how the exponent, in the sta-ble clustering hypothesis, depends on the initial conditions,characterized by the exponent n of the power spectrum ofdensity fluctuations, and on the expansion rate parametrizedby the constant κ , which we recall is κ = Γ p /gn = (cid:18) H a πGρ (cid:19) / . (19)Note that as κ → γ sc = 0, i.e. a correlationfunction which is predicted to be flat over an arbitrarilylarge distance. Indeed, in this static limit, the stable clus-tering hypothesis corresponds to stability also in comovingcoordinates, and therefore to a constant x min and ξ ( x min ).Clearly, however, we do not expect in this case that the for-mation of non-linear structure at a given scale will leave ap-proximately unperturbed the structures formed previouslyat smaller scales. The numerical study in Joyce & Sicard (2011), as most pre-vious studies of the three models belonging to the class weare considering, exploits fully the very attractive propertyof these models that they admit “exact” numerical integra-tion: the force on particles, whether given in the form (7) or(8), is such that the motion (with or without damping) ad-mits an analytical solution other than at particle crossings.An event driven algorithm can then be employed, in whichthe determination of the time of the next “collision”, andwhich pair of particles it involves, requires only the solutionof algebraic equations. As described by Noullez et al. (2003)the integration can be sped up optimally using a “heap”structure.For the previously considered cases, corresponding to κ = 0 , , √
3, the algebraic equations involved are, respec-tively, quadratic, quintic and cubic. For a generic κ , however,the equations are not polynomial, and the numerical cost ofsolving them increases and makes the code, which is already In Joyce & Sicard (2011), as in various other papers on thesemodels, units were used in which 2 gn = 1, and thus Γ = κ/ √ (cid:13) , 000–000 D. Benhaiem, M. Joyce and F. Sicard expensive numerically, even more time consuming . As wewish to study fully a broad range of κ and initial conditions,we have chosen here instead to use a particle-mesh (PM)code of which the numerical speed can be greatly enhancedby the use of FFT techniques and parallelization . In orderto be sure that we are not, as a result, losing the advantagesof the precision of the integration accessible in these 1Dmodels, we calibrate, as discussed below, our choice of thenew discreteness parameters introduced by the PM code bycomparing our results with those obtained, for κ = 0 , , √ n ( x ) is defined at each time stepby a “clouds in cell” interpolation of the particle positions.The modified 1D Poisson equation d φd x = 2 g [ n ( x ) − n ] (20)is then solved for the potential in Fourier space using anFFT, using the solver fftw3 . An inverse FFT then deter-mines the force on the grid, which is then interpolated backonto the particle positions. The latter are then advanced us-ing a leapfrog algorithm, with a simple adaptative timestep:it is chosen at each time by imposing that the mean distancetravelled by particles during a time step be equal to a chosenfraction η of the PM grid.Compared to the exact code, our PM code therefore in-troduces resolution effects controlled by two dimensionlessparameters: the size of the PM grid compared to the latticespacing in the initial conditions, which we denote by ǫ , andthe parameter η controlling the size of the time steps. Wehave made an extensive study of the effects introduced bythis finite resolution of the PM code, and in particular on thedetermination of the two point correlation function which isthe quantity which interests us here. Shown in Fig. 2 is, forexample, a comparison of the results for runs from identicalinitial conditions, for the case n = 2 and κ = 1, obtainedfrom the exact code and the PM code, for two different val-ues of ǫ , and η sufficiently small (of order unity) so thatconvergence is obtained with respect to it. We find, very rea-sonably, that excellent agreement is obtained provided ǫ ischosen smaller than x min , while a larger ǫ than this minimalscale leads to a visible deviation of the correlation functionfrom its correct value at scales below ǫ . We have chosen ournumerical parameters in the simulations reported below inorder to obtain such an indiscernible difference in the cor-relation function between our code and the exact code, forthe cases κ = 0 , , √
3. All the results presented here arefor systems with N = 10 particles (in periodic boundaryconditions).For initial conditions with power law power spectra P ( k ) ∼ k n (at small k , a cut-off at large k is always im-plicitly assumed), we expect hierarchical structure forma-tion in d dimensions for n in the range − d < n n − d the variance of mass fluctuations diverges The equation for the crossing times is of the form 1 + Az β + Bz ( β +1) / = 0 where z = e Γ τ , β = q κ , and A and B areconstants. See, e.g., http://astro.uchicago.edu/ ∼ andrey/Talks/PM/pm.pdfand references therein. ξ ( x ,t ) x/L exact simulation ξ (x,t) for ε =0.76 ξ (x,t) for ε =0.012 Figure 2.
Two point correlation functions obtained by evolvingfrom the same initial condition (with n = 2) and κ = 1, for threedifferent cases: using the “exact” code (solid line) and using ourPM code with the two indicated values of the PM grid spacing. at large scales, while for n > n = 0 , ,
4. Point processeswith such power spectra are easy to set up: n = 0 is obtainedby randomly distributing points in the interval (giving apower spectrum P ( k ) = 1 /n ), the case n = 2 by applyingsmall random and uncorrelated displacements to a regularlattice, and n = 4 in the same way but with the additionalconstraint that neighbouring pairs have equal and oppositedisplacements (see Gabrielli & Joyce (2008) and referencestherein). Our results are summarized in Fig. 3. Each point on theplot corresponds to a single simulation with given n and κ , where the latter varies in the range 0 —2 .
5. The figureshows then the value of the exponent γ of the power law re-gion of the correlation function. In each case the given valueis determined using a linear regression in the range between x min and x max , the latter being determined themselves byeye (cf. Fig. 1). The uncertainty in the fitted exponent as-sociated with the (by eye) estimation of the end-points forthe fit is very small (of order a few percent at most) inmost parts of the parameter space. However, as we will nowdiscuss and explain, at lower values of κ and n , where theextent of the power law region is much smaller, the resul-tant uncertainty becomes more significant. The continuouscurves plotted correspond to the predictions of the stableclustering hypothesis, Eq. (18).Our results show excellent agreement with the stableclustering prediction for a very large part of the exploredparameter space of initial conditions and cosmology. Thisregion can be well characterized simply by the condition c (cid:13) , 000–000 xponents of non-linear clustering in one dimension γ κ n=0 γ n=0 γ SC n=2 γ n=2 γ SC n=4 γ n=4 γ SC Figure 3.
Exponents of non-linear power law clustering γ measured in simulations with different initial conditions(parametrized by n ) and expansion rates (parametrized by κ ).The solid lines are the predictions of stable clustering, Eq. (18). γ sc ( n, κ ) ∼ > .
2, i.e., good agreement with stable cluster-ing is observed simply in the region where the predictedexponent be larger than some critical value. In the rest ofthe parameter space (i.e. for γ sc ( n, κ ) ∼ < .
2) there is cleardisagreement with the stable clustering prediction, and themeasured exponents lie in a very narrow range, between 0 . .
15, or a little smaller. It is precisely in this part of theparameter space, however, that there is also a considerablescatter about the curves of the measured exponents. Thisscatter is in fact just a measure of the uncertainty in the de-termination of these exponents, which is most difficult (forreasons we explain in detail below) in the region where theexponents become small. Thus, our results are quite con-sistent with the hypothesis that the exponent in the region γ sc ( n, κ ) ∼ < . γ sc ( n, κ ) = γ . Rather than a line, there may of coursebe a small region in parameter space in which the exponentvaries in a continuous manner from its stable clustering valueto a limiting universal value of γ ≈ . . In the language of statistical physics our result could be de-scribed as an evidence for an “out of equilibrium phase transi-tion”: as control parameters for the system are changed the (outof equilibrium) behaviour of the system changes in a discontinu-ous manner. With the precision of our current results we cannotdetermine with confidence whether there is really such a transi-tion or rather a continuous “crossover” between the two “phases”. κ nUniversality Stable Clustering Figure 4.
Schematic representation of our results: the space ofinitial conditions ( n ) and cosmology ( κ ) breaks into two parts,with a boundary defined approximately by γ sc ( n, κ ) = 0 . We have found numerically that the stable clustering hy-pothesis predicts the exponents of non-linear clustering veryaccurately at sufficiently large n and κ , i.e., if the spectrumof initial fluctuations is “sufficiently blue” and the expansionrate is sufficiently fast. More precisely it appears that thecriterion for it to work is that γ sc , which is a monotonicallyincreasing function of both n and κ , be larger than somecritical value. It is not difficult to give a simple physicalexplanation for this behaviour, as we now explain.Let us consider two overdense regions, initially (in thelinear regime) of comoving size L and L > L . Self-similarity of their evolution implies that there is simply atime delay ∆ τ between every stage of their evolution givenby R s (∆ τ ) = ( L /L ). Using Eq. (16) we obtain∆ τ = 1 + n α ln( L /L ) . (21)Let us suppose that the fluctuation 1 goes non-linear andapproximately virializes at some given time τ v , when its sizeis L v . When fluctuation 2 virializes, at time τ v = τ v + ∆ τ ,its size is L v = ( L /L ) L v . If in the interval between τ v and τ v , structure 1 behaves to a very good approximation asan isolated structure, its size in comoving coordinates willdecrease by a factor e − τ/ . It follows that the relativesize of the two structures at the time when 2 virializes is (cid:18) L L (cid:19) = (cid:18) L L (cid:19) − n α (cid:18) L L (cid:19) (22)where α is the growing mode exponent given in (15). Quitesimply, the faster the expansion rate, or the larger is n , themore a given scale which has entered the non-linear regimecan “shrink” relative to a larger scale in the time beforethe latter goes non-linear. Further we note, using (15), that c (cid:13) , 000–000 D. Benhaiem, M. Joyce and F. Sicard
Eq. (22) can be written (cid:18) L L (cid:19) = (cid:18) L L (cid:19) − γ sc1 − γ sc (cid:18) L L (cid:19) (23)i.e., the dependence on ( n, κ ) is in fact completely deter-mined through the exponent γ sc itself.The predicted stable clustering exponent is therefore di-rectly related to a physical quantity which we would expectnaturally to control the validity of the assumption madein deriving it: the more “concentrated” are the pre-existingvirialized substructures inside a larger structure when it col-lapses, the better should become the approximation that thisstructure will behave simply as a collection of sub-structureswhich are not disrupted by the virialization and subsequentevolution of the larger structure. Indeed, in the limit that γ sc approaches unity, any structure which collapses and virial-izes will see the substructures it contains essentially as pointparticles. As γ sc decreases, on the other hand, we expectthat the interaction between structures can lead to theirdisruption, and in particular that “mergers” of substruc-tures become much more probable. Given that the value ofthe exponent γ sc characterizes precisely the relative “con-densation” of scales, it is very natural that a critical valueof the exponent should characterize the breakdown of stableclustering.What our results do not allow us to conclude, as we havediscussed, is what precisely happens in the regime wherestable clustering breaks down: our numerical results do notallow us to clearly distinguish between the possibility of anabrupt (discontinuous) transition to a region in which thereis a truly universal exponent, or a smoother transition to aregion in which the exponent depends only very weakly onthe relevant parameters ( n, κ ). Further we do not (currently)have a model to explain the exponent (or narrow range ofexponents) observed in this part of the parameter space.Extending the analysis given just above to determinethe scalings of relative sizes of structures, it is simple to un-derstand why it is precisely the region of small γ sc in whichthe numerical results for the exponent of the power law inthe correlation function are most noisy. The degree of preci-sion in the measured exponent is essentially just a functionof the range of scale over which the power-law behaviour ex-tends, i.e., it depends on the ratio x min /x max accessible in anumerical simulation. The scale x max corresponds approxi-mately at any time to the size of the largest (approximately)virialized scale. Let use denote by x i the comoving size of the first structure which virializes in the simulation, at the time τ i at which it virializes. Self-similarity and stable clusteringthen imply the temporal evolutions x min x i ∼ e − τ − τ ) / x max x i ∼ R s ( τ ) R s ( τ ) ∼ e − αn +1 Γ( τ − τ ) . (24)We can infer, using Eq. (21), that, the range over whichpower law scaling is expected for τ > τ at the end of thesimulation can be expressed asln (cid:18) x max x min (cid:19) = 11 − γ sc ln L f L i ! = 11 − γ sc ln (cid:18) N f N i (cid:19) (25)where L i and L f are, respectively, the initial comoving sizeof the first and last scale to virialize during the evolution, and N f ( N i ) are the number of particles they contain, re-spectively. Thus, while for γ sc closer to unity this range ofscales is large, as γ sc decreases it contracts and is simply oforder the ratio L f L i for γ sc ≈ γ .For our simulations, with N = 10 particles, taking N i ∼
10, or a little larger, and N f ∼ − , we can, forsmall γ sc , access power law clustering over approximatelytwo orders of magnitude, which leads in practice to the ob-served order of the uncertainly in the exponent ( ∼ . Let us compare our results to those in three dimensions.Numerical simulations have been done by various groups totest the stable clustering hypothesis, via the study of theEdS model starting from power-law initial conditions (see,e.g. Efstathiou et al. (1988); Padmanabhan et al. (1996);Colombi et al. (1996); Jain & Bertschinger (1996, 1998);Ma & Fry (2000b); Smith et al. (2003)). In this case, corre-sponding to κ = 1 in our notation, the predicted exponent,in three dimensions, of the two point correlation function is γ = 3(3 + n ) / (5 + n ). The range of power spectra whichhas been probed numerically is − < n
1. Further, moregenerally, the stable clustering hypothesis leads to scalingrelations for higher order correlation functions, which havebeen tested by some authors (Colombi et al. 1996; Ma & Fry2000b).Until the most recent and extensive study ofSmith et al. (2003), previous works had concluded that thestable clustering prediction worked quite well, with, in somecases, indications of small deviations close to the limits ofnumerical resolution. This last study, on the other hand,reports measurements of the exponent γ which are clearlydiscrepant with the stable clustering prediction: for n = − γ = 0 .
77 rather than γ = 1; for n = − . γ = 0 .
91 ratherthan γ = 1 .
29; for n = − γ = 1 .
26 rather than γ = 1 . n = 0, γ = 1 .
49 rather than γ = 1 . qualitatively different fromthose we have found in our 1D models, for any κ . Whilethe measured γ increases as n increases, just as in one di-mension, the discrepancy with stable clustering in three di-mensions is associated with an exponent γ which is smaller than the predicted stable clustering exponent. Further thereis no indication — in this limited range of n probed — thatthe discrepancy with respect to stable clustering is increas-ing, nor that the associated discrepancy manifests itself asan (exact or approximate) universality of the measured ex-ponent.These differences between 1D and 3D results could, ofcourse, be simply a reflection of the difference of the grav-itational clustering in the two cases. One difference whichmay be essential, as discussed in Joyce & Sicard (2011), isthat in one dimension there are no tidal forces, and a virial-ized structure may only be perturbed by another structureactually crossing it (which is, on the other hand, much morelikely than in three dimensions). However, one of the mainmotivations for this study is precisely that the reliability ofresults from 3D simulations is not clear: the combination ofthe fact that the range of comoving scales probed grows inproportion to N / , and the necessity to introduce a (rel- c (cid:13) , 000–000 xponents of non-linear clustering in one dimension atively large) smoothing parameter to regularize the smallscale divergence in the 3D gravitational force, mean that theexponents (in existing studies) are measured over at mosta little more than one order of magnitude. Further thereare consideable discrepancies between some of the studies.Jain & Bertschinger (1998) find, for example, good agree-ment with stable clustering for the case n = − .Let us focus a little more on these differences between1D and 3D simulations linked to the small scale smoothingof the force. In 1D, as discussed, the equations of motion canbe integrated exactly in absence of smoothing, so that thereis in practice no lower limit on spatial resolution other thanthat arising from the finite particle density, of order the scale x min in the correlation function. We note that, following theanalysis above, we infer that, if stable clustering holds, wehave ln (cid:18) x min x i (cid:19) = γ sc − γ sc ln L f L i ! (26)in one dimension, while the same arguments applied to threedimensions giveln (cid:18) x min x i (cid:19) = γ sc − γ sc ln L f L i ! (27)The scale x i , the comoving size of the first structure whichvirializes, is of order the initial interparticle spacing λ . Inone dimension, even when we use, as here, a PM code, thesmoothing can, without excessive numerical cost, be chosenso that it is at all times considerably smaller than the scale x min , so that we can be absolutely confident that it plays norole in the evolution. In three dimensions, on the contrary,cosmological simulations of a reasonable size (i.e., so thatthe range of scales L f L i which go non-linear cover at least adecade or two) are typically performed with a smoothing ofbetween one tenth and one hundredth of the interparticlespacing λ . Thus 3D simulations cannot in practice simu-late systems manifesting stable clustering only over a verylimited range of length scales, while remaining in the regimewhere ε ≪ x min . While the use of a cut-off ε > x min does notnecessarily imply that the clustering above the scale ε is notaccurately reproduced, it is quite possible that this could bethe case if great care is not taken in the choice of numericalparameters (see Knebe et al. (2000); Joyce & Sylos Labini(2012)). In this respect even when a relatively large smooth-ing is employed, two body relaxation may play a role in dis-rupting the desired collisionless evolution of N -body simula-tions, and such effects can potentially lead to a “pollution”of larger scales than ε (Knebe et al. 2000). In the 1D sys-tem, in contrast, there is in fact no two body collisionality This notable difference is attributed by the latter, on the ba-sis of a qualitative analysis, to the fact that Jain & Bertschinger(1998) evolve to a time at which the amplitude of the power at thescale of the box is higher, leading possibly to poor modelling of thecoupling to the “missing” long wavelength modes. Another no-table difference in the simulations is, however, that a considerablysmaller smoothing length is employed by Jain & Bertschinger(1998). As discussed below, a larger smoothing could be respon-sible for suppressing power on smaller scales, leading potentiallyto an effective exponent which is lowered. analagous to that in three dimensions, and collisional relax-ation (of which the mechanism is not yet fully understood)is extremely slow compared to that induced by 3D two bodyrelaxation(see Joyce & Worrakitpoonpon (2010) and refer-ences therein). Simple estimates show that such relaxationshould play absolutely no role in our 1D simulations, whilethe same is not true in typical 3D simulations. Our studythus suggests that particular attention should be paid tothis point in the analysis of 3D simulations.We note also that simulations of clustering in a 3Duniverse have been performed for the case of a static uni-verse, for the cases n = 0 (Bottaccio et al. 2002) and n = 2(Baertschiger et al. 2007a,b, 2008). The results are in thiscase qualitatively very consistent with those observed hereand in Joyce & Sicard (2011) for the static ( κ = 0) limit.Just as in the 1D simulations self-similarity is observed, andfurther a correlation function which appears to be the inde-pendent of the initial conditions, with a very shallow expo-nent in the inner part, γ ≈ . Finally let us discuss the question of “universality”. A strik-ing feature of our 1D results is that the measured exponentcharacterising clustering appears to become independent ofinitial conditions and cosmology when the stable clusteringapproximation breaks down — more specifically when theexponent n in the initial power spectrum is sufficiently smallfor any given expansion rate parameter κ . In three dimen-sions, following notably the work of Navarro et al. (1996,1997), much emphasis has been placed on a “universality”of cosmological gravitational clustering which manifests it-self as an apparent independence of halo profiles of initialconditions and cosmology. The range over which such a uni-versality might apply is not, however, clear: studies such asthat of Knollmann et al. (2008) show that, even in the EdScosmology, the inner exponent of halos may vary sensiblyfor power law initial conditions as a function of n .Does the apparent universality we observe in 1D cor-respond to a universality in halo profiles? As discussed inJoyce & Sicard (2011), one can identify and extract halosby prescriptions analogous to those used in three dimen-sions. However, the objects thus extracted are simply so“clumpy”, at any scale significantly above the scale x min ,that it is not possible to approximate the run of density bya smooth profile: indeed it can be verified directly that thepower law exponent (over several decades) in the correla-tion function corresponds to a truly scale invariant fractaldistribution of the matter in the corresponding range. If,however, the extracted halos are smoothed (by hand) overa scale much larger than x min , the exponent characterisingthe decay of density about the centre of the halo, should beapproximately the same as that measured in the correlationfunction. This is true because the latter measures, by defi-nition, the mean density as a function of distance about anypoint, which in a simple fractal type distribution is (mod-ulo fluctuations) the same about any point. Thus, in the 1Dmodel, the universality we observe here is indeed expectedto be associated to a universality of inner exponents of halos.Extrapolated to 3D, our results thus suggest that therelative smoothness of halo profiles is an artefact of insuffi- c (cid:13) , 000–000 D. Benhaiem, M. Joyce and F. Sicard cient spatial resolution, and that the minimal scale of sub-structure is limited only by such resolution. In this case, forpower law initial conditions, the observed power law slopeof the correlation function should in fact coincide with theinner exponent of the measured halo profiles. We note thatthis is completely different to what is predicted usually inthree dimensions on the basis of a smooth halo model (seee.g. Ma & Fry (2000b,c,a); Yano & Gouda (2000)). Up tonow 3D numerical studies have not, to our knowledge, ad-dressed this point directly (by comparing measured slopesin halos and two point correlation functions, in scale freecosmologies). We note that the results of Knollmann et al.(2008) on inner halo profiles appear, however, to be veryconsistent with a behaviour completely analogous to whatwe have seen in our 1D models: in their study of power lawinitial conditions in an EdS cosmology, for n varying in arange between − . − .
75, inner slopes for halo pro-files are estimated using different fitting procedures. For thesimplest such procedure it is found (see Fig. 2) that the mea-sured exponent varies from a value in good agreement withstable clustering for the largest value of n , but appears todeviate, as n decreases, towards an exponent which is larger than predicted by stable clustering, tending rapidly towardsan asymptotic exponent a little larger than unity, and con-sistent with that of Navarro et al. (1997). This suggests thatin the 3D EdS model, the universality of halos would trans-late to an insensitivity to initial conditions below an n oforder −
2. Given that the range of effective exponents of ini-tial conditions of typical cosmological (e.g. ΛCDM) initialconditions falls in or close to this range, this would appearto be very consistent with the approximate universality ofhalo profile exponents observed in these kind of simulations.
Our 1D study thus motivates further careful numerical studyof non-linear clustering in the pure matter ( κ = 1) EdSmodel in three dimensions, extending recent studies suchas Smith et al. (2003) and Knollmann et al. (2008). In thesame way as has been done here, the study of this model canbe extended to the full family of EdS models of which onelimit is a static universe. The goal of such a study would beto see if, in this space of models, and for power law initialconditions, a similar structure is found to that we have seenin one dimension. Our 1D study suggests that particular at-tention needs to be applied to controlling for the robustnessof measure of exponents in the two point correlation functionto small scale force smoothing, particularly when an expo-nent below the predicted stable clustering value is found (asin Smith et al. (2003)). Further, study of the relation be-tween exponents measured from the two point correlationfunction and those measured from halos, and their variationas a function of n and κ , should clarify whether the qualita-tive nature of clustering is indeed like that usually assumedin theoretical modelling (smooth halo models) or instead likethat in the 1D models (hierarchical fractal clustering).If the framework suggested by the 1D models — of a“universal” region in the space of initial conditions of clus-tering, delimited by a “critical” value of the predicted sta-ble clustering exponent — is correct, and applies also inthree dimensions, the theoretical problem remains of un-derstanding both the corresponding value of the “universal exponent”, and why this universality applies in the specificregion where it does. Further study of the 1D model maythen help to throw further light on this. As the static modelis in the relevant part of the parameter space, it may sufficeto study this particular limit, i.e., the region in which thereis universality is where the expansion of the universe has lit-tle effect on the non-linear clustering. We note that a veryrecent study in Schulz et al. (2012) shows that there is anapparent universality in the properties of the “halos” formedfrom the 1D collapse of a single structure with a range ofcold initial conditions, and argues that this universality maybe related to that observed numerically in three dimensions.The 1D halos in Schulz et al. (2012), however, are smooth,very different to the non-linear structures we observe here,which are very clumpy even in the static limit of the 1Dmodel. As discussed above, however, in this limit the rangeof non-linear clustering explored by our simulations is stillquite modest, and it may be that larger simulations mightestablish a link to the study of Schulz et al. (2012).With respect to this last point — the effect of varyingparticle number — we note that we have not discussed herethe question of whether the dynamics of these simulations inthe non-linear regime is representative of the fluid or VlasovPoisson limit. In three dimensions this is a crucial questionwhich has been the subject of discussion and some contro-versy (Splinter et al. 1998; Knebe et al. 2000; Power et al.2003; Joyce et al. 2008; Romeo et al. 2008), and one mightexpect that these 1D models, with their large accessible spa-tial resolution, could help to throw light on this question.In N body simulations this question can in principle be di-rectly probed by comparing simulations in which the particlenumber N is varied while keeping the lower cut-off scale tothe initial density fluctuations fixed, or, alternatively, by in-creasing particle density in units of the grid of our PM code.Further, in contrast to three dimensions, it should be feasibleto address this question much more directly by comparisonof the results of the integration of the N body system with adirect numerical integration of the Vlasov-Poisson equationsthemselves. This is another interesting future direction forwork on these models.We are indebted to B. Marcos for numerous discussionsand suggestions. We also thank S. Colombi, B. Miller, J.Morand, J.-L. Rouet, F. Sylos Labini, P. Viot and T. Wor-rakitpoonpon for useful conversations. REFERENCES
Aurell E., Fanelli D., 2002, Astron. Astrophys., 395, 399Baertschiger T., Joyce M., Gabrielli A., Sylos Labini F.,2007a, Phys. Rev., E75, 021113Baertschiger T., Joyce M., Gabrielli A., Sylos Labini F.,2007b, Phys. Rev., E76, 011116Baertschiger T., Joyce M., Sylos Labini F., Marcos B.,2008, Phys. Rev., E77, 051114Binney J., 2004, Mon. Not. R. Astron. Soc., 350, 939Bottaccio M., Pietronero L., Amici A., Miocchi P., Ca-puzzo Dolcetta R., Montuori M., 2002, Physica A, 305,247Colombi S., Bouchet F. R., Hernquist L., 1996, Astrophys.J., 465, 14 c (cid:13) , 000–000 xponents of non-linear clustering in one dimension Efstathiou G., Frenk C. S., White S. D. M., Davis M., 1988,Mon. Not. R. Astron. Soc., 235, 715Ferreira P. G., Joyce M., 1998, Phys. Rev. D., 58, 023503Gabrielli A., Joyce M., 2008, Phys. Rev., E77, 031139Gabrielli A., Joyce M., Sicard F., 2009, Phys. Rev. E, 80,041108Gabrielli A., Sylos Labini F., Joyce M., Pietronero L., 2005,Statistical Physics for Cosmic Structures. SpringerJain B., Bertschinger E., 1996, Astrophys. J., 456, 43Jain B., Bertschinger E., 1998, Astrophys. J., 509, 517Joyce M., Marcos B., Baertschiger T., 2008, Mon. Not. R.Astron. Soc.Joyce M., Sicard F., 2011, M. Not. R. Astron. Soc., 413,1439Joyce M., Sylos Labini F., (to appear) 2012, Mon. Not. R.Astron. Soc.Joyce M., Worrakitpoonpon T., 2010, Journal of StatisticalMechanics: Theory and Experiment, 10, 12Joyce M., Worrakitpoonpon T., 2011, Phys. Rev. E, 84,011139Knebe A., Kravtsov A., Gottl¨ober S., Klypin A., 2000,Mon. Not. Roy. Astron. Soc., 317, 630Knollmann S. R., Power C., Knebe A., 2008, Mon. Not. R.Astron. Soc., 385, 545Ma C.-P., Fry J. N., 2000a, Astrophys. J., 543, 503Ma C.-P., Fry J. N., 2000b, Astrophy. J. Lett., 531, L87Ma C.-P., Fry J. N., 2000c, Astrophys. J. Lett., 538, L107Melott A. L., 1982, Phys. Rev. Lett., 48, 894Melott A. L., 1983, Astrophys. J., 264, 59Miller B., Rouet J., 2002, Phys. Rev., E65, 056121Miller B., Rouet J., 2006, C. R. Phys., 7, 383Miller B. N., Rouet J., 2010a, ArXiv e-printsMiller B. N., Rouet J., 2010b, ArXiv e-printsMiller B. N., Rouet J., Le Guirriec E., 2007, Phys. Rev. E,76, 036705Navarro J. F., Frenk C. S., White S. D. M., 1996, Astro-phys. J., 462, 563Navarro J. F., Frenk C. S., White S. D. M., 1997, Astro-phys. J., 490, 493Noullez A., Aurell E., Fanelli D., 2003, J. Comp. Phys.,186, 697Padmanabhan T., Cen R., Ostriker J. P., Summers F. J.,1996, Astrophys. J., 466, 604Peebles P. J. E., 1980, The Large-Scale Structure of theUniverse. Princeton University PressPower C., Navarro J. F., Jenkins A., Frenk C. S., WhiteS. D. M., Springel V., Stadel J., Quinn T., 2003, Mon.Not. R. Astron. Soc., 338, 14Romeo A. B., Agertz O., Moore B., Stadel J., 2008, Astro-phys. J., 686, 1Schulz A. E., Dehnen W., Jungman G., Tremaine S., 2012,ArXiv e-printsSmith R. E., Peacock J. A., Jenkins A., White S. D. M.,Frenk C. S., Pearce F. R., Thomas P. A., Efstathiou G.,Couchman H. M. P., 2003, Mon. Not. R. Astron. Soc.,341, 1311Splinter R. J., Melott A. L., Shandarin S. F., Suto Y., 1998,Astrophys. J., 497, 38Valageas P., 2006, Phys. Rev.., E74, 016606Yano T., Gouda N., 1998, Astrophy. J. Supp., 118, 267Yano T., Gouda N., 2000, Astrophys. J., 539, 493 c (cid:13)000