© ESO, 2016

1. Introduction

For several decades, observations have shown that galaxies of all types have significant angular momentum whose origin remains a fascinating open theoretical problem (for a review, see, e.g. Romanowsky & Fall 2012). Globular clusters have also been observed to have net rotation more recently (Hénault-Brunet et al. 2012; Bellazzini et al. 2012; Kacharov et al. 2014). The most popular theory to account for angular momentum in virialised structures is the so-called tidal torque theory, in which the virialised structure gains angular momentum by the action of the torque that is caused by the tidal fields generated by surrounding structures (Peebles 1969). Here, we explore a distinct mechanism which, despite its simplicity, appears not to have been considered in the literature, apart from in one recent study of cold spherical collapse (Worrakitpoonpon 2015): generation of angular momentum by ejection of matter in violent relaxation. Indeed violent relaxation of self-gravitating structures is characterised by very large amplitude variations of the mean field potential of a structure, which can give enough energy to particles to escape, even if all the mass is initially bound. The amount of mass ejected depends strongly on the initial conditions, varying from zero to as much as 40% for highly uniform and completely cold initial conditions (Joyce et al. 2009)¹. If, in addition, the mass distribution is not spherically symmetric during the violent phase of the relaxation, this ejected mass could be expected to carry at least some angular momentum – and “generate” an (equal and opposite) amount of angular momentum in the remaining bound mass. In other words the two components of the mass – the part that is ejected and the other part that remains bound – can exert a net torque on one another during the violent relaxation leading to a net angular momentum for both of them. Further, given that violent relaxation starting from a wide range of initially cold conditions is often characterised by a very strong breaking of spherical symmetry, leading to triaxial relaxed structures (e.g. Merritt & Aguilar 1985; Barnes et al. 2009; Benhaiem & Sylos Labini 2015; Sylos Labini et al. 2015) the effect might be far from negligible.

Using numerical simulations, we explore here how much angular momentum can be generated by this mechanism in a virialised structure. More specifically, we consider an initially isolated cold distribution of matter in open boundary conditions and without expansion, and a range of initial spatial distributions, which are very close to spherically symmetric: (i) particles distributed around a centre following a mean density profile, which decays as a power law of the radial distance; or (ii) particles distributed with uniform mean density inside an ellipsoidal region. These are initial conditions that have been extensively studied in the literature (see references below) to investigate the processes that are involved in the formation of galaxies and other astrophysical structures. We note that, in a cosmological context, the simulations can be taken to represent the evolution, in physical coordinates, of a single isolated overdensity with the chosen initial profile. On the other hand, the relation of such simulations to the more general case of the evolution of an overdensity in an expanding universe (which cannot necessarily be well-approximated as isolated) is a more complex issue and will be discussed in detail in a forthcoming publication. Indeed, the collapse in an expanding universe involves several additional parameters that can play a role in the dynamics (e.g. how precisely the background density is modelled) and additional numerical issues (notably concerning the control of numerical accuracy without energy conservation – see Joyce & Sylos Labini 2012; Sylos Labini 2013b).

In the present context the fluctuations that break spherical symmetry in the initial conditions are evidently of central importance for the phenomenon we are studying. For our first class of initial conditions spherical symmetry is broken only by the finite particle number fluctuations, while for the second class these fluctuations break the residual reflection symmetries. We find that almost all these initial conditions – except where there is negligible mass ejection – indeed lead to a measurable net angular momentum in the relaxed virialised structure. The angular momentum is larger by about an order of magnitude in the second class and in many cases is sufficiently large to suggest that the mechanism could potentially account for angular momenta observed in astrophysical structures. Indeed, the typical measured values of the spin parameters of galaxies (see e.g. Hernandez et al. 2007) are of the same order of magnitude as the ones that we find in some of our simulations.

The paper is organised as follows. We first present the details of the numerical simulations and initial conditions in Sect. 2. The in Sect. 3 we then present our results, first describing the relevant features of the evolution qualitatively and then giving the quantitative results. In Sect. 4 we summarize our results and conclude.

2. Numerical simulations

2.1. Initial conditions

The detail of the initial conditions we study are as follows:

• N particles distributed randomly inside a sphere of radius R₀, following the radius-dependent density profile ρ(r) ∝ r^{− α}, where r is the radial distance from the centre and α is a constant. We refer to these as “spherical initial conditions”. Initially, the particles have vanishing velocities. We report here the results for the range of α with 0 ≤ α ≤ 2. We restrict ourselves to this range since, in a previous study (Sylos Labini 2012; Sylos Labini et al. 2015), we observed that for α> 2, there is negligible (≪1%) mass ejection. As shown in the same study, and except for α very close to zero, the evolution leads to virial equilibria that are very non-spherically symmetric and typically triaxial. We varied the number of particles N from 10³ to 10⁵.

