The evolution of self-gravitating systems, and long-range interacting systems more generally, from initial configurations far from dynamical equilibrium is often described as a simple two phase process: a first phase of violent relaxation bringing it to a quasi-stationary state in a few dynamical times, followed by a slow adiabatic evolution driven by collisional processes. In this context the complex spatial structure evident, e.g., in spiral galaxies is understood either in terms of instabilities of quasi-stationary states, or a result of dissipative non-gravitational interactions. We illustrate here, using numerical simulations, that purely self-gravitating systems evolving from quite simple initial configurations can in fact give rise easily to structures of this kind of which the lifetime can be large compared to the dynamical characteristic time, but short compared to the collisional relaxation time scale. More specifically, for a broad range of non-spherical and non-uniform rotating initial conditions, gravitational relaxation gives rise quite generically to long-lived non-stationary structures of a rich variety, characterized by spiral-like arms, bars and even ring-like structures in special cases. These structures are a feature of the intrinsically out-of-equilibrium nature of the system's collapse, associated with a part of the system's mass while the bulk is well virialized. They are characterized by predominantly radial motions in their outermost parts, but also incorporate an extended flattened region which rotates coherently about a well virialized core of triaxial shape with an approximately isotropic velocity dispersion. We discuss the possible relevance of these simple toy models to the observed structure of real galaxies emphasizing the difference between dissipative and dissipationless disc formation.

Galaxy velocity maps often show the typical pattern of a rotating disk, consistent with the dynamical model where emitters rotate in circular orbits around the galactic center. The simplest template used to fit these maps consists in the rotating disk model (RDM) where the amplitude of circular velocities is fixed by the observed velocity profile along the kinematic axis. A more sophisticated template is the rotating tilted-ring model (RTRM) that takes into account the presence of warps and allows a radius-dependent orientation of the kinematic axis. In both cases, axisymmetry is assumed and residuals between the observed and the model velocity fields are interpreted as noncircular motions. We show that if a galaxy is not axisymmetric, there is an intrinsic degeneracy between a rotational and a radial velocity field. We then introduce a new galaxy template, the radial ellipse model (REM), that is not axisymmetric and has a purely radial velocity field with an amplitude that is correlated with the major axis of the ellipse. We show that best fits to the observed two-dimensional velocity fields of 28 galaxies extracted from the THINGS sample with both the REM and the RDM give residuals with similar amplitudes, where the REM residuals trace nonradial motions. Best fits obtained with the RTRM, because of its larger number of free parameters, give the smallest residuals: however, we argue that this does not necessarily imply that the RTRM gives the most accurate representation of a galaxy velocity field. Instead, we show that this method is not able to disentangle between circular and radial motions for the case of nonaxisymmetric systems. Finally, we consider the physical motivation of the REM, and discuss how the interpretation of galactic dynamics changes if one assumes that the main component of a galaxy velocity field is modeled as a RDM/RTRM or as a REM.

We describe how a simple class of out of equilibrium, rotating and asymmetrical mass distributions evolve under their self-gravity to produce a quasi-planar spiral structure surrounding a virialized core, qualitatively resembling a spiral galaxy. The spiral structure is transient, but can survive tens of dynamical times, and further reproduces qualitatively noted features of spiral galaxies as the predominance of trailing two-armed spirals and large pitch angles. As our models are highly idealized, a detailed comparison with observations is not appropriate, but generic features of the velocity distributions can be identified to be potential observational signatures of such a mechanism. Indeed, the mechanism leads generically to a characteristic transition from predominantly rotational motion, in a region outside the core, to radial ballistic motion in the outermost parts. Such radial motions are excluded in our Galaxy up to 15 kpc, but could be detected at larger scales in the future by GAIA. We explore the apparent motions seen by external observers of the velocity distributions of our toy galaxies, and find that it is difficult to distinguish them from those of a rotating disc with sub-dominant radial motions at levels typically inferred from observations. These simple models illustrate the possibility that the observed apparent motions of spiral galaxies might be explained by non-trivial non-stationary mass and velocity distributions without invoking a dark matter halo or modification of Newtonian gravity. In this scenario the observed phenomenological relation between the centripetal and gravitational acceleration of the visible baryonic mass could have a simple explanation.

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.

