We explore the formation and relaxation of so-called quasi-stationary states (QSS) for particle distributions in three dimensions interacting via an attractive radial pair potential $V(r \rightarrow \infty) \sim 1/r^\gamma$ with $\gamma > 0$, and either a soft-core or hard-core regularization at small $r$. In the first part of the paper we generalize, for any spatial dimension $d \geq 2$, Chandrasekhar's approach for the case of gravity to obtain analytic estimates of the rate of collisional relaxation due to two body collisions. The resultant relaxation rates indicate an essential qualitative difference depending on the integrability of the pair force at large distances: for $\gamma >d-1$ the rate diverges in the large particle number $N$ (mean field) limit, unless a sufficiently large soft core is present; for $\gamma < d-1$, on the other hand, the rate vanishes in the same limit even in the absence of any regularization. In the second part of the paper we compare our analytical predictions with the results of extensive parallel numerical simulations in $d=3$, for a range of different exponents $\gamma$ and soft cores leading to the formation of QSS. We find, just as for the previously well studied case of gravity (which we also revisit), excellent agreement between the parametric dependence of the observed relaxation times and our analytic predictions. Further, as in the case of gravity, we find that the results indicate that, when large impact factors dominate, the appropriate cut-off is the size of the system (rather than, for example, the mean inter-particle distance). Our results provide strong evidence that the existence of QSS is robust only for long-range interactions with a large distance behavior $\gamma < d-1$; for $\gamma \geq d-1$ the existence of such states will be conditioned strongly on the short range properties of the interaction.

We present a rigorous study of the classical scattering for anytwo-body inter-particle potential of the form $v(r)=g/r^\gamma$, with$\gamma\textgreater{}0$, for repulsive ($g\textgreater{}0$) and attractive ($g\textless{}0$)interactions. We give a derivation of the complete power series of thedeflection angle in terms of the impact factor for the weak scatteringregime (large impact factors) as well as the asymptotic expressionsfor the hard scattering regime (small impact factors). We see a verydifferent qualitative and quantitative behavior depending whether theinteraction is repulsive or attractive. In the latter case, thefamilies of trajectories depend also strongly on the value of$\gamma$. We also study carefully the modifications of the resultswhen a regularization is introduced in the potential at small scales.We check and illustrate all the results with the exact integration ofthe equations of motion.

We study non-linear gravitational clustering from cold gaussian power-law initial conditions in a family of scale-free EdS models, characterized by a free parameter $\kappa$ fixing the ratio between the mass driving the expansion and the mass which clusters. As in the "usual" EdS model, corresponding to $\kappa=1$, self-similarity provides a powerful instrument to delimit the physically relevant clustering resolved by a simulation. Likewise, if stable clustering applies, it implies scale-free non-linear clustering. We derive the corresponding exponent $\gamma_{sc} (n, \kappa)$ of the two point correlation function. We then report the results of extensive N-body simulations, of comparable size to those previously reported in the literature for the case $\kappa=1$, and performed with an appropriate modification of the GADGET2 code. We observe in all cases self-similarity in the two point correlations, down to a lower cut-off which decreases monotonically in time in comoving coordinates. The self-similar part of the non-linear correlation function is fitted well in all cases by a single power-law with an exponent in good agreement with $\gamma_{sc} (n, \kappa)$. Our results thus indicate that stable clustering provides an excellent approximation to the non-linear correlation function over the resolved self-similar scales, at least down to $\gamma_{sc} (n, \kappa) \approx 1$, corresponding to the case $n=-2$ for $\kappa=1$. We conclude, in contrast notably with the results of Smith et al. (2003), that a clear identification of the breakdown of stable clustering in self-similar models - and the possible existence of a "universal" region in which non-linear clustering becomes independent of initial conditions - remains an important open problem, which should be addressed further in significantly larger simulations.

Systems with long-range interactions, such as self-gravitating clusters and magnetically confined plasmas, do not relax to the usual Boltzmann-Gibbs thermodynamic equilibrium, but become trapped in quasi-stationary states (QSS) the life time of which diverges with the number of particles. The QSS are characterized by the lack of ergodicity which can result in a symmetry broken QSS starting from a spherically symmetric particle distribution. We will present a theory which allows us to quantitatively predict the instability threshold for spontaneous symmetry breaking for a class of d-dimensional self-gravitating systems.