• N particles are distributed randomly in a prolate ellipsoidal region. We define a₃ to be the largest semi-principal axis and a₁ = a₂ the smaller ones. The three eigenvalues of the inertia tensor are simply $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{5}}\mathit{M}\mathrm{(}{\mathit{a}}_{\mathit{j}}^{\mathrm{2}}\mathrm{+}{{\mathit{a}}_{\mathit{k}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}& & \end{array}$(1)where M is the total mass of the system and i ≠ j ≠ k and i,j,k = 1,...,3: from the definition of the semi-principal axes we have Λ₁ ≥ Λ₂ ≥ Λ₃. The principal axis corresponding to the eigenvalue Λ₁ is oriented in the direction of the shortest semi-principal axis a₁, while the principal axis associated with Λ₃ is in the direction of the longest one a₃. The particle number N again spans the range from 10³ to 10⁵. The shape at any time may then be characterised by the parameter²$\begin{array}{ccc}\mathit{\iota }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Lambda }}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathrm{\Lambda }}_{\mathrm{3}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{-}\mathrm{1.}& & \end{array}$(2)We refer to these as “ellipsoidal initial conditions”, and report here results for the value of ι in the range from ι(0) = 0.01 to ι(0) = 0.25. The mass ejection and amplification of the spherical symmetry that breaks during the evolution from these initial conditions has been studied in Benhaiem & Sylos Labini (2015).

2.2. Numerical simulations

We have used the N-body code Gadget2 (Springel 2005). All results presented here are for simulations in which energy is conserved to within one tenth of a percent over the timescale evolved, with maximal deviations at any time of less than half a percent (see Benhaiem & Sylos Labini 2015; Sylos Labini et al. 2015, for more details). This accuracy was attained using values of the essential numerical parameters in the Gadget2 code [0.025 for the η parameter, which controls the time-step, and a force accuracy of α_F = 0.001] in the range of typically used values for this code. In the specific cases of spherical initial conditions with α = 0, which is singular in the limit N → ∞, the treatment is as detailed in Joyce et al. (2009; see also Worrakitpoonpon 2015). We also performed extensive tests of the effect of varying the force-smoothing parameter ε, and found that we can obtain very stable results, providing ε is significantly smaller than the minimal size reached by the whole structure during collapse³. We discuss the dependence on particle number N of our results below.

As noted above, our simulations are performed with open boundary conditions and in a non-expanding background, but can be taken to represent the evolution of an isolated overdensity in an expanding universe, i.e. as overdensity which at the initial time has a large density constrast and it is at rest, in physical coordinates. In this respect we could perform the simulations here using the expanding universe version of the code and periodic boundary conditions. The size of the periodic box relative to that of the initial sphere would then fix the relative overdensity represented by the sphere at the initial time (and thus also the time scale for virialisation compared to the Hubble time). Such a simulation gives, as discussed in detail and studied at length in Joyce & Sylos Labini (2012), Sylos Labini (2013b), identical results in physical coordinates, modulo finite-size corrections suppressed by the ratio of the size of the structure to the size of the periodic box, and possible numerical effects (including the effect of smoothing). Indeed, as discussed in Joyce & Sylos Labini (2012), Sylos Labini (2013b), the differences between such simulations and those in open non-expanding space provides information about these kind of numerical effects, and indeed a tool to better control expanding universe simulations. Our “direct” simulations in physical coordinates are preferable because they are simpler and, notably, far easier to control for accuracy (via energy conservation).

Given that we are interested here in the angular momentum, we also test our simulations for their overall conservation and, as detailed below, compare the measured angular momentum of the final bound structure with the numerical error. We will also check systematically for the accuracy of conservation of the total linear momentum.

3. Results

3.1. Qualitative description of mechanism

Numerical studies of evolution from some of these initial conditions, or very similar ones with small but non-zero virial ratios, have been reported extensively in the literature (see, e.g. Henon 1973; van Albada 1982; Aarseth et al. 1988; Theis & Spurzem 1999; Boily et al. 2002; Joyce et al. 2009; Sylos Labini 2012, 2013a; Benhaiem & Sylos Labini 2015; Sylos Labini et al. 2015), with an emphasis, in particular, on the study of the shape and profile of the virialised structure. Because the system is initially cold, in all cases it undergoes a strong collapse in a time of order 1/$\sqrt{\mathit{G}{\mathit{\rho }}_{\mathrm{0}}}$, where ρ₀ is the initial mean density, followed by a re-expansion, which leads rapidly to virialisation of most of the mass. The degree of violence of the collapse, which can be characterised approximately by the maximal contraction the system undergoes, varies with the initial condition. The most extreme case is the spherical case (α = 0), whose collapse is singular in the limit N → ∞.

While both early studies, and most studies since, have focused on the properties of the virialised structure (e.g. its density profile, shape, and velocity distributions) the phenomenon of mass (and energy) ejection has only been considered in detail more recently: for the case of cold uniform spherical conditions (i.e. the case α = 0 here) in Joyce et al. (2009); for the case α = 0 with non-zero initial velocities in Sylos Labini (2012); and for ellipsoidal initial conditions in Benhaiem & Sylos Labini (2015)⁴. Essentially, the amount of mass ejected depends on the degree of violence of the relaxation, which can be quantified roughly by the maximal contraction attained by the system during the collapse phase. The reason for this correlation between the ejection of mass and the violence of collapse is easier to understand when the mechanism for ejection is studied (Joyce et al. 2009; Carucci et al. 2014; Benhaiem & Sylos Labini 2015; Sylos Labini et al. 2015): The ejected particles are essentially those whose fall times to the centre of the structure during collapse are longest. As a consequence, they pass through the centre of the whole structure when it begins to re-expand, which means that they travel into a deeper potential than the one they travel out of. As a result they pick up an energy “kick”. The magnitude of this energy gain is related directly to the strength of the potential at this time, which depends essentially on the extent of the contraction.

For the generation of angular momentum, the second crucial ingredient is the breaking of spherical symmetry, in addition to mass ejection. Indeed mass ejection can occur in a purely radial gravitational field, but will not then generate angular momentum. That cold initial conditions in the class we are considering lead to virialised states, which are very far from spherically symmetric, was observed first by Merritt & Aguilar (1985) and extensively studied in the literature (see, e.g. Aguilar & Merritt 1990; Theis & Spurzem 1999; Boily & Athanassoula 2006; Barnes et al. 2009). This instability of initially spherical systems to relaxation to non-spherical viral equilibria is usually referred to as “radial orbit instability”, because of its apparent link to an instability of equilibrium systems with purely radial orbits, originally proved by Antonov (1961), Fridman et al. (1984). In two recent papers we have studied in detail the development of this asymmetry during the collapse phase for both the classes of initial conditions we study in this paper (spherical initial conditions in Sylos Labini et al. 2015, and ellipsoidal initial conditions in Benhaiem & Sylos Labini 2015). These studies reveal that the same process involved in the energy ejection, which leads to mass ejection, plays a crucial role for many initial conditions in amplifying the symmetry breaking. Indeed, among the late arriving particles, there is also a spread in arrival times as a function of angle, which leads to an energy injection, and thus a spatial distribution of mass, which is also a function of angle. The effect is maximal in enhancing the final asymmetry in the range α ∈ [ 0.5,1.5 ]. It is, on the other hand, relatively suppressed in the limit α = 0 because the particles’ fall times, and consequently their energy changes, are much less strongly correlated with their initial radial position in this case, leading to a “washing out” of the effect. The development of the asymmetry in the mass distribution, and thus in the gravitational potential, is the combined effect of the growth of the initial finite N fluctuations that break spherical symmetry through gravitational instability, in this case known as the Lin-Mestel-Shu instability (Lin et al. 1965), which is then amplified by the energy ejection to particles, and which also leads to mass ejection.

A schematic representation of the mechanism is given in Fig. 1. Any given realisation of our initial conditions has a longest semi-principal axis, along which particles are, on average, further from the orthogonal plane. In the case of the spherical initial conditions, the non-zero effective ellipticity is a finite N effect, while in the ellipsoidal initial conditions, it is dialed by the parameter ι(0). In the first phase of collapse, this initial asymmetry is amplified by the (gravitational) instability of Lin et al. (1965) with the latest collapse occurring along the longest axis. Particles arriving from the corresponding directions pass through the centre as the structure is already re-expanding in the other directions, leading to a greater energy injection to them. As the potential they are traveling is not spherically symmetric, the particles also gain traverse velocities and thus angular momenta with respect to the centre. As there are also fluctuations breaking rotational symmetry, which grow in the course of the collapse, both the radial and transverse velocities will vary as a function of direction and, for example, those ejected in opposite directions will have slightly different angular momentum. The ejected outgoing particles can then carry net non-zero angular momentum, leaving behind the opposite angular momentum in the particles which virialise in the bound structure.

3.2. Measurement of angular momentum

We now focus on the evolution of the angular momentum during the relaxation. More specifically we decompose the total angular momentum L_T at any time into two components $\begin{array}{ccc}{{L}}_{\mathrm{T}}\mathrm{=}{{L}}_{\mathrm{b}}\mathrm{+}{{L}}_{\mathrm{f}}\mathit{,}& & \end{array}$(3)where L_b (L_f) is the total angular momentum of the “bound” (“free”) particles, i.e. which have negative (positive) energy at the given time. In practice, particles’ energies almost never change sign more than once, i.e. the asymptotically ejected particles are those which are tagged as “free” when their energy becomes positive. Further we consider the decomposition of L_b as $\begin{array}{ccc}{{L}}_{\mathrm{b}}\mathrm{=}{{L}}_{\mathrm{b}}^{\mathrm{com}}\mathrm{+}{{L}}_{\mathrm{b}}^{\mathrm{p}}\mathit{,}& & \end{array}$(4)where ${{L}}_{\mathrm{b}}^{\mathrm{com}}\mathrm{=}{\mathit{M}}_{\mathrm{b}}{{R}}_{\mathrm{b}}\mathrm{\times }{{V}}_{\mathrm{b}}$ is the angular momentum of the centre of mass of the bound particles, located at R_b and moving with velocity V_b, and ${{L}}_{\mathrm{b}}^{\mathrm{p}}$ is the angular momentum of the bound particles with respect to their centre of mass. It is the latter which interests us most but we also monitor ${{L}}_{\mathrm{b}}^{\mathrm{com}}$ and L_T, measured with respect to the fixed origin at the centre of mass of the initial configuration, to assess the accuracy of the simulation for which the total angular momentum should be conserved (and equal to zero).

We measure time in units of $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\sqrt{\frac{\mathrm{3}\mathit{\pi }}{\mathrm{32}\mathit{G}{\overline{\mathit{\rho }}}_{\mathrm{0}}}}\mathit{,}& & \end{array}$(5)where ${\overline{\mathit{\rho }}}_{\mathrm{0}}$ is the total initial mass density, i.e. the total mass divided by the volume of the initial sphere or ellipsoid. This corresponds to the time for a particle that is initially at the outer periphery of a cold sphere with this mass density to fall to its centre in the continuum approximation (i.e. taking N → ∞ and keeping the initial mass density profile fixed) and without shell-crossing.

We measure angular momentum in the natural units given by $\begin{array}{ccc}{\mathit{L}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\sqrt{\mathrm{\left|}{\mathit{E}}_{\mathrm{0}}\mathrm{\right|}}}\mathit{,}& & \end{array}$(6)where M is the total (initial) mass of the system, and E₀ its total (initial) energy (equal to its potential energy). Further, we report results for the angular momentum of the bound mass ${{L}}_{\mathrm{b}}^{\mathrm{p}}$, given in terms of the so-called spin parameter λ, as defined by Peebles (1969), Knebe & Power (2008): $\begin{array}{ccc}\mathit{\lambda }\mathrm{=}\frac{\mathrm{\left|}{{L}}_{\mathrm{b}}^{\mathrm{p}}\mathrm{\right|}}{{\mathit{L}}_{\mathrm{0}}^{\mathrm{b}}}\mathit{,}& & \end{array}$(7)where ${\mathit{L}}_{\mathrm{0}}^{\mathrm{b}}$ is $\begin{array}{ccc}{\mathit{L}}_{\mathrm{0}}^{\mathit{b}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{b}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\sqrt{\mathrm{\left|}{\mathit{E}}_{\mathrm{b}}\mathrm{\right|}}}\mathit{,}& & \end{array}$(8)where M_b is the bound mass and E_b (W_b) the total (gravitational potential) energy of this mass (with respect to its center of mass)⁵.

3.3. Results for a chosen initial condition

First we present results for the case of spherical initial conditions with α = 1 and N = 10⁵ particles. We consider this case since its evolution is typical of those we observe for all our different initial conditions. For this case, the panels of Fig. 2 show as a function of time (i) the fraction f₊ of the initial mass in free particles and the gravitational radius defined as ${\mathit{R}}_{\mathrm{g}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{b}}^{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{\left|}{\mathit{W}}_{\mathrm{b}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right|}}$; (ii) the flattening ratio ι₈₀ as defined by (2) with the subscript indicating that the inertia tensor was calculated using the particles with energy less than 20% of the energy of the most bound particle; (iii) the magnitude of the angular momenta | L_T |, $\mathrm{\left|}{{L}}_{\mathrm{b}}^{\mathrm{com}}\mathrm{\right|}$, $\mathrm{\left|}{{L}}_{\mathrm{b}}^{\mathrm{p}}\mathrm{\right|}$ (on a logarithmic scale); and (iv) the spin parameter λ.

Fig. 1 Schematic representation of the generation of angular momentum in the ejected particles. The upper panel shows schematically symmetry breaking in the initial cloud. When particles are ejected, both their radial and traverse components generally have some dependence on angle, leading to a net angular momentum being carried away by the ejected particles.

Fig. 2 Evolution from spherical initial condition with α = 1, and N = 105: top left: the fraction of particles with positive energy f+ (full black line) and the gravitational radius Rg (dotted red line) in units of the initial sphere radius R0; top right: the flatness ratio ι80 for the 80% most bounded particles; bottom left: the total angular momentum LT (full black line) with respect to the initial centre of mass, the angular momentum of the centre of mass of the bound particles ${\mathit{L}}_{\mathrm{b}}^{\mathrm{com}}$ (red dashed line) and the angular momentum of the bound particles ${\mathit{L}}_{\mathrm{b}}^{\mathit{p}}$ (blue dashed-dotted line) with respect to their centre of mass (all in units of L0 defined in Eq. (6)); bottom right: the spin parameter λ.

As described above, the first plot illustrates clearly that at τ_c ≈ 1, the system reaches its maximal contraction, and the strong energy injection, which leads to particle ejection, occurs in a very short interval around this time, which corresponds approximately to the dispersion in fall times of the mass. Likewise, the second panel, that shows the behaviour of ι₈₀, we see that the breaking of spherical symmetry grows monotonically from the beginning and is then strongly amplified around the maximal collapse (as described in detail in Sylos Labini et al. 2015).

As anticipated, and as can be seen from the lower two plots, the bound mass indeed moves towards a stationary state with non-zero angular momentum (relative to its centre of mass). This final angular momentum of the bound mass, albeit small in the natural units, is clearly not numerical in origin; it is about three orders larger than the final total angular momentum, which gives a measure of the violation of angular momentum conservation that is due to numerical error. Further we have performed several tests, varying the essential numerical parameters, and found our results to be stable. On the other hand, the total angular momentum L_b of the bound mass (with respect to the initial centre of mass of the whole system) is dominated by its translational motion: the bound mass “recoils” with a net linear momentum, compensating for that of the ejected mass along an axis which is off-set from the initial centre of mass. We give some further details below on this ejected linear momentum, which is measured numerically to a precision of order 10^-6 (i.e. the measured momentum change of these particles is 10⁶ times larger than the change in the total momentum due to numerical error).

3.4. Generation of angular momentum during collapse

It is interesting to follow the dynamics leading to the generation of net angular momentum in the collapse phase. During this phase the mass which is finally ejected has not yet undergone the energy boost which leads to this ejection, but the gravitational potential has already developed very significant asymmetry, which leads the particles to have non-radial velocities, and thus non-zero angular momentum (while the total angular momentum remains zero).

Figure 3 shows, for the case α = 1 and N = 10⁵, the measured dispersion of the transverse velocities in radial shells, as a function of radius, and at different times during the evolution. The velocity is normalised in units of ${\mathit{v}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{0}}}$. We observe that the transverse velocities are already at t = 0.5τ_c comparable with, and at t = 0.75τ_c even larger than, their magnitude in the final virialised structure. We note that these are the transverse velocities produced by the growth of the initial density fluctuations that are breaking spherically symmetry through gravitational instability, well before the energy injection to particles that fall around t ≈ τ_c leading to the mass (and angular momentum) ejection. For longer time scales, i.e. t>τ_c, the outer shell, corresponding to very weakly bound particles with energy very close to zero, is expanding since these particles are still travelling outwards on their large orbits. On the other hand, the inner shell has reached a stationary state, as can be seen by comparing the velocity profiles at t = 2.5τ_c and at t = 5τ_c.