The determination of the resolution of cosmological N-body simulations, i.e., the range of scales in which quantities measured in them represent accurately the continuum limit, is an important open question. We address it here using scale-free models, for which self-similarity provides a powerful tool to control resolution. Such models also provide a robust testing ground for the so-called stable clustering approximation, which gives simple predictions for them. Studying large N-body simulations of such models with different force smoothing, we find that these two issues are in fact very closely related: our conclusion is that the accuracy of two point statistics in the non-linear regime starts to degrade strongly around the scale at which their behaviour deviates from that predicted by the stable clustering hypothesis. Physically the association of the two scales is in fact simple to understand: stable clustering fails to be a good approximation when there are strong interactions of structures (in particular merging) and it is precisely such non-linear processes which are sensitive to fluctuations at the smaller scales affected by discretisation. Resolution may be further degraded if the short distance gravitational smoothing scale is larger than the scale to which stable clustering can propagate. We examine in detail the very different conclusions of studies by Smith et al. (2003) and Widrow et al. (2009} and find that the strong deviations from stable clustering reported by these works are the results of over-optimistic assumptions about scales resolved accurately by the measured power spectra, and the reliance on Fourier space analysis. We emphasise the much poorer resolution obtained with the power spectrum compared to the two point correlation function..

We study the collapse of an isolated, initially cold, irregular (but almost spherical) and (slightly) inhomogeneous cloud of self-gravitating particles. The cloud is driven towards a virialized quasi-equilibrium state by a fast relaxation mechanism, occurring in a typical time $\tau_c$, whose signature is a large change in the particle energy distribution. Post-collapse particles are divided into two main species: bound and free, the latter being ejected from the system. Because of the initial system's anisotropy, the time varying gravitational field breaks spherical symmetry so that the ejected mass can carry away angular momentum and the bound system can gain a non-zero angular momentum. In addition, while strongly bound particles form a compact core, weakly bound ones may form, in a time scale of the order of $\tau_c$, several satellite sub-structures. These satellites have a finite lifetime that can be longer than $\tau_c$ and generally form a flattened distribution. Their origin and their abundance are related to the amplitude and nature of initial density fluctuations and to the initial cloud deviations from spherical symmetry, which are both amplified during the collapse phase. Satellites show a time dependent virial ratio that can be different from the equilibrium value $b\approx -1$: although they are bound to the main virialized object, they are not necessarily virially relaxed.

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.

An isolated, initially cold and ellipsoidal cloud of self-gravitating particles represents a relatively simple system to study the effects of the deviations from spherical symmetry in the mechanism of violent relaxation. Initial deviations from spherical symmetry are shown to play a dynamical role that is equivalent to that of density fluctuations in the case of an initially spherical cloud. Indeed, these deviations control the amount of particles energy change and thus determine the properties of the final energy distribution, particularly the appearance of two species of particles: bound and free. Ejection of mass and energy from the system together with the formation of a density profile decaying as $\rho(r) \sim r^{-4}$ and a Keplerian radial velocity dispersion profile, are the prominent features similar to those observed after the violent relaxation of spherical clouds. In addition, we find that ejected particles are characterized by highly non-spherical shapes, whose features can be traced in the initial deviations from spherical symmetry that are amplified during the dynamical evolution: particles can indeed form anisotropic configurations, like bars and/or disks, even though the initial cloud was very close to spherical.

Initially cold and spherically symmetric self-gravitating systems may give rise to a virial equilibrium state which is far from spherically symmetric, and typically triaxial. We focus here on how the degree of symmetry breaking in the final state depends on the initial density profile. We note that the most asymmetric structures result when, during the collapse phase, there is a strong injection of energy preferentially into the particles which are localized initially in the outer shells. These particles are still collapsing when the others, initially located in the inner part, are already re-expanding; the motion of particles in a time varying potential allow them to gain kinetic energy --- in some cases enough to be ejected from the system. We show that this mechanism of energy gain amplifies the initial small deviations from perfect spherical symmetry due to finite $N$ fluctuations. This amplification is more efficient when the initial density profile depends on radius, because particles have a greater spread of fall times compared to a uniform density profile, for which very close to symmetric final states are obtained}. These effects lead to a distinctive correlation of the orientation of the final structure with the distribution of ejected mass, and also with the initial (very small) angular fluctuations.

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.

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.