Systems with long range interactions present generically the formation of quasi-stationary long-lived non-equilibrium states. These states relax to Boltzmann equilibrium following a dynamics which is not well understood. In this paper we study this process in two-dimensional inhomogeneous self-gravitating systems. Using the Chandrasekhar -- or local -- approximation we write a simple approximate kinetic equation for the relaxation process, obtaining a Fokker -- Planck equation for the velocity distribution with explicit analytical diffusion coefficients. Performing molecular dynamics simulations and comparing them with the evolution predicted by the Fokker -- Planck equation, we observe a good agreement with the model for all the duration of the relaxation, from the formation of the quasi-stationary state to thermal equilibrium. We observe however an overestimate or underestimate of the relaxation rate of the particles with the slower or larger velocities respectively. It is due to systematic errors in estimating the velocities of the particles at the moment of the collisions, inherent to the Chandrasekhar approximation when applied to inhomogeneous systems. Theory and simulations give a scaling of the relaxation time proportional to the number of particles in the system.

We report results of a study of the Newtonian dynamics of N self-gravitating particles which start in a quasi-uniform spherical configuration, without initial velocities. These initial conditions would lead to a density singularity at the origin at a finite time when N \rightarrow \infty, but this singularity is regulated at any finite N (by the associated density fluctuations). While previous studies have focussed on the behaviour as a function of N of the minimal size reached during the contracting phase, we examine in particular the size and energy of the virialized halo which results. We find the unexpected result that the structure decreases in size as N increases, scaling in proportion to N^{-1/3}, a behaviour which is associated with an ejection of kinetic energy during violent relaxation which grows in proportion to N^{1/3}. This latter scaling may be qualitatively understood, and if it represents the asymptotic behaviour in N implies that this ejected energy is unbounded above. We discuss also tests we have performed which indicate that this ejection is a mean-field phenomenon (i.e. a result of collisionless dynamics).

"Quasi-stationary" states are approximately time-independent out of equilibrium states which have been observed in a variety of systems of particles interacting by long-range interactions. We investigate here the conditions of their occurrence for a generic pair interaction V(r \rightarrow \infty) \sim 1/r^a with a > 0, in d>1 dimensions. We generalize analytic calculations known for gravity in d=3 to determine the scaling parametric dependences of their relaxation rates due to two body collisions, and report extensive numerical simulations testing their validity. Our results lead to the conclusion that, for a < d-1, the existence of quasi-stationary states is ensured by the large distance behavior of the interaction alone, while for a > d-1 it is conditioned on the short distance properties of the interaction, requiring the presence of a sufficiently large soft-core in the interaction potential.

We formalize a classification of pair interactions based on the convergence properties of the {\it forces} acting on particles as a function of system size. We do so by considering the behavior of the probability distribution function (PDF) P(F) of the force field F in a particle distribution in the limit that the size of the system is taken to infinity at constant particle density, i.e., in the "usual" thermodynamic limit. For a pair interaction potential V(r) with V(r) \rightarrow \infty) \sim 1/r^a defining a {\it bounded} pair force, we show that P(F) converges continuously to a well-defined and rapidly decreasing PDF if and only if the {\it pair force} is absolutely integrable, i.e., for a > d-1, where d is the spatial dimension. We refer to this case as {\it dynamically short-range}, because the dominant contribution to the force on a typical particle in this limit arises from particles in a finite neighborhood around it. For the {\it dynamically long-range} case, i.e., a \leq d-1, on the other hand, the dominant contribution to the force comes from the mean field due to the bulk, which becomes undefined in this limit. We discuss also how, for a \leq d-1 (and notably, for the case of gravity, a=d-2) P(F) may, in some cases, be defined in a weaker sense. This involves a regularization of the force summation which is generalization of the procedure employed to define gravitational forces in an infinite static homogeneous universe. We explain that the relevant classification in this context is, however, that which divides pair forces with a > d-2 (or a < d-2), for which the PDF of the {\it difference in forces} is defined (or not defined) in the infinite system limit, without any regularization. In the former case dynamics can, as for the (marginal) case of gravity, be defined consistently in an infinite uniform system.