Fig. 3 Average in radial shells of the square of the transverse component of particle velocities (in units in which ${\mathit{v}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{0}}}$ is unity), as a function of radius, at each of the indicated times, for a simulation with N = 105 particles of spherical initial conditions with α = 1. Note that already well before the system reaches it maximal contraction, at t ≈ 1, the Poissonian fluctuations breaking spherical symmetry have been amplified to produce tangential velocities larger than in the final virialized state.

Figure 4 shows the distribution P( | ℓ |) of the absolute values of the particle angular momentum ℓ with respect to the centre of mass of bound particles, at the indicated times (before and after the collapse). In line with the results on the transverse velocities, we observe that significant angular momentum is already generated before the completion of the collapse phase. Further, for t>τ_cP( | ℓ |) rapidly approaches the distribution characterising the final virialised state.

Fig. 4 Distribution of the absolute values of the particle angular momentum ℓ (in units of L0) with respect to the centre of mass of bound particles P( | ℓ |), for α = 1 and N = 105 at different times before and after the collapse.

The upper panel of Fig. 5 shows the temporal evolution of the average value of the modulus of the particle angular momentum for the bound particles (with respect to their centre of mass, in units of L₀). In line with what would be expected from the plots of the transverse velocities, we observe a constant growth during the collapse phase, followed by a very significant boost around the time of maximal contraction, when the energy ejection occurs. Thus we see the angular momentum that is generated by the collapse receiving a very significant boost in the final phase of the collapse. This is further quantified by the lower panel of Fig. 5, which shows the evolution, as a function of time, of $\begin{array}{ccc}{\mathit{l}}_{\mathrm{b}}^{\mathrm{p}}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}}{{\sum }_{\mathit{i}}\mathrm{|}\mathit{\ell }{\mathrm{|}}_{\mathit{i}}}\mathit{,}& & \end{array}$(9)i.e. the modulus of the angular momentum of bound particles, normalised to the sum of the moduli of the angular momenta of the same particles. The dashed horizontal line corresponds to the value $\sqrt{{\mathit{N}}_{\mathrm{+}}}\mathit{/}{\mathit{N}}_{\mathrm{-}}\mathrm{\approx }{\mathit{f}}_{\mathrm{+}}\mathit{/}\sqrt{\mathit{N}}$, which is the amplitude of the bound angular momentum which would be expected if the ejection operated as a simple random removal of particles without any modification of their angular momentum. We observe that the final value is, in fact, at least several times larger. Thus, during the final phase of collapse, when the large energy changes that gives rise to ejection are occurring, the ejected and bound particles exert a torque on one another that very significantly amplifies their final (equal and opposite) angular momenta.

