Simulations of purely self-gravitating N-body systems are often used in astrophysics and cosmology to study the collisionless limit of such systems. Their results for macroscopic quantities should then converge well for sufficiently large N. Using a study of the evolution from a simple space of spherical initial conditions - including a region characterised by so-called "radial orbit instability" - we illustrate that the values of N at which such convergence is obtained can vary enormously. In the family of initial conditions we study, good convergence can be obtained up to a few dynamical times with N $ \sim 10^3$ - just large enough to suppress two body relaxation - for certain initial conditions, while in other cases such convergence is not attained at this time even in our largest simulations with N $\sim 10^5$. The qualitative difference is due to the stability properties of fluctuations introduced by the N-body discretisation, of which the initial amplitude depends on N. We discuss briefly why the crucial role which such fluctuations can potentially play in the evolution of the N-body system could, in particular, constitute a serious problem in cosmological simulations of dark matter.

During the violent relaxation of a self-gravitating system a significant fraction of its mass may be ejected. If the time varying gravitational field also breaks spherical symmetry this mass can potentially carry angular momentum. Thus starting initial configurations with zero angular momentum can in principle lead to a bound virialized system with non-zero angular momentum. We explore here, using numerical simulations, how much angular momentum can be generated in a virialized structure in this way, starting from configurations of cold particles which are very close to spherically symmetric. For initial configurations in which spherical symmetry is broken only by the Poissonian fluctuations associated with the finite particle number $N$, with $N$ in range $10^3$ to $10^5$, we find that the relaxed structures have standard "spin" parameters $\lambda \sim 10^{-3}$, and decreasing slowly with $N$. For slightly ellipsoidal initial conditions, in which the finite-$N$ fluctuations break the residual reflection symmetries, we observe values $\lambda \sim 10^{-2}$, of the same order of magnitude as those reported for elliptical galaxies. The net angular momentum vector is typically aligned close to normal to the major semi-axis of the triaxial relaxed structure, and also with that of the ejected mass. This simple mechanism may provide an alternative, or complement, to "tidal torque theory" for understanding the origin of angular momentum in astrophysical structures.

We study, using $N$-body simulation, the shape evolution in gravitational collapse of cold uniform spherical system. The central interest is on how the deviation from spherical symmetry depends on particle number $N$. By revisit of the spherical collapse model, we hypothesize that the departure from spherical symmetry is regulated by the finite-$N$ density fluctuation. Following this assumption, the estimate of the flattening of relaxed structures is derived to be $N^{-1/3}$. In numerical part, we find that the virialized states can be characterized by the core-halo structures and the flattenings of the cores fit reasonably well with the prediction. Moreover the results from large $N$ systems suggest the divergence of relaxation time to the final shapes with $N$.We also find that the intrinsic shapes of the cores are considerably diverse as they vary from nearly spherical, prolate, oblate or completely triaxial in each realization. When $N$ increases, this variation is suppressed as the final shapes do not differ much from initial symmetry. In addition, we observe the stable rotation of the virialized states. Further investigation reveals that the origin of this rotation is related in some way to the initial density fluctuation.

We study quasi-stationary states (QSS) resulting from violent relaxation in the one-dimensional self-gravitating "sheet model", revisiting in particular the question of the adequacy of the theory of Lynden-Bell (LB) to describe them. For "waterbag" initial conditions characterized by a single phase space density, the prediction of this theory is, in this model, a function of only one parameter, which can conveniently be chosen to be the ratio of the energy to that in the degenerate limit. Studying a class of such initial conditions in which the shape of the initial waterbag is varied, we find that the LB predictions are reasonably good always in the low energy region, while at higher energies (i.e. in the non-degenerate limit) they are generally not even qualitatively correct, although certain initial conditions can still be found where they are as good as at low energy. We find notably that, in line with what has been observed by Levin et al. in some other models, when LB theory does not work the QSS are always characterized by the presence of a {\it degenerate} core, which these authors explain as the result of dynamical resonances. In short LB theory appears to be a good approximation only when violent relaxation is sufficiently "gentle", and otherwise a degenerate core-halo structure results.

We revisit the issue of relaxation to thermal equilibrium in the so-called "sheet model", i.e., particles in one dimension interacting by attractive forces independent of their separation. We show that this relaxation may be very clearly detected and characterized by following the evolution of order parameters defined by appropriately normalized moments of the phase space distribution which probe its entanglement in space and velocity coordinates. For a class of quasi-stationary states which result from the violent relaxation of rectangular waterbag initial conditions, characterized by their virial ratio R_0, we show that relaxation occurs on a time scale which (i) scales approximately linearly in the particle number N, and (ii) shows also a strong dependence on R_0, with quasi-stationary states from colder initial conditions relaxing much more rapidly. The temporal evolution of the order parameter may be well described by a stretched exponential function. We study finally the correlation of the relaxation times with the amplitude of fluctuations in the relaxing quasi-stationary states, as well as the relation between temporal and ensemble averages.