When an open system of classical point particles interacting by Newtonian gravity collapses and relaxes violently, an arbitrary amount of energy may in principle be carried away by particles which escape to infinity. We investigate here, using numerical simulations, how this released energy and other related quantities (notably the binding energy and size of the virialized structure) depends on the initial conditions, for the one parameter family of starting configurations given by randomly distributing N cold particles in a spherical volume. Previous studies have established that the minimal size reached by the system scales approximately as N^{-1/3}, a behaviour which follows trivially when the growth of perturbations (which regularize the singularity of the cold collapse in the infinite N limit) are assumed to be unaffected by the boundaries. Our study shows that the energy ejected grows approximately in proportion to N^{1/3}, while the fraction of the initial mass ejected grows only very slowly with N, approximately logarithmically, in the range of N simulated. We examine in detail the mechanism of this mass and energy ejection, showing explicitly that it arises from the interplay of the growth of perturbations with the finite size of the system. A net lag of particles compared to their uniform spherical collapse trajectories develops first at the boundaries and then propagates into the volume during the collapse. Particles in the outer shells are then ejected as they scatter through the time dependent potential of an already re-expanding central core. Using modified initial configurations we explore the importance of fluctuations at different scales, and discreteness (i.e. non-Vlasov) effects in the dynamics.

The effects of discreteness arising from the use of the N-body method on the accuracy of simulations of cosmological structure formation are not currently well understood. After a discussion of how the relevant discretisation parameters introduced should be extrapolated to recover the Vlasov-Poisson limit, we study numerically, and with analytical methods we have developed recently, the central issue of how finite particle density affects the precision of results. In particular we focus on the power spectrum at wavenumbers around and above the Nyquist wavenumber, in simulations in which the force resolution is taken smaller than the initial interparticle spacing. Using simulations of identical theoretical initial conditions sampled on four different "pre-initial" configurations (three different Bravais lattices, and a glass) we obtain a {\it lower bound} on the real discreteness error. With the guidance of our analytical results, we establish with confidence that the measured dispersion is not contaminated either by finite box size effects or by subtle numerical effects. Our results show notably that, at wavenumbers {\it below} the Nyquist wavenumber, the dispersion increases monotonically in time throughout the simulation, while the same is true above the Nyquist wavenumber once non-linearity sets in. For normalizations typical of cosmological simulations, we find lower bounds on errors at the Nyquist wavenumber of order of a percent, and larger above this scale. The only way this error may be reduced below these levels at these scales, and indeed convergence to the physical limit firmly established, is by extrapolation, at fixed values of the other relevant parameters, to the regime in which the mean comoving interparticle distance becomes less than the force smoothing scale.

Discreteness effects are a source of uncontrolled systematic errors of N-body simulations, which are used to compute the evolution of a self-gravitating fluid. We have already developed the so-called "Particle Linear Theory" (PLT), which describes the evolution of the position of self-gravitating particles located on a perturbed simple cubic lattice. It is the discrete analogue of the well-known (Lagrangian) linear theory of a self-gravitating fluid. Comparing both theories permits to quantify precisely discreteness effects in the linear regime. It is useful to develop the PLT also for other perturbed lattices because they represent different discretizations of the same continuous system. In this paper we detail how to implement the PLT for perturbed cubic Bravais lattices (simple, body and face-centered) in a cubic simulation box. As an application, we will study the discreteness effects -- in the linear regime -- of N-body simulations for which initial conditions have been set-up using these different lattices.

We present the problematic of controlling the discreteness effects in cosmological N-body simulations. We describe a perturbative treatment which gives an approximation describing the evolution under self-gravity of a lattice perturbed from its equilibrium, which allows to trace the evolution of the fully discrete distribution until the time when particles approach one another ("shell-crossing"). Perturbed lattices are typical initial conditions for cosmological N-body simulations and thus we can describe precisely the early time evolution of these simulations. A quantitative comparison with fluid Lagrangian theory permits to study discreteness effects in the linear regime of the simulations. We show finally some work in progress about quantifying discreteness effects in the non-perturbative (highly non-linear) regime of cosmological N-body simulations by evolving different discretizations of the same continuous density field.