Fig. 5 For the case α = 1 and N = 105, upper panel: evolution of the average of the modulus of the angular momenta of bound particles (with respect to their centre of mass) as a function of time. The lower panel: quantity in Eq. (9) and the dashed line the estimated value expected if the ejection were a random sampling of the total mass.

3.5. Results for different initial conditions (fixed N)

In the above, for simplicity we have focused on the single case of spherical initial conditions with α = 1, and we have shown results for N = 10⁵. The qualitative behaviours we have discussed in this case, for what concerns the generation of angular momentum, are shared by all our other initial conditions. Quantitatively, there is a variation of the final angular momentum, which depends notably on the details of the mass ejection and on the degree of symmetry breaking during the collapse. The precision of the simulations also varies from case to case; simulations for the spherical initial conditions with α = 0 are particularly delicate since the collapse is the most extreme in this case, with a singular behaviour in the limit N → ∞. Correspondingly, in this case the measured angular momentum is only a few times larger that the numerical error in the conservation of the total angular momentum, while, as shown above (see Fig. 2), it is almost three orders of magnitude larger for the case α = 1.

Figure 6 shows the “final” values of f₊ (top panel), ι₈₀ (third panel) and λ (bottom panel), i.e. at t = 5τ_c, for the spherical initial conditions as a function of α. In addition, in the second panel, we also show f_E, the “fraction of ejected energy”, defined by f_E = (E_b − E₀) /E₀ = −E_p/E₀, where E_p is (by energy conservation) the total ejected energy in the frame in which the particles are initially at rest, i.e. the sum of the kinetic energy of the centre of mass of the bound mass and the kinetic energy of the ejected particles. Figure 7 shows the same four quantities for the ellipsoidal initial conditions as a function of ι(0).

