{'fraction_non_alphanumeric': 0.034315720126334555, 'fraction_numerical': 0.014202949098025539, 'mean_word_length': 4.54648038688877, 'pattern_counts': {'":': 1, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exponents of non-linear clustering in scale-free one dimensional cosmological simulations

One dimensional versions of cosmological N-body simulations have been shown to share many qualitative behaviours of the three dimensional problem. They can resolve a large range of time and length scales, and admit exact numerical integration. We use such models to study how non-linear clustering depends on initial conditions and cosmology. More specifically, we consider a family of models which, like the 3D EdS model, lead for power-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 study how the corresponding exponent \gamma depends on the initial conditions, characterized by the exponent n of the power spectrum of initial fluctuations, and on a single parameter \kappa controlling the rate of expansion. The space of initial conditions/cosmology divides very clearly into two parts: (1) a region in which \gamma depends strongly on both n and \kappa and where it agrees very well with a simple generalisation of the so-called stable clustering hypothesis in three dimensions, and (2) a region in which \gamma is more or less independent of both the spectrum and the expansion of the universe. We explain the observed location of the boundary in (n, \kappa) space dividing the"stable clustering"region from the"universal"region. We compare and contrast our findings to results in three dimensions, and discuss in particular the light they may throw on the question of"universality"of non-linear clustering in this context.


INTRODUCTION
Cosmological N body simulations are the primary instrument used to make theoretical predictions for structure formation in current models of the universe. Analytical understanding of the non-linear regime, crucial for many nontrivial tests of these models, remains very poor despite the many phenomenological results derived in simulations. Better understanding of the physics of non-linear clustering could help to better constrain and control cosmological simulations, which despite their ever increasing power and sophistication are still subject to questions about their accuracy and resolution. In this work we use simplified one dimensional (1D) "toy models" of the full three dimensional (3D) problem to study non-linear structure formation. More specifically we study, using a family of such models, how the properties of non-linear clustering depend on both initial conditions and cosmology. Our main result is that, for cold initial conditions, we observe that the large phase space of initial conditions and cosmology we explore breaks up into two parts: on the one hand, a region in which there is an apparent "universality" of the non-linear clustering, i.e., in which it depends very weakly both on initial conditions and on the cosmological expansion; and on the other hand, a region where the non-linear clustering depends strongly on both the initial conditions and expansion history. In the latter region, the exponent characterizing the two point correlations is in excellent agreement with a simple analytical prediction which follows from the so-called stable clustering hypothesis, appropriately generalized to this class of models. Further the apparently well localized boundary dividing the phase space into these two disjoint regions can in fact be characterized using the stable clustering prediction, in terms of a critical value of the corresponding exponent. The nonlinear clustering appears to be strictly universal in the limit in which there is no expansion, but also very well approximated for steeper initial spectra when there is expansion. We compare our results with those in three dimensions, and discuss how our results may provide a framework for understanding the much discussed apparent "universality" of halo profiles in this context.
The exploration of 1D models as a tool for understanding the physics of non-linear structure formation in a cosmological context goes back at least as far as the work of Melott (1982Melott ( , 1983 for the case of hot dark matter models. Different initial conditions and variants of the model have been discussed by a number of other authors (Yano & Gouda 1998;Miller & Rouet 2002;Aurell & Fanelli 2002;Miller & Rouet 2006;Valageas 2006;Miller et al. 2007;Gabrielli et al. 2009;Miller & Rouet 2010a). There is also an extensive literature on the statistical mechanics of finite one dimensional self-gravitating systems (see e.g. Joyce & Worrakitpoonpon (2011) and references therein), as well as some studies which use this case to explore 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 that in Joyce & Sicard (2011). This latter paper showed, for the three 1D models studied in the previous literature, the very strong qualitative similarities in the temporal development of clustering to that in the 3D case. For cold cosmological like initial conditions clustering is hierachical and driven by the linear amplification of fluctuations, leading, for the case of power law initial conditions, to self similar evolution of the non-linear correlations. The exponent characterising the non-linear clustering was, further, found to be in very good agreement with that predicted by the so-called "stable clustering" hypothesis. In this paper we extend this study to a one parameter family of models which allows continuous interpolation between the static model and the two expanding models previously considered. This allows us to probe fully the extent of validity of stable clustering, where it breaks down and what happens in this case. The study reveals that there is a fairly abrupt switchover from the validity of stable clustering to an apparent universality in the clustering.
In the next section we define the class of 1D models we study, discussing also their relation to 3D cosmological models. This presentation improves and simplifies that given in Joyce & Sicard (2011). For the family of scale-free models we focus on, we derive, more rigorously than in Joyce & Sicard (2011), the prediction for the behaviour of an isolated structure and use it to derive the exponents predicted for the non-linear two point correlation when the stable clustering hypothesis is made. In the following section we discuss our numerical simulations, and in particular how we calibrate them using "exact" simulations. The next section gives our results, and in the last section we discuss their interpretation and possible relation to the three dimensional problem.

A FAMILY OF 1D SCALE-FREE MODELS
The class of 1D models we study are defined as follows. We take the equations of motion in comoving coordinates of 3D cosmological N body simulations, and simply replace the 3D Newtonian gravitational interaction with the 1D gravitational interaction. This gives 1D equations of motion where xi are the particle positions on the line. The superscript 'J' in the force term indicates that the sum, which extends over an infinite distribution of masses with non-zero mean mass density, is regularized, just as in 3D, by subtraction of the contribution of the mean mass density, i.e., the force term is sourced only by fluctuations about the mean density. As we will discuss further below this sum can be written 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 direct relation between the 1D and 3D models, one can consider the "particles" in the 1D system to represent infinite parallel sheets of zero thickness embedded in 3D. In this case one has the identification g = 2πGΣ where Σ is the surface mass density of the sheets. The mean (comoving) mass density ρ0 of the corresponding 3D universe and the number density n0 of the 1D system are then related as ρ0 = Σn0, and therefore gn0 = 2πGρ0 . (2) For any given 3D cosmological model which provides a scale factor a(t) and associated Hubble constant H(t) we then have a unique corresponding 1D model. The case a(t) = 1 (and H(t) = 0) on the other hand, defines a model analogous to a 3D universe without expansion (which is well defined but does not correspond to a cosmological model derived from general relativity). The one parameter family of models we will study is a continuous interpolation between the 1D model obtained using the 3D EdS cosmology, and the static model. To define them precisely, and motivate their choice, it is convenient to change time variable in Eq. (1) defining τ = dt a 3/2 . This gives the equations of motion in the form i.e., in which all effects of the cosmology appear only as a fluid type damping term, with For a 3D EdS universe with comoving mass density ρ0, i.e., a = (t/t0) 2/3 where t0 = 1/ √ 6πGρ0, we obtain a Γ which is independent of time and given by where we have used Eq.
(2) to derive the last equality. In a ΛCDM cosmology, on the other hand, Γ is approximately constant through the matter dominated era and increases slightly at late times. We consider here simply the set of 1D models in which Γ takes any constant value, and the focus of our study is the dependence of non-linear clustering on this parameter, as well as on the initial conditions. Γ can simply be considered as a control parameter for understanding the role of expansion in the determination of the properties of the non-linear clustering. While most previous studies have considered either the EdS model (also known as the "quintic" model for reasons which will be recalled below), or the static model, one other model in this family, corresponding to the case Γ = √ gn0, introduced in Miller & Rouet (2002), has been 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 obtained by taking the functional dependence of the 3D EdS expansion law, but allowing a freedom in the normalisation of the expansion rate to the matter density, i.e., we can obtain these models taking H 2 = κ 2 8πGρ0/3a 3 where κ is a positive constant. Eq. (4) then still holds, and thus a = e 2Γτ , but instead of (5) we have, using Eq. (2), that Γ = κ gn0/3 .
In terms of 3D cosmological models the case κ > 1, i.e., Γ > Γ EdS , can be considered to be a model in which there is not just "ordinary" matter but an additional pressureless component of the energy density which does not cluster, e.g., a homogeneous scalar field with appropriate equation of state (see e.g. Ferreira & Joyce (1998) and references therein). Equivalently one can consider that the κ = 1 models are obtained using the "correct" normalisation for the Hubble constant, H 2 = 8πGρ0/3, but choosing the 1D coupling freely instead of imposing the relation Eq.
(2). It is important to note that, given that 3D Hubble law is imposed "by hand" on the 1D model, there is in principle no correct value for Γ. Our point of view here is simply to use Γ (or κ) as a control parameter in trying to understand the role of expansion in non-linear structure formation. In other words we are seeking to improve our qualitative understanding only through the study of these toy models.

STABLE CLUSTERING PREDICTIONS
In the context of the problem of structure formation the family of models we study has the interest of being, like EdS models and static models in 3D, scale-free: the expansion introduces no additional time or length scales. This has as a consequence that, for cold particles with initial density fluctuations with a power law spectrum, one expects to obtain, at sufficiently long times, temporal evolution which is "self-similar", i.e., the evolution in time of the spatial correlation functions is equivalent to a simple time dependent rescaling of the spatial variables. The reason why this occurs is that there is, in this case, only a single physically relevant length scale in the problem: the scale at which fluctuations are of order one, of which the evolution is predicted by linear theory 1 . Further if the clustering in the strongly non-linear regime does not depend on this scale one would expect to obtain a scale-free behaviour of spatial correlations in this regime. Such power law behaviour has been observed in both expanding (see, e.g., Efstathiou et al. (1988); Smith et al. (2003) and references therein) and static 3D simulations (Baertschiger et al. 2007a for a range of power law initial conditions. 3D simulations, however, can show such behaviour over a very limited spatial range, making it difficult to establish if it is associated with a true scale invariance of the non-linear clustering. The study of 1D models in the family we are considering (Miller & Rouet 2006;Miller et al. 2007;Miller & Rouet 2010a;Joyce & Sicard 2011), with much greater spatial resolution and numerical accuracy, shows very convincing evidence that these models do indeed give rise to power-law clustering indicative of truly scale-invariant clustering.
A central problem of non-linear structure formation is then that of understanding how the exponent (or possibly exponents) characterizing the non-linear clustering is determined. In 3D EdS cosmology Peebles (see Peebles (1980)) proposed many years ago a simple explanation of how such scale invariant clustering might arise: if highly non-linear structures (of all different sizes) are supposed to evolve essentially independently after their formation, they will virialize and remain stable in physical coordinates. From this one may derive an exponent which depends only on the exponent in the power spectrum of initial density fluctuations. We derive here now, more rigorously than the heuristic derivation given in Joyce & Sicard (2011), the prediction of such a stable clustering hypothesis for the 1D models we study. To do so we simply derive and then analyse the equations of motion of a finite sub-system.

Equations of motion of a finite subsystem
In the simulations which we consider here, we start, as in 3D cosmological simulations, from a particle configuration generated by applying small displacements to an infinite perfect lattice configuration. At any moment in time such a configuration can be fully specified by giving the displacements ui of particles from the sites of a perfect lattice. The position of the i-th particle with respect to some arbitrary origin, is xi = iλ + ui, where λ = n −1 0 . In a particle labelling with xi xi−1 for all i, the displacements are uniquely defined, and such that no two particles cross when they are applied to the lattice. In this case it has been shown (Gabrielli et al. 2009) that the force on the particle i, given by the sum on the right side of Eq. (3), with an appropriate specification of the regularisation, has the simple exact expression where u is the average particle displacement. This result has been shown for an infinite system with statistically translationally invariant displacements in Gabrielli et al. (2009), with the sole assumption of decay of correlations of the displacements at asymptotically large scales. The same result has been shown to hold in Miller & Rouet (2010b) for the case of an infinite periodic system, using an appropriate regularisation of an Ewald summation. In this case, as we will now show explicitly, the expression (7) can be written is the number of particles in the periodic cell to the left (right) of particle i, and x p is the average particle position, i.e., centre of mass of the particles in the cell 2 .
Let us consider now any finite subsystem S of the infinite system defined simply by the N particles, with i = 1...N say, in some interval [x0, x0 + L] at a given instant in time, Substituting ui = xi − iλ in (7), and noting that the number of particles in the interval to the left (right) of the ith particle is given by NL(i) = i − 1 (NR(i) = N − i), we obtain that the force on the particle i may be written (9) where x S is simply the position of the centre of mass of the subsystem, and u S the average displacement of the particles in the subsystem. If the subsystem chosen is the periodic cell of an infinite periodic system, we have u = u S and thus obtain (8). For any generic subsystem S, on the other hand, the last term on the right in (9) is the same for all particles in the subsystem. Defining yi = xi − x S we thus obtain the equation of motion for any particle in the subsystem, relative to the CM of the subsystem, as where the sum extends over the particles in the subsystem. This remains valid for as long as no other particle outside the system crosses the particles at its extremities. Eq. (10) are thus exactly the equations of motion of a subsystem as long as it is isolated in this sense.

Definition of "quasi-physical" coordinates
For Γ = 0 we can transform to the coordinates in which Eq. (10) becomes These coordinates have been chosen in order to remove the damping term in the equations of motion, and all dependence on the expansion now appears only in a time dependent force term, of which the amplitude decays as 1/t 2 . At sufficiently long times, therefore, an isolated subsystem of the infinite 1D "expanding" model has, to an arbitarily good approximation, the equations of motion of a strictly isolated non-expanding 1D self-gravitating system. We will refer to the coordinates in (11) as quasi-physical coordinates: they are the closest analogy we have to what are called physical coordinates in the three dimensional problem (i.e. r = a(t) x, t). In these coordinates the 3D cosmological equations of motion are equivalent to those for an infinite Newtonian self-gravitating system expanding about a centre in an infinite (static) space. Correspondingly, any subsystem (e.g. galaxy) taken far from other mass in an expanding universe will evolve in these coordinates like an isolated Newtonian system (see e.g. Joyce & Sylos Labini (2012) for a derivation). The additional term in 1/t 2 in the 1D model, which we cannot get rid of by a coordinate transformation, appears because we have imposed in the derivation of the model the 3D expansion law rather than using the analogy of this law in 1D (which gives a collapse in a finite time, qualitatively different to that of the 3D expansion which we wish to model). It is important to note that, re-expressed in terms of the three dimensional scale factor a(t), the transformation in (11) from the "comoving" coordinates yi to the "quasiphysical" coordinates in (11) is ri = e 2Γτ /3 yi ≡ a 1/3 yi, and not ri = ayi as one might naively expect 3 .

Scaling of a finite virialized subsystem
The second term on the right hand side of Eq. (10) may be rewritten simply as [ κ 2 9 + 1]gn0yi The first term is of order gN ∼ gnSL where nS is the average comoving density of the subsystem. Thus if we consider any isolated subsystem which is already highly non-linear (with nS ≫ n0) we are necessarily in the regime where the second term may be neglected. In this case the isolated subsystem will thus behave, in quasi-physical coordinates, like a truly isolated 1D self-gravitating system. Just as in 3D, such a system, for large N , in the collisionless limit, attains on a few dynamical time scales (∼ 1/ √ gnS ) a stationary virialized state (see e.g. Joyce & Worrakitpoonpon (2011) and references therein), i.e., a system which becomes time independent in the quasiphysical coordinates ri. More generally, even if not strictly virialized, it rapidly attains a well defined average size about which virial oscillations persist.

Linear amplification and self-similarity
In order to derive the exponent in what follows we will also use the self-similar behaviour of the evolution, which, for the two point correlation function ξ(x, τ ) means that in the spatial range where the self-similarity is observed 4 . The function Rs(τ ) is derived from linear theory, which describes a self-similar evolution from power law initial conditions once the growing mode dominates. In the 1D model 3 We note that Yano & Gouda (1998) give an incorrect generalization of the stable clustering hypothesis to the 1D EdS model assuming that stability will be attained in the coordinates ax. As a result they arrive at the incorrect conclusion from their numerical study that stable clustering is not observed in the model. 4 We follow the standard notation here in which ξ(x) denotes the reduced two point correlation function, i.e., ξ(x) = n −2 0 n(0)n(x) − 1 where n(x) is the local particle density and · · · denotes the average over the statistically translationally invariant point process. See e.g. Peebles (1980) Power law fit ξ(x,t) the linear theory growth of perturbations can be easily derived by using the expression for the force Eq. (7) in terms of the displacements ui. Setting, without loss of generality, the mean displacement u to zero, we have simply which is characterized by the growing mode (15) For κ = 1 we thus recover the same growth law as in the standard 3D EdS model. Assuming a power law power spectrum of density fluctuations P (k) ∝ k n (and thus ξ ∝ 1/x n+1 at large x) the evolution is self-similar with

Exponent of non-linear two point correlation function
Shown in Fig. 1 is a typical two point correlation ξ(x) obtained from the evolution of power law initial conditions in the class of models we study. Between a lower cut-off scale xmin, and an upper cut-off scale xmax a simple power law behaviour of the correlation function is observed ξ ∼ x −γ . The exponent γ may be written as The temporal evolution of the correlation function is found to be self-similar, down to the scale xmin, below which it is broken. The scale xmax is in the range in which selfsimilarity applies, and thus scales in proportion to Rs, with ξ (xmax) a fixed amplitude corresponding to the transition between the linear and non-linear regime.
If we make now the hypothesis that the structures giving rise to the correlation measured at the lower cut-off are stable in quasi-physical coordinates, i.e., these structures behave like isolated approximately virialized structures, we have that xmin ∼ a −1/3 . Further, since for ξ ≫ 1, ξ(x) mesures simply the average density at distance x from an occupied point normalized to the mean density, we have that ξ (xmin) ∼ a 1/3 . It follows that the exponent tends asymptotically (at sufficiently long times) to This result is identical to that derived in Joyce & Sicard (2011) using more heuristic arguments based on the scaling of the energy 5 . It tells us how the exponent, in the stable clustering hypothesis, depends on the initial conditions, characterized by the exponent n of the power spectrum of density fluctuations, and on the expansion rate parametrized by the constant κ, which we recall is Note that as κ → 0 we obtain γsc = 0, i.e. a correlation function which is predicted to be flat over an arbitrarily large distance. Indeed, in this static limit, the stable clustering hypothesis corresponds to stability also in comoving coordinates, and therefore to a constant xmin and ξ (xmin). Clearly, however, we do not expect in this case that the formation of non-linear structure at a given scale will leave approximately unperturbed the structures formed previously at smaller scales.

Method and calibration
The numerical study in Joyce & Sicard (2011), as most previous studies of the three models belonging to the class we are considering, exploits fully the very attractive property of these models that they admit "exact" numerical integration: the force on particles, whether given in the form (7) or (8), is such that the motion (with or without damping) admits an analytical solution other than at particle crossings. An event driven algorithm can then be employed, in which the determination of the time of the next "collision", and which pair of particles it involves, requires only the solution of 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, 1, √ 3, the algebraic equations involved are, respectively, quadratic, quintic and cubic. For a generic κ, however, the equations are not polynomial, and the numerical cost of solving them increases and makes the code, which is already expensive numerically, even more time consuming 6 . As we wish 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 enhanced by the use of FFT techniques and parallelization 7 . In order to be sure that we are not, as a result, losing the advantages of the precision of the integration accessible in these 1D models, we calibrate, as discussed below, our choice of the new discreteness parameters introduced by the PM code by comparing our results with those obtained, for κ = 0, 1, √ 3, using the "exact" code.
Our PM code defines a regular grid on the simulation box, on which a density field n(x) is defined at each time step by a "clouds in cell" interpolation of the particle positions. The modified 1D Poisson equation is then solved for the potential in Fourier space using an FFT, using the solver fftw3 8 . An inverse FFT then determines the force on the grid, which is then interpolated back onto the particle positions. The latter are then advanced using a leapfrog algorithm, with a simple adaptative timestep: it is chosen at each time by imposing that the mean distance travelled by particles during a time step be equal to a chosen fraction η of the PM grid. Compared to the exact code, our PM code therefore introduces resolution effects controlled by two dimensionless parameters: the size of the PM grid compared to the lattice spacing in the initial conditions, which we denote by ǫ, and the parameter η controlling the size of the time steps. We have made an extensive study of the effects introduced by this finite resolution of the PM code, and in particular on the determination of the two point correlation function which is the quantity which interests us here. Shown in Fig. 2 is, for example, a comparison of the results for runs from identical initial conditions, for the case n = 2 and κ = 1, obtained from the exact code and the PM code, for two different values of ǫ, and η sufficiently small (of order unity) so that convergence is obtained with respect to it. We find, very reasonably, that excellent agreement is obtained provided ǫ is chosen smaller than xmin, while a larger ǫ than this minimal scale leads to a visible deviation of the correlation function from its correct value at scales below ǫ. We have chosen our numerical parameters in the simulations reported below in order to obtain such an indiscernible difference in the correlation function between our code and the exact code, for the cases κ = 0, 1, √ 3. All the results presented here are for systems with N = 10 5 particles (in periodic boundary conditions).
For initial conditions with power law power spectra P (k) ∼ k n (at small k, a cut-off at large k is always implicitly assumed), we expect hierarchical structure formation in d dimensions for n in the range −d < n 4 (Peebles 1980): for n −d the variance of mass fluctuations diverges 6 The equation for the crossing times is of the form 1 + Az β + Bz (β+1)/2 = 0 where z = e Γτ , β = 1 + 24 κ 2 , and A and B are constants. 7 See, e.g., http://astro.uchicago.edu/∼andrey/Talks/PM/pm.pdf and references therein. x/L exact simulation ξ(x,t) for ε=0.76 ξ(x,t) for ε=0.012 Figure 2. Two point correlation functions obtained by evolving from the same initial condition (with n = 2) and κ = 1, for three different cases: using the "exact" code (solid line) and using our PM code with the two indicated values of the PM grid spacing.
at large scales, while for n > 4 linear theory breaks down (and the dynamics at large scales will be determined by the power at the ultra-violet scale). We thus consider (as in Joyce & Sicard (2011)) the cases n = 0, 2, 4. Point processes with such power spectra are easy to set up: n = 0 is obtained by randomly distributing points in the interval (giving a power spectrum P (k) = 1/n0), the case n = 2 by applying small random and uncorrelated displacements to a regular lattice, and n = 4 in the same way but with the additional constraint that neighbouring pairs have equal and opposite displacements (see Gabrielli & Joyce (2008) and references therein).

Results
Our results are summarized in Fig. 3. Each point on the plot corresponds to a single simulation with given n and κ, where the latter varies in the range 0 -2.5. The figure shows then the value of the exponent γ of the power law region of the correlation function. In each case the given value is determined using a linear regression in the range between xmin and xmax, the latter being determined themselves by eye (cf. Fig. 1). The uncertainty in the fitted exponent associated with the (by eye) estimation of the end-points for the fit is very small (of order a few percent at most) in most parts of the parameter space. However, as we will now discuss and explain, at lower values of κ and n, where the extent of the power law region is much smaller, the resultant uncertainty becomes more significant. The continuous curves plotted correspond to the predictions of the stable clustering hypothesis, Eq. (18).
Our results show excellent agreement with the stable clustering prediction for a very large part of the explored parameter space of initial conditions and cosmology. This region can be well characterized simply by the condition γsc(n, κ) ∼ > 0.2, i.e., good agreement with stable clustering is observed simply in the region where the predicted exponent be larger than some critical value. In the rest of the parameter space (i.e. for γsc(n, κ) ∼ < 0.2) there is clear disagreement with the stable clustering prediction, and the measured exponents lie in a very narrow range, between 0.2 and 0.15, or a little smaller. It is precisely in this part of the parameter space, however, that there is also a considerable scatter about the curves of the measured exponents. This scatter is in fact just a measure of the uncertainty in the determination of these exponents, which is most difficult (for reasons we explain in detail below) in the region where the exponents become small. Thus, our results are quite consistent with the hypothesis that the exponent in the region γsc(n, κ) ∼ < 0.2 is "universal", in the sense that it is independent of initial conditions and cosmology. A schematic representation of this result is shown in Fig. 4: the parameter space of initial conditions and cosmology breaks up into a "stable clustering" region, and a region of "universality". The contour defining the two disjoint regions is simply given by γsc(n, κ) = γ0. Rather than a line, there may of course be a small region in parameter space in which the exponent varies in a continuous manner from its stable clustering value to a limiting universal value of γ0 ≈ 0.15 9 . 9 In the language of statistical physics our result could be described as an evidence for an "out of equilibrium phase transition": as control parameters for the system are changed the (out of equilibrium) behaviour of the system changes in a discontinuous manner. With the precision of our current results we cannot determine with confidence whether there is really such a transition or rather a continuous "crossover" between the two "phases".

Interpretation of 1D results
We have found numerically that the stable clustering hypothesis predicts the exponents of non-linear clustering very accurately at sufficiently large n and κ, i.e., if the spectrum of initial fluctuations is "sufficiently blue" and the expansion rate is sufficiently fast. More precisely it appears that the criterion for it to work is that γsc, which is a monotonically increasing function of both n and κ, be larger than some critical value. It is not difficult to give a simple physical explanation for this behaviour, as we now explain. Let us consider two overdense regions, initially (in the linear regime) of comoving size L 0 1 and L 0 2 > L 0 1 . Selfsimilarity of their evolution implies that there is simply a time delay ∆τ12 between every stage of their evolution given by Rs(∆τ12) = (L 0 2 /L 0 1 ). Using Eq. (16) we obtain ∆τ12 = 1 + n 4Γα ln(L 0 2 /L 0 1 ) .
Let us suppose that the fluctuation 1 goes non-linear and approximately virializes at some given time τ v 1 , when its size is L v 1 . When fluctuation 2 virializes, at time τ v 2 = τ v 1 + ∆τ12, its size is L v 2 = (L 0 2 /L 0 1 )L v 1 . If in the interval between τ v 1 and τ v 2 , structure 1 behaves to a very good approximation as an isolated structure, its size in comoving coordinates will decrease by a factor e −2Γ∆τ /3 . It follows that the relative size of the two structures at the time when 2 virializes is where α is the growing mode exponent given in (15). Quite simply, the faster the expansion rate, or the larger is n, the more a given scale which has entered the non-linear regime can "shrink" relative to a larger scale in the time before the latter goes non-linear. Further we note, using (15), that Eq. (22) can be written i.e., the dependence on (n, κ) is in fact completely determined through the exponent γsc itself.
The predicted stable clustering exponent is therefore directly related to a physical quantity which we would expect naturally to control the validity of the assumption made in deriving it: the more "concentrated" are the pre-existing virialized substructures inside a larger structure when it collapses, the better should become the approximation that this structure will behave simply as a collection of sub-structures which are not disrupted by the virialization and subsequent evolution of the larger structure. Indeed, in the limit that γsc approaches unity, any structure which collapses and virializes will see the substructures it contains essentially as point particles. As γsc decreases, on the other hand, we expect that the interaction between structures can lead to their disruption, and in particular that "mergers" of substructures become much more probable. Given that the value of the exponent γsc characterizes precisely the relative "condensation" of scales, it is very natural that a critical value of the exponent should characterize the breakdown of stable clustering.
What our results do not allow us to conclude, as we have discussed, is what precisely happens in the regime where stable clustering breaks down: our numerical results do not allow us to clearly distinguish between the possibility of an abrupt (discontinuous) transition to a region in which there is a truly universal exponent, or a smoother transition to a region in which the exponent depends only very weakly on the relevant parameters (n, κ). Further we do not (currently) have a model to explain the exponent (or narrow range of exponents) observed in this part of the parameter space.
Extending the analysis given just above to determine the scalings of relative sizes of structures, it is simple to understand why it is precisely the region of small γsc in which the numerical results for the exponent of the power law in the correlation function are most noisy. The degree of precision in the measured exponent is essentially just a function of the range of scale over which the power-law behaviour extends, i.e., it depends on the ratio xmin/xmax accessible in a numerical simulation. The scale xmax corresponds approximately at any time to the size of the largest (approximately) virialized scale. Let use denote by xi the comoving size of the first structure which virializes in the simulation, at the time τi at which it virializes. Self-similarity and stable clustering then imply the temporal evolutions (24) We can infer, using Eq. (21), that, the range over which power law scaling is expected for τ > τ0 at the end of the simulation can be expressed as where L 0 i and L 0 f are, respectively, the initial comoving size of the first and last scale to virialize during the evolution, and N f (Ni) are the number of particles they contain, respectively. Thus, while for γsc closer to unity this range of scales is large, as γsc decreases it contracts and is simply of For our simulations, with N = 10 5 particles, taking Ni ∼ 10, or a little larger, and N f ∼ 10 3 − 10 4 , we can, for small γsc, access power law clustering over approximately two orders of magnitude, which leads in practice to the observed order of the uncertainly in the exponent (∼ 0.05).

Comparison with 3D simulations: correlation analysis
Let us compare our results to those in three dimensions. Numerical simulations have been done by various groups to test the stable clustering hypothesis, via the study of the EdS 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; Ma & Fry (2000b); Smith et al. (2003)). In this case, corresponding to κ = 1 in our notation, the predicted exponent, in three dimensions, of the two point correlation function is γ3 = 3(3 + n)/(5 + n). The range of power spectra which has been probed numerically is −3 < n 1. Further, more generally, the stable clustering hypothesis leads to scaling relations for higher order correlation functions, which have been tested by some authors (Colombi et al. 1996;Ma & Fry 2000b). Until the most recent and extensive study of Smith et al. (2003), previous works had concluded that the stable clustering prediction worked quite well, with, in some cases, indications of small deviations close to the limits of numerical resolution. This last study, on the other hand, reports measurements of the exponent γ which are clearly discrepant with the stable clustering prediction: for n = −2, γ = 0.77 rather than γ3 = 1; for n = −1.5, γ = 0.91 rather than γ3 = 1.29; for n = −1, γ = 1.26 rather than γ3 = 1.5; for n = 0, γ = 1.49 rather than γ3 = 1.8. These results are clearly qualitatively different from those we have found in our 1D models, for any κ. While the measured γ increases as n increases, just as in one dimension, the discrepancy with stable clustering in three dimensions is associated with an exponent γ which is smaller than the predicted stable clustering exponent. Further there is no indication -in this limited range of n probed -that the discrepancy with respect to stable clustering is increasing, nor that the associated discrepancy manifests itself as an (exact or approximate) universality of the measured exponent.
These differences between 1D and 3D results could, of course, be simply a reflection of the difference of the gravitational clustering in the two cases. One difference which may be essential, as discussed in Joyce & Sicard (2011), is that in one dimension there are no tidal forces, and a virialized structure may only be perturbed by another structure actually crossing it (which is, on the other hand, much more likely than in three dimensions). However, one of the main motivations for this study is precisely that the reliability of results from 3D simulations is not clear: the combination of the fact that the range of comoving scales probed grows in proportion to N 1/3 , and the necessity to introduce a (rel-atively large) smoothing parameter to regularize the small scale divergence in the 3D gravitational force, mean that the exponents (in existing studies) are measured over at most a little more than one order of magnitude. Further there are consideable discrepancies between some of the studies. Jain & Bertschinger (1998) find, for example, good agreement with stable clustering for the case n = −2 with simulations of a similar size to those of Smith et al. (2003) 10 .
Let us focus a little more on these differences between 1D and 3D simulations linked to the small scale smoothing of the force. In 1D, as discussed, the equations of motion can be integrated exactly in absence of smoothing, so that there is in practice no lower limit on spatial resolution other than that arising from the finite particle density, of order the scale xmin in the correlation function. We note that, following the analysis above, we infer that, if stable clustering holds, we have in one dimension, while the same arguments applied to three dimensions give The scale xi, the comoving size of the first structure which virializes, is of order the initial interparticle spacing λ. In one dimension, even when we use, as here, a PM code, the smoothing can, without excessive numerical cost, be chosen so that it is at all times considerably smaller than the scale xmin, so that we can be absolutely confident that it plays no role in the evolution. In three dimensions, on the contrary, cosmological simulations of a reasonable size (i.e., so that the range of scales L 0 f L 0 i which go non-linear cover at least a decade or two) are typically performed with a smoothing of between one tenth and one hundredth of the interparticle spacing λ. Thus 3D simulations cannot in practice simulate systems manifesting stable clustering only over a very limited range of length scales, while remaining in the regime where ε ≪ xmin. While the use of a cut-off ε xmin does not necessarily imply that the clustering above the scale ε is not accurately reproduced, it is quite possible that this could be the case if great care is not taken in the choice of numerical parameters (see Knebe et al. (2000); Joyce & Sylos Labini (2012)). In this respect even when a relatively large smoothing is employed, two body relaxation may play a role in disrupting the desired collisionless evolution of N -body simulations, and such effects can potentially lead to a "pollution" of larger scales than ε (Knebe et al. 2000). In the 1D system, in contrast, there is in fact no two body collisionality analagous to that in three dimensions, and collisional relaxation (of which the mechanism is not yet fully understood) is extremely slow compared to that induced by 3D two body relaxation(see Joyce & Worrakitpoonpon (2010) and references therein). Simple estimates show that such relaxation should play absolutely no role in our 1D simulations, while the same is not true in typical 3D simulations. Our study thus suggests that particular attention should be paid to this point in the analysis of 3D simulations. We note also that simulations of clustering in a 3D universe have been performed for the case of a static universe, for the cases n = 0 (Bottaccio et al. 2002) and n = 2 (Baertschiger et al. 2007a. The results are in this case qualitatively very consistent with those observed here and in Joyce & Sicard (2011) for the static (κ = 0) limit. Just as in the 1D simulations self-similarity is observed, and further a correlation function which appears to be the independent of the initial conditions, with a very shallow exponent in the inner part, γ ≈ 0.3.

Comparison with 3D simulations: halos and universality
Finally let us discuss the question of "universality". A striking feature of our 1D results is that the measured exponent characterising clustering appears to become independent of initial conditions and cosmology when the stable clustering approximation breaks down -more specifically when the exponent n in the initial power spectrum is sufficiently small for any given expansion rate parameter κ. In three dimensions, following notably the work of Navarro et al. (1996Navarro et al. ( , 1997, much emphasis has been placed on a "universality" of cosmological gravitational clustering which manifests itself as an apparent independence of halo profiles of initial conditions and cosmology. The range over which such a universality might apply is not, however, clear: studies such as that of Knollmann et al. (2008) show that, even in the EdS cosmology, the inner exponent of halos may vary sensibly for power law initial conditions as a function of n. Does the apparent universality we observe in 1D correspond to a universality in halo profiles? As discussed in Joyce & Sicard (2011), one can identify and extract halos by prescriptions analogous to those used in three dimensions. However, the objects thus extracted are simply so "clumpy", at any scale significantly above the scale xmin, that it is not possible to approximate the run of density by a smooth profile: indeed it can be verified directly that the power law exponent (over several decades) in the correlation function corresponds to a truly scale invariant fractal distribution of the matter in the corresponding range. If, however, the extracted halos are smoothed (by hand) over a scale much larger than xmin, the exponent characterising the decay of density about the centre of the halo, should be approximately the same as that measured in the correlation function. This is true because the latter measures, by definition, the mean density as a function of distance about any point, which in a simple fractal type distribution is (modulo fluctuations) the same about any point. Thus, in the 1D model, the universality we observe here is indeed expected to be associated to a universality of inner exponents of halos.
Extrapolated to 3D, our results thus suggest that the relative smoothness of halo profiles is an artefact of insuffi-cient spatial resolution, and that the minimal scale of substructure is limited only by such resolution. In this case, for power law initial conditions, the observed power law slope of the correlation function should in fact coincide with the inner exponent of the measured halo profiles. We note that this is completely different to what is predicted usually in three dimensions on the basis of a smooth halo model (see e.g. Ma & Fry (2000b,c,a); Yano & Gouda (2000)). Up to now 3D numerical studies have not, to our knowledge, addressed this point directly (by comparing measured slopes in halos and two point correlation functions, in scale free cosmologies). We note that the results of Knollmann et al. (2008) on inner halo profiles appear, however, to be very consistent with a behaviour completely analogous to what we have seen in our 1D models: in their study of power law initial conditions in an EdS cosmology, for n varying in a range between −0.5 and −2.75, inner slopes for halo profiles are estimated using different fitting procedures. For the simplest such procedure it is found (see Fig. 2) that the measured exponent varies from a value in good agreement with stable clustering for the largest value of n, but appears to deviate, as n decreases, towards an exponent which is larger than predicted by stable clustering, tending rapidly towards an asymptotic exponent a little larger than unity, and consistent with that of Navarro et al. (1997). This suggests that in the 3D EdS model, the universality of halos would translate to an insensitivity to initial conditions below an n of order −2. Given that the range of effective exponents of initial conditions of typical cosmological (e.g. ΛCDM) initial conditions falls in or close to this range, this would appear to be very consistent with the approximate universality of halo profile exponents observed in these kind of simulations.

Future studies
Our 1D study thus motivates further careful numerical study of non-linear clustering in the pure matter (κ = 1) EdS model in three dimensions, extending recent studies such as Smith et al. (2003) and Knollmann et al. (2008). In the same way as has been done here, the study of this model can be extended to the full family of EdS models of which one limit is a static universe. The goal of such a study would be to see if, in this space of models, and for power law initial conditions, a similar structure is found to that we have seen in one dimension. Our 1D study suggests that particular attention needs to be applied to controlling for the robustness of measure of exponents in the two point correlation function to small scale force smoothing, particularly when an exponent below the predicted stable clustering value is found (as in Smith et al. (2003)). Further, study of the relation between exponents measured from the two point correlation function and those measured from halos, and their variation as a function of n and κ, should clarify whether the qualitative nature of clustering is indeed like that usually assumed in theoretical modelling (smooth halo models) or instead like that 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 clustering, delimited by a "critical" value of the predicted stable clustering exponent -is correct, and applies also in three dimensions, the theoretical problem remains of understanding both the corresponding value of the "universal exponent", and why this universality applies in the specific region where it does. Further study of the 1D model may then help to throw further light on this. As the static model is in the relevant part of the parameter space, it may suffice to study this particular limit, i.e., the region in which there is universality is where the expansion of the universe has little effect on the non-linear clustering. We note that a very recent study in Schulz et al. (2012) shows that there is an apparent universality in the properties of the "halos" formed from the 1D collapse of a single structure with a range of cold initial conditions, and argues that this universality may be 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 1D model. As discussed above, however, in this limit the range of non-linear clustering explored by our simulations is still quite modest, and it may be that larger simulations might establish a link to the study of Schulz et al. (2012).
With respect to this last point -the effect of varying particle number -we note that we have not discussed here the question of whether the dynamics of these simulations in the non-linear regime is representative of the fluid or Vlasov Poisson limit. In three dimensions this is a crucial question which has been the subject of discussion and some controversy (Splinter et al. 1998;Knebe et al. 2000;Power et al. 2003;Joyce et al. 2008;Romeo et al. 2008), and one might expect that these 1D models, with their large accessible spatial resolution, could help to throw light on this question. In N body simulations this question can in principle be directly probed by comparing simulations in which the particle number N is varied while keeping the lower cut-off scale to the initial density fluctuations fixed, or, alternatively, by increasing particle density in units of the grid of our PM code. Further, in contrast to three dimensions, it should be feasible to address this question much more directly by comparison of the results of the integration of the N body system with a direct numerical integration of the Vlasov-Poisson equations themselves. This is another interesting future direction for work on these models.
We are indebted to B. Marcos for numerous discussions and suggestions. We also thank S. Colombi, B. Miller, J. Morand, J.-L. Rouet, F. Sylos Labini, P. Viot and T. Worrakitpoonpon for useful conversations.