In two recent articles a detailed study has been presented of the out of equilibrium dynamics of an infinite system of self-gravitating points initially located on a randomly perturbed lattice. In this article we extend the treatment of the early time phase during which strong non-linear correlations first develop, prior to the onset of ``self-similar'' scaling in the two point correlation function. We establish more directly, using appropriate modifications of the numerical integration, that the development of these correlations can be well described by an approximation of the evolution in two phases: a first perturbative phase in which particles' displacements are small compared to the lattice spacing, and a subsequent phase in which particles interact only with their nearest neighbor. For the range of initial amplitudes considered we show that the first phase can be well approximated as a transformation of the perturbed lattice configuration into a Poisson distribution at the relevant scales. This appears to explain the ``universality'' of the spatial dependence of the asymptotic non-linear clustering observed from both shuffled lattice and Poisson initial conditions.

We apply a recently developed perturbative formalism which describes the evolution under their self-gravity of particles displaced from a perfect lattice to quantify precisely, up to shell crossing, the effects of discreteness in dissipationless cosmological N-body simulations. We give simple expressions, explicitly dependent on the particle density, for the evolution of power in each mode as a function of red-shift. For typical starting red-shifts the effect of finite particle number is to {\it slow down} slightly the growth of power compared to that in the fluid limit (e.g. by about ten percent at half the Nyquist frequency), and to induce also dispersion in the growth as a function of direction at a comparable level. In the limit that the initial red-shift tends to infinity, at fixed particle density, the evolution in fact diverges from that in the fluid limit (described by the Zeldovich approximation). Contrary to widely held belief, this means that a simulation started at a red-shift much higher than the red-shift of shell crossing actually gives a worse, rather than a better, result. We also study how these effects are modified when there is a small-scale regularization of the gravitational force. We show that such a smoothing may reduce the anisotropy of the discreteness effects, but it then {\it increases} their average effect. This behaviour illustrates the fact that the discreteness effects described here are distinct from those usually considered in this context, due to two-body collisions. Indeed the characteristic time for divergence from the collisionless limit is proportional to $N^{2/3}$, rather than $N/ \log N$ in the latter case.

The relation between the results of cosmological N-body simulations, and the continuum theoretical models they simulate, is currently not understood in a way which allows a quantification of N dependent effects. In this first of a series of papers on this issue, we consider the quantification of such effects in the initial conditions of such simulations. A general formalism developed in [1] allows us to write down an exact expression for the power spectrum of the point distributions generated by the standard algorithm for generating such initial conditions. Expanded perturbatively in the amplitude of the input (i.e. theoretical, continuum) power spectrum, we obtain at linear order the input power spectrum, plus two terms which arise from discreteness and contribute at large wavenumbers. For cosmological type power spectra, one obtains as expected, the input spectrum for wavenumbers k smaller than that characteristic of the discreteness. The comparison of real space correlation properties is more subtle because the discreteness corrections are not as strongly localised in real space. For cosmological type spectra the theoretical mass variance in spheres and two point correlation function are well approximated {\it above} a finite distance. For typical initial amplitudes this distance is a few times the inter-particle distance, but it diverges as this amplitude (or, equivalently, the initial red-shift of the cosmological simulation) goes to zero, at fixed particle density. We discuss briefly the physical significance of these discreteness terms in the initial conditions, in particular with respect to the definition of the continuum limit of N-body simulations.

We study the statistics of the force felt by a particle in the class of spatially correlated distribution of identical point-like particles, interacting via a $1/r^2$ pair force (i.e. gravitational or Coulomb), and obtained by randomly perturbing an infinite perfect lattice. In the first part we specify the conditions under which the force on a particle is a well defined stochastic quantity. We then study the small displacements approximation, giving both the limitations of its validity, and, when it is valid, an expression for the force variance. In the second part of the paper we extend to this class of particle distributions the method introduced by Chandrasekhar to study the force probability density function in the homogeneous Poisson particle distribution. In this way we can derive an approximate expression for the probability distribution of the force over the full range of perturbations of the lattice, i.e., from very small (compared to the lattice spacing) to very large where the Poisson limit is recovered. We show in particular the qualitative change in the large-force tail of the force distribution between these two limits. Excellent accuracy of our analytic results is found on detailed comparison with results from numerical simulations. These results provide basic statistical information about the fluctuations of the interactions (i) of the masses in self-gravitating systems like those encountered in the context of cosmological N-body simulations, and (ii) of the charges in the ordered phase of the One Component Plasma.