Fig. 6 Final values of f+(fp), fE, ι80 and λ for the spherical initial conditions as a function of α, for N = 104. The points are the average values for Ns = 20 realisations for each α ≠ 0 and for 47 realisations in the case α = 0; the indicated error bar is the corresponding standard deviation σ (for λ we show the error on the mean $\mathit{\sigma }\mathit{/}\sqrt{{\mathit{N}}_{\mathrm{s}}}$). The lower panel also shows (red dashed points) the angular momentum of the bound mass normalised to L0. For the smallest values of α, the angular momentum generated actually decreases but, because of the increasing mass (and energy) ejection, the characteristic angular momentum of the final bound mass decreases even faster.

Fig. 7 Final values of f+(fp), fE, ι80 and λ for ellipsoidal initial conditions, as a function of ι(0). The points are averages for 20 realisations and the error bar, the corresponding standard deviation (for λ we show the error on the mean). As in the previous figure the lower panel also displays (red points, dashed line) the final angular momentum normalised to L0).

Comparing these two figures, the most striking result is that there is an order of magnitude difference in the final angular momentum (as characterised by the spin parameter λ) for most of the spherical initial conditions with α ≠ 0, and the ellipsoidal initial conditions with ι(0) ≠ 0, with a sharp interpolation between the two classes around the case α = 0 (which indeed is also the limit ι(0) → 0 of the ellipsoidal initial conditions). In both figures, we also show in the lower panel the final angular momentum normalised to L₀. Compared with λ, we see that around α = 0, where the relaxation is most violent, there is a significant difference. This arises simply from the fact that there is a considerable ejected mass and energy in this case. Indeed the relative amplitude of the two quantities can be expressed as $\begin{array}{ccc}\frac{\mathit{\lambda }}{{\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}}\mathrm{=}{\left(\frac{\mathit{M}}{{\mathit{M}}_{\mathrm{b}}}\right)}^{\mathrm{5}\mathit{/}\mathrm{2}}{\left(\frac{\mathrm{\left|}{\mathit{E}}_{\mathrm{b}}\mathrm{\right|}}{\mathrm{\left|}{\mathit{E}}_{\mathrm{0}}\mathrm{\right|}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{f}}_{\mathrm{E}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{+}}{\mathrm{)}}^{\mathrm{5}\mathit{/}\mathrm{2}}}\mathrm{·}& & \end{array}$(10)Thus, the final spin includes an amplification which comes from the change in the characteristic scale of angular momentum that results from the mass and energy ejection (which leads to a more bound, but less massive, structure than the initial condition).

Fig. 8 Final spin parameter λ, and the modulus of final total linear momentum of the bound particles μ, in all our (322) realizations of initial conditions (both spherical and ellipsoidal) with N = 104. μ is normalized in units in which GM3/R0 = 1. The dashed line is λ ∝ μ.

It is clear, however, from the two plots that this amplification from the normalisation is not the main factor explaining the much larger final angular momenta from the ellipsoidal initial conditions. In fact, the main effect, which amplifies the spin in the ellipsoidal initial conditions compared to the spherical ones, comes clearly from the essential difference in the initial conditions: the stronger initial breaking of the spherical symmetry. This initial asymmetry leads both to a more anisotropic collapse, compared to that in the spherical case and, at the same time, the ejection of a comparable amount of mass to that in the small α spherical models. As described in detail in Benhaiem & Sylos Labini (2015), the particles on the stretched axis of the initial condition fall slightly later and receive a large energy boost. We note that it is only for ι₈₀(0) < 0.15 that the final ι(t) is linearly proportional to the initial one. Instead, for ι₈₀(0) > 0.15, substructures form during the collapse, leading to a substantial difference in the shape of the virialised structure, see Benhaiem & Sylos Labini (2015). Similarly the percentage of ejected particles decreases linearly with ι₈₀(0) only for ι₈₀(0) < 0.15.

Besides the angular momentum itself, the generation of a net relative motion of the centre of masses of the final bound and ejected mass is also a direct measure of symmetry-breaking. We would thus expect to possibly find a similar trend in the final linear momentum of the ejected (or bound) particles. Figure 8 shows a plot, for all our simulations with N = 10⁴ particles, of both λ and μ, the modulus of the “final” linear momentum of the ejected particles (equal and opposite to that of the bound mass up to a numerical precision that is smaller by several orders of magnitude in all cases.) We observe that there is indeed a clear correlation between the two quantities. Furthermore, it appears highly consistent with a roughly linear relation, as one might expect, given that they are both indirect measures of an initially small break in spherical symmetry.

Finally, we have also studied the scale-dependence of the normalised spin parameter, defined as (Bullock et al. 2001) $\begin{array}{ccc}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\sqrt{\mathrm{2}\mathit{GM}\mathrm{(}\mathit{r}{\mathrm{)}}^{\mathrm{3}}\mathit{r}}}\mathit{,}& & \end{array}$(11)where ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ is the angular momentum of bound particles with distance <r from their centre of mass, and M(r) is the mass enclosed in a sphere of radius r. The results for five different sets of simulations with α = 0,0.5,1,1.5,2, and N = 10⁵ are shown in Fig. 9. The profiles are similar in all cases, but their amplitude depends on α as highlighted in Fig. 6. The fact that λ′(r) is larger for α = 0 than for α> 0 is related to the fact that orbits are slightly more radial in the former case.

Fig. 9 Spin parameter λ′(r) (see Eq. (11)), averaged over 20 realisations, for α = 0,0.5,1,1.5,2, and N = 105. The distance is normalised to the structure size Rmax.

3.6. Correlation of angular momentum with final spatial distribution

As analysed in detail in Benhaiem & Sylos Labini (2015), Sylos Labini et al. (2015), for both classes of initial conditions, there is a correlation between the direction in which mass is ejected and the longest axis of both the initial condition and the final bound mass. The latter of these two are strongly correlated because it is the growth of the initial fluctuations breaking spherical symmetry, amplified by the mass ejection, which leads to the final spatial asymmetry. The ejection of mass is amplified along these same directions because they are the particles which fall latest and, as a result, pick up a large energy kick. Given that we have seen that most of the final angular momentum is generated at the time of ejection, we would expect that it to preferably be aligned orthogonal to the preferential axis for ejection.

Figures 10 and 11 show histograms for 20 realisations of each indicated initial condition, of the modulus of the cosine of the angle between the final angular momentum of the bound mass, ${{L}}_{\mathrm{b}}^{\mathrm{p}}$, and the eigenvector corresponding to the longest axis of the final bound mass (blank histograms) and the final ejected mass (filled histograms). For the ellipsoidal initial conditions – which very strongly single out an initial preferred axis and which remains strongly correlated with the final longest axis – the anticipated correlation is very clearly present. In the spherical initial conditions, the correlation is weaker but nevertheless visible, except perhaps in the case α = 0. Indeed in this latter case, as discussed in Sylos Labini et al. (2015), the symmetry-breaking in the final state is, in fact, much weaker, as is the correlation between the final and initial asymmetry. In this case, with the fall time of all particles being the same in the limit N → ∞, there is a much weaker correlation between the initial radial position of a particle and the energy change it undergoes in the violent collapse.

Fig. 10 Histograms of the cosine of the angle well after the collapse between ${{L}}_{\mathrm{b}}^{\mathrm{p}}$ and the eigenvector that corresponds to the longest principal axis of the final bound mass (blank histograms), and of the ejected mass (filled histograms) for 20 realisations of spherical initial conditions for different α indicated in the caption.

Fig. 11 Pdf of the cosine of the angle well after the collapse between ${{L}}_{\mathrm{b}}^{\mathrm{p}}$ and the eigenvector corresponding to the longest principal axis of the final bound (blank histograms), and of the ejected mass (filled histograms) for 20 realisations of ellipsoidal initial conditions for different ι(0), indicated in the caption.

3.7. Dependence on N

The N dependence of the angular momentum that was generated by mass ejection has been explored in Worrakitpoonpon (2015) for the case α = 0, and evidence for a monotonic decrease ~N^{− β}, with a best fit around β = 1/3, was found. In the case of simulations with N = 10⁴ and N = 10⁵, the results we reported above are consistent with this finding, while for initial conditions with α> 0, an apparent decrease of the average ${\mathit{L}}_{\mathrm{p}}^{\mathrm{b}}$ as N grows is also observed. To explore this further, we performed a more detailed study of the case α = 1, performing a larger number of simulations for a greater number of values of N in the range 10³ to 10⁵. In this case, we also evolved each simulation for a longer time to study the stability of our “final” measures of the different quantities and to check for the possible role of finite N effects.

Figure 12 shows the result for the measured angular momentum ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}$, at the indicated times, in a series of simulations for the different indicated values of N. As anticipated, in the cases where different (and longer) times have been analysed, the measured angular momentum shows no significant dependence on the time: after a few dynamical times there is negligible interaction between the ejected and bound mass and the angular momentum of each is constant. On the other hand, as a function of N, there is a clear monotonic decrease, which is well fitted by a power-law dependence, ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}\mathrm{\propto }{\mathit{N}}^{\mathrm{-}\mathit{\gamma }}$, with γ ≈ 0.4.

Fig. 12 Logarithm of the final angular momentum of bound virialised mass ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}$ (in units of L0) as a function of particle number N, for spherical initial conditions with α = 1. The point for the cases N = 105 corresponds to the results reported above, at t = 5, while for the other N there is data at three times, extending to t ≈ 14. Each point is an average over M realisations, with M = 90 for N = 1000, M = 50 for N = 2000, M = 40 for N = 4000, M = 30 for N = 8000, M = 6 for N = 16 000, M = 4 for N = 32 000, and M = 4 for N = 105. The error bars indicate the corresponding estimated error on the mean. The continuous line (∝ N-0.4) is a best fit to the data, while the dashed lines correspond to N− 1/2 and N− 1/3.