We present a perturbative treatment of the evolution under their mutual self-gravity of particles displaced off an infinite perfect lattice, both for a static space and for a homogeneously expanding space as in cosmological N-body simulations. The treatment, analogous to that of perturbations to a crystal in solid state physics, can be seen as a discrete (i.e. particle) generalization of the perturbative solution in the Lagrangian formalism of a self-gravitating fluid. Working to linear order, we show explicitly that this fluid evolution is recovered in the limit that the initial perturbations are restricted to modes of wavelength much larger than the lattice spacing. The full spectrum of eigenvalues of the simple cubic lattice contains both oscillatory modes and unstable modes which grow slightly faster than in the fluid limit. A detailed comparison of our perturbative treatment, at linear order, with full numerical simulations is presented, for two very different classes of initial perturbation spectra. We find that the range of validity is similar to that of the perturbative fluid approximation (i.e. up to close to ``shell-crossing''), but that the accuracy in tracing the evolution is superior. The formalism provides a powerful tool to systematically calculate discreteness effects at early times in cosmological N-body simulations.

We investigate the possibility of generating initial conditions for cosmological N-body simulations by simulating a system whose correlations at thermal equilibrium approximate well those of cosmological density perturbations. The system is an appropriately modified version of the standard ``one component plasma'' (OCP). We show first how a well-known semi-analytic method can be used to determine the potential required to produce the desired correlations, and then verify our results for some cosmological type spectra with simulations of the full molecular dynamics. The advantage of the method, compared to the standard one, is that it gives by construction an accurate representation of both the real and reciprocal space correlation properties of the theoretical model. Furthermore the distributions are also statistically homogeneous and isotropic. We discuss briefly the modifications needed to implement the method to produce configurations appropriate for large N-body simulations in cosmology, and also the generation of initial velocities in this context.

We apply a simple linearization, well known in solid state physics, to approximate the evolution at early times of cosmological N-body simulations of gravity. In the limit that the initial perturbations, applied to an infinite perfect lattice, are at wavelengths much greater than the lattice spacing $l$ the evolution is exactly that of a pressureless self-gravitating fluid treated in the analagous (Lagrangian) linearization, with the Zeldovich approximation as a sub-class of asymptotic solutions. Our less restricted approximation allows one to trace the evolution of the discrete distribution until the time when particles approach one another (i.e. ``shell crossing''). We calculate modifications of the fluid evolution, explicitly dependent on $l$ i.e. discreteness effects in the N body simulations. We note that these effects become increasingly important as the initial red-shift is increased at fixed $l$. The possible advantages of using a body centred cubic, rather than simple cubic, lattice are pointed out.

A well known argument in cosmology gives that the power spectrum (or structure function) $P(k)$ of mass density fluctuations produced from a uniform initial state by physics which is causal (i.e. moves matter and momentum only up to a finite scale) has the behaviour $P(k) \propto k^4$ at small $k$. Noting the assumption of analyticity at $k=0$ of $P(k)$ in the standard derivation of this result, we introduce a class of solvable one dimensional models which allows us to study the relation between the behaviour of $P(k)$ at small $k$ and the properties of the probability distribution $f(l)$ for the spatial extent $l$ of mass and momentum conserving fluctuations. We find that the $k^4$ behaviour is obtained in the case that the first {\it six} moments of $f(l)$ are finite. Interestingly the condition that the fluctuations be localised - taken to correspond to the convergence of the first two moments of $f(l)$ - imposes only the weaker constraint $P(k) \propto k^n$ with $n$ anywhere in the range $0< n \leq 4$. We interpret this result to suggest that the causality bound will be loosened in this way if quantum fluctuations are permitted.