Figure 13 shows the result for the same data as in the previous figure and, for the different N shown, for the flatness ratio ι₈₀. In contrast, in this case we see that the value measured, for all but the largest value of N, shows visible dependence on time. More specifically the bound mass evolves towards a more spherical distribution in time, at a rate which is clearly faster for the smaller N. This kind of N dependence of the timescale of the evolution of the virialised structure is clearly qualitatively coherent with collisional relaxation and, indeed, a quantitative analysis of this type of behaviour from similar initial conditions in Theis & Spurzem (1999) shows that the timescale of this evolution is, indeed, in agreement with that predicted by two-body collisionality. At the very largest N, on the other hand, the result appears to converge to a time-independent value. In the small range of (larger) N, in which the time-dependence that results from two-body collisionality is small, the final value appears to be consistent with an N-independent behaviour, but clearly extrapolation to larger N would be needed to reach a firmer conclusion. However, we underline that we do not have this problem for the above determination of the scaling with N of ${\mathit{L}}_{\mathrm{p}}^{\mathrm{b}}$ since this quantity is essentially frozen, as we have seen, at about t ≈ 2, and thus not modified by the two-body collisionality at longer times. We note that Joyce et al. (2009) report tests for collisionality during the collapse phase, for the very extreme case of a flat (α = 0) profile, and conclude, in line with the study of Aarseth et al. (1988), for N> 10³ that collisional effects during the collapse phase are indeed negligible. This result, consistent with theoretical prediction for the collisional timescale, is borne out by the observation that the evolution of the macroscopic features of the collapse, i.e. potential energy, viral ratio, size of structure and, in particular, the fraction of the mass ejected, are stable when N increases and, likewise, when the gravitational force-softening is varied, provided it is sufficiently small to resolve the mean gravitational field.

Fig. 13 Flatness ratio ι80 as a function of particle number N, for the same simulations in the previous figure (α = 1) for the indicated values of N. For this quantity the observed N dependence is due to two-body relaxation, which make the system evolve in time towards a more spherical quasi-equilibrium.

Given that, as we have emphasized, the breaking of spherical symmetry in the initial conditions is due to finite N fluctuations, the fact we find that the final angular momentum ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}$ decreases as N increases, and in principle goes to zero as N → ∞, is what one would naively expect. A similar decrease in amplitude, with an exponent of about the same value, has been observed for the case α = 0 in Worrakitpoonpon (2015). In our simulations of initial conditions from elliptical conditions, we have not been able to detect a clear trend for ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}$ to decrease as N increases within the significant scatter of different realisations at given N. This is probably because of a weaker N-dependence in this case, as the spherical symmetry is “mostly” broken by ι(0) ≠ 0, which does not depend on N. In contrast, the lack of evidence for any N dependence in the final ι – once N is sufficiently large so that collisional relaxation plays no role on the time scales probed – may be indicative of the true physical behaviour: in the limit of exact spherical symmetry of a perfectly cold system, it is expected that symmetry will be broken because of a so-called radial orbit instability. In future work we will study more fully the N-dependence of both the breaking of symmetry and that of the angular momentum of the ejected particles, and their connection to one another.

4. Discussion and conclusions

The origin of the angular momentum in astrophysical structures – notably that measured observationally in galaxies – is a fascinating open problem in astrophysics and cosmology. In this article, we have discussed and studied a way, which is different to commonly considered ones, such as tidal torques, in which angular momentum can be generated by self-gravity only. This can occur, in principle, for initial conditions which evolve sufficiently violently, with large variations of the gravitational mean field, so that mass is ejected during relaxation towards a virialised state. Using numerical simulations, we have shown for two simple broad classes of initial conditions of this kind that significant angular momentum can, indeed, be generated in this way. In most cases this angular momentum is orders of magnitude larger than the angular momentum induced by numerical errors.

For illustrative purposes and for simplicity, we have considered initial conditions very close to spherical symmetry. Indeed, in all our initial conditions rotational symmetry is only fully broken by the Poissonian density fluctuations that are associated with the finite number of particles. Interestingly we find that, despite the “artificiality” of our initial conditions, which are not motivated by any detailed specific physical model, we obtain angular momenta that are of comparable size, i.e. only smaller by a factor of two in some cases, than those typically estimated observationally for galaxies. Indeed we have obtained values of the normalised “spin parameter” up to λ ~ 0.02, while, for example, Hernandez et al. (2007) estimate its average value in the 11 597 spirals and elliptical galaxies observed by the Sloan Digital Sky Survey to be λ₀ = 0.04 ± 0.005. Clearly, while galaxy formation, in particular, involves much more complicated non-gravitational processes than is described by the simple models we have considered here and, moreover, our specific initial conditions are quite ad hoc (and not directly motivated by a cosmological model) it is intriguing that the order of magnitude of the values we obtain are in line with those observed.

We note that studies of dark matter halos in cosmological N-body simulations (see e.g. Vitvitska et al. 2002; Bullock et al. 2001; Bett et al. 2007) give rise to similar values of spin parameter (λ ≈ 0.04), i.e. of the same order as those observed in elliptical galaxies, and as those we find in our simulations in certain cases. As dark matter halos of standard cold dark matter cosmologies form through a hierarchical process rather than in the kind of monolithic collapse we have studied, the mechanism we have studied has likely no relevance in this particular context. Indeed it is the strong violence of the collapse leading to mass ejection which is crucial, and the formation of halos in a cold dark matter cosmology setting is not of this kind. Nevertheless, as discussed in Carucci et al. (2014), Samsing (2015), significant mass ejection can occur in a cosmological setting when halos merge, and it would be interesting to investigate whether this could also lead to the generation of angular momentum at a significant level, compared to the processes of mass accretion, for example, which has been argued in Vitvitska et al. (2002), to account for the spin of dark matter halos.

Just as in cosmological simulations we underline that, essentially, we are simulating the collisionless regime of the gravitational dynamics (apart from the two-body collisional effects which we observed to be present at longer times in smaller N simulations). In so far as this is true, the particle number N is relevant, in principle, in the properties of the final state because it fixes the amplitude of the initial fluctuations. Indeed these fluctuations are, for our initial conditions in this paper, simply Poissonian, and thus their amplitude at any scale (and all their statistical properties) are regulated by this single parameter. Alternatively one could consider setting up an initial particle configuration with the same average density profile but with statistical fluctuations with non-trivial correlation, described for example by non-trivial two-point correlation properties. In this case one could then, in principle, vary N while keeping the relevant fluctuations fixed and obtain results independent of N. As the fluctuations play an essential role in the symmetry-breaking, we expect that the detailed nature of these initial fluctuations may have an impact on the properties of the final structure. We postpone a detailed study of this interesting issue for a future work. We will also explore whether cold but more irregular/clumpy initial conditions can produce values of the spin of the order of magnitude to those observed for galaxies. Furthermore, we will explore whether the kind of correlation we discussed briefly in the last part of the paper – between the shape of the asymmetric virialized structure and its angular momentum – could be used to find observational evidence for or against the role of very violent relaxation in generating the structure of galaxies.

Numerical simulations have been run on the Cineca PLX cluster (project ISCRA BSS-GC), and on the HPC resources of The Institute for Scientific Computing and Simulation financed by Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX- 29-01) overseen by the French National Research Agency (ANR), as part of the Investissements d’Avenir programme.

References

All Figures

	Fig. 1 Schematic representation of the generation of angular momentum in the ejected particles. The upper panel shows schematically symmetry breaking in the initial cloud. When particles are ejected, both their radial and traverse components generally have some dependence on angle, leading to a net angular momentum being carried away by the ejected particles.
In the text

Fig. 2 Evolution from spherical initial condition with α = 1, and N = 105: top left: the fraction of particles with positive energy f+ (full black line) and the gravitational radius Rg (dotted red line) in units of the initial sphere radius R0; top right: the flatness ratio ι80 for the 80% most bounded particles; bottom left: the total angular momentum LT (full black line) with respect to the initial centre of mass, the angular momentum of the centre of mass of the bound particles ${\mathit{L}}_{\mathrm{b}}^{\mathrm{com}}$ (red dashed line) and the angular momentum of the bound particles ${\mathit{L}}_{\mathrm{b}}^{\mathit{p}}$ (blue dashed-dotted line) with respect to their centre of mass (all in units of L0 defined in Eq. (6)); bottom right: the spin parameter λ.

In the text

Fig. 3 Average in radial shells of the square of the transverse component of particle velocities (in units in which ${\mathit{v}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{0}}}$ is unity), as a function of radius, at each of the indicated times, for a simulation with N = 105 particles of spherical initial conditions with α = 1. Note that already well before the system reaches it maximal contraction, at t ≈ 1, the Poissonian fluctuations breaking spherical symmetry have been amplified to produce tangential velocities larger than in the final virialized state.

In the text

	Fig. 4 Distribution of the absolute values of the particle angular momentum ℓ (in units of L0) with respect to the centre of mass of bound particles P( | ℓ |), for α = 1 and N = 105 at different times before and after the collapse.
In the text

	Fig. 5 For the case α = 1 and N = 105, upper panel: evolution of the average of the modulus of the angular momenta of bound particles (with respect to their centre of mass) as a function of time. The lower panel: quantity in Eq. (9) and the dashed line the estimated value expected if the ejection were a random sampling of the total mass.
In the text

Fig. 6 Final values of f+(fp), fE, ι80 and λ for the spherical initial conditions as a function of α, for N = 104. The points are the average values for Ns = 20 realisations for each α ≠ 0 and for 47 realisations in the case α = 0; the indicated error bar is the corresponding standard deviation σ (for λ we show the error on the mean $\mathit{\sigma }\mathit{/}\sqrt{{\mathit{N}}_{\mathrm{s}}}$). The lower panel also shows (red dashed points) the angular momentum of the bound mass normalised to L0. For the smallest values of α, the angular momentum generated actually decreases but, because of the increasing mass (and energy) ejection, the characteristic angular momentum of the final bound mass decreases even faster.

In the text

	Fig. 7 Final values of f+(fp), fE, ι80 and λ for ellipsoidal initial conditions, as a function of ι(0). The points are averages for 20 realisations and the error bar, the corresponding standard deviation (for λ we show the error on the mean). As in the previous figure the lower panel also displays (red points, dashed line) the final angular momentum normalised to L0).
In the text

	Fig. 8 Final spin parameter λ, and the modulus of final total linear momentum of the bound particles μ, in all our (322) realizations of initial conditions (both spherical and ellipsoidal) with N = 104. μ is normalized in units in which GM3/R0 = 1. The dashed line is λ ∝ μ.
In the text

	Fig. 9 Spin parameter λ′(r) (see Eq. (11)), averaged over 20 realisations, for α = 0,0.5,1,1.5,2, and N = 105. The distance is normalised to the structure size Rmax.
In the text

	Fig. 10 Histograms of the cosine of the angle well after the collapse between ${{L}}_{\mathrm{b}}^{\mathrm{p}}$ and the eigenvector that corresponds to the longest principal axis of the final bound mass (blank histograms), and of the ejected mass (filled histograms) for 20 realisations of spherical initial conditions for different α indicated in the caption.
In the text

	Fig. 11 Pdf of the cosine of the angle well after the collapse between ${{L}}_{\mathrm{b}}^{\mathrm{p}}$ and the eigenvector corresponding to the longest principal axis of the final bound (blank histograms), and of the ejected mass (filled histograms) for 20 realisations of ellipsoidal initial conditions for different ι(0), indicated in the caption.
In the text

Fig. 12 Logarithm of the final angular momentum of bound virialised mass ${\mathit{L}}_{\mathrm{b}}^{\mathrm{p}}$ (in units of L0) as a function of particle number N, for spherical initial conditions with α = 1. The point for the cases N = 105 corresponds to the results reported above, at t = 5, while for the other N there is data at three times, extending to t ≈ 14. Each point is an average over M realisations, with M = 90 for N = 1000, M = 50 for N = 2000, M = 40 for N = 4000, M = 30 for N = 8000, M = 6 for N = 16 000, M = 4 for N = 32 000, and M = 4 for N = 105. The error bars indicate the corresponding estimated error on the mean. The continuous line (∝ N-0.4) is a best fit to the data, while the dashed lines correspond to N− 1/2 and N− 1/3.

In the text

	Fig. 13 Flatness ratio ι80 as a function of particle number N, for the same simulations in the previous figure (α = 1) for the indicated values of N. For this quantity the observed N dependence is due to two-body relaxation, which make the system evolve in time towards a more spherical quasi-equilibrium.
In the text