Generation of angular momentum in cold gravitational collapse
^{1}
Istituto dei Sistemi Complessi Consiglio Nazionale delle Ricerche, via dei
Taurini 19,
00185
Rome,
Italy
email:
Francesco.SylosLabini@roma1.infn.it
^{2}
UPMC Univ. Paris 06, UMR 7585, LPNHE, 75005
Paris,
France
^{3}
CNRS IN2P3, UMR 7585, LPNHE, 75005
Paris,
France
^{4}
Centro Studi e Ricerche Enrico Fermi, Via Panisperna 89 A,
Compendio del Viminale, 00184
Rome,
Italy
^{5}
INFN Unit Rome 1, Dipartimento di Fisica, Universitá di Roma
Sapienza, Piazzale Aldo Moro
2, 00185
Roma,
Italy
^{6}
Faculty of Science and Technology, Rajamangala University of
Technology Suvarnabhumi, Nonthaburi
Campus, 11000
Nonthaburi,
Thailand
Received: 16 June 2015
Accepted: 4 November 2015
During the violent relaxation of a selfgravitating system, a significant fraction of its mass may be ejected. If the timevarying 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 virialised system with nonzero angular momentum. Using numerical simulations we explore here how much angular momentum can be generated in a virialised structure in this way, starting from configurations of cold particles that are very close to spherically symmetric. For the 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 λ ~ 10^{3}, and decreasing slowly with N. For slightly ellipsoidal initial conditions, in which the finiteN fluctuations break the residual reflection symmetries, we observe values λ ~ 10^{2}, i.e. 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 semiaxis of the triaxial relaxed structure and of the ejected mass. This simple mechanism may provide an alternative, or complement, to the socalled tidal torque theory for understanding the origin of angular momentum in astrophysical structures.
Key words: methods: numerical / galaxies: elliptical and lenticular, cD / galaxies: formation
© 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énaultBrunet 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 socalled 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 selfgravitating 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)^{1}. 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 wellapproximated 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_{0}, following the radiusdependent 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 nonspherically symmetric and typically triaxial. We varied the number of particles N from 10^{3} to 10^{5}.

N particles are distributed randomly in a prolate ellipsoidal region. We define a_{3} to be the largest semiprincipal axis and a_{1} = a_{2} the smaller ones. The three eigenvalues of the inertia tensor are simply (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 semiprincipal axes we have Λ_{1} ≥ Λ_{2} ≥ Λ_{3}. The principal axis corresponding to the eigenvalue Λ_{1} is oriented in the direction of the shortest semiprincipal axis a_{1}, while the principal axis associated with Λ_{3} is in the direction of the longest one a_{3}. The particle number N again spans the range from 10^{3} to 10^{5}. The shape at any time may then be characterised by the parameter^{2}(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 Nbody 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 timestep, 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 forcesmoothing 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^{3}. 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 nonexpanding 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 finitesize 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 nonexpanding 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 nonzero 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/, where ρ_{0} is the initial mean density, followed by a reexpansion, 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 nonzero initial velocities in Sylos Labini (2012); and for ellipsoidal initial conditions in Benhaiem & Sylos Labini (2015)^{4}. 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 reexpand, 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 nonspherical 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 LinMestelShu 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 semiprincipal axis, along which particles are, on average, further from the orthogonal plane. In the case of the spherical initial conditions, the nonzero 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 reexpanding 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 nonzero 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 (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 (4)where is the angular momentum of the centre of mass of the bound particles, located at R_{b} and moving with velocity V_{b}, and 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 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 (5)where 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 shellcrossing.
We measure angular momentum in the natural units given by (6)where M is the total (initial) mass of the system, and E_{0} its total (initial) energy (equal to its potential energy). Further, we report results for the angular momentum of the bound mass , given in terms of the socalled spin parameter λ, as defined by Peebles (1969), Knebe & Power (2008): (7)where is (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)^{5}.
3.3. Results for a chosen initial condition
First we present results for the case of spherical initial conditions with α = 1 and N = 10^{5} 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 ; (ii) the flattening ratio ι_{80} 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} , , (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. 

Open with DEXTER 
Fig. 2 Evolution from spherical initial condition with α = 1, and N = 10^{5}: top left: the fraction of particles with positive energy f_{+} (full black line) and the gravitational radius R_{g} (dotted red line) in units of the initial sphere radius R_{0}; top right: the flatness ratio ι_{80} for the 80% most bounded particles; bottom left: the total angular momentum L_{T} (full black line) with respect to the initial centre of mass, the angular momentum of the centre of mass of the bound particles (red dashed line) and the angular momentum of the bound particles (blue dasheddotted line) with respect to their centre of mass (all in units of L_{0} defined in Eq. (6)); bottom right: the spin parameter λ. 

Open with DEXTER 
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 ι_{80}, 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 nonzero 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 offset 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^{6} 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 nonradial velocities, and thus nonzero angular momentum (while the total angular momentum remains zero).
Figure 3 shows, for the case α = 1 and N = 10^{5}, 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 . 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 is unity), as a function of radius, at each of the indicated times, for a simulation with N = 10^{5} 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. 

Open with DEXTER 
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>τ_{c}P(  ℓ ) rapidly approaches the distribution characterising the final virialised state.
Fig. 4 Distribution of the absolute values of the particle angular momentum ℓ (in units of L_{0}) with respect to the centre of mass of bound particles P(  ℓ ), for α = 1 and N = 10^{5} at different times before and after the collapse. 

Open with DEXTER 
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_{0}). 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 (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 , 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 = 10^{5}, 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. 

Open with DEXTER 
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^{5}. 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), ι_{80} (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_{0}) /E_{0} = −E_{p}/E_{0}, 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_{+}(f_{p}), f_{E}, ι_{80} and λ for the spherical initial conditions as a function of α, for N = 10^{4}. The points are the average values for N_{s} = 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 ). The lower panel also shows (red dashed points) the angular momentum of the bound mass normalised to L_{0}. 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. 

Open with DEXTER 
Fig. 7 Final values of f_{+}(f_{p}), f_{E}, ι_{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 L_{0}). 

Open with DEXTER 
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_{0}. 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 (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 = 10^{4}. μ is normalized in units in which GM^{3}/R_{0} = 1. The dashed line is λ ∝ μ. 

Open with DEXTER 
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 ι_{80}(0) < 0.15 that the final ι(t) is linearly proportional to the initial one. Instead, for ι_{80}(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 ι_{80}(0) only for ι_{80}(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 symmetrybreaking. 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^{4} 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 scaledependence of the normalised spin parameter, defined as (Bullock et al. 2001) (11)where 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^{5} 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 = 10^{5}. The distance is normalised to the structure size R_{max}. 

Open with DEXTER 
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, , 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 symmetrybreaking 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 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. 

Open with DEXTER 
Fig. 11 Pdf of the cosine of the angle well after the collapse between 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. 

Open with DEXTER 
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^{4} and N = 10^{5}, the results we reported above are consistent with this finding, while for initial conditions with α> 0, an apparent decrease of the average 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^{3} to 10^{5}. 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 , 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 powerlaw dependence, , with γ ≈ 0.4.
Fig. 12 Logarithm of the final angular momentum of bound virialised mass (in units of L_{0}) as a function of particle number N, for spherical initial conditions with α = 1. The point for the cases N = 10^{5} 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 = 10^{5}. 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}. 

Open with DEXTER 
Figure 13 shows the result for the same data as in the previous figure and, for the different N shown, for the flatness ratio ι_{80}. 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 twobody collisionality. At the very largest N, on the other hand, the result appears to converge to a timeindependent value. In the small range of (larger) N, in which the timedependence that results from twobody collisionality is small, the final value appears to be consistent with an Nindependent 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 since this quantity is essentially frozen, as we have seen, at about t ≈ 2, and thus not modified by the twobody 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^{3} 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 forcesoftening 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 twobody relaxation, which make the system evolve in time towards a more spherical quasiequilibrium. 

Open with DEXTER 
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 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 to decrease as N increases within the significant scatter of different realisations at given N. This is probably because of a weaker Ndependence 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 socalled radial orbit instability. In future work we will study more fully the Ndependence 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 selfgravity 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} = 0.04 ± 0.005. Clearly, while galaxy formation, in particular, involves much more complicated nongravitational 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 Nbody 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 twobody 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 nontrivial correlation, described for example by nontrivial twopoint 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 symmetrybreaking, 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 BSSGC), 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 ANR10EQPX 2901) overseen by the French National Research Agency (ANR), as part of the Investissements d’Avenir programme.
Mergers of virialized structures can also lead to ejection of a significant mass, in particular when the two structures have very different mass (Carucci et al. 2014).
Note that these parameters are generally defined as a function of the semiprincipal axes a_{1},a_{2},a_{3}, rather than as a function of the eigenvalues (see Eq. (1)). However for small deformations of a perfect sphere, which is the case we consider here, these definitions are almost equivalent. The advantage of this definition of the parameters using the eigenvalues is that it can be used to characterise any distribution.
Mass ejection in mergers of two virialised structures is studied in detail in Carucci et al. (2014), Samsing (2015).
Another commonly used normalization for the spin is that introduced in Bullock et al. (2001), defined by with with where W_{b} is the potential energy. For a virialized structure, with E_{b} = W_{b}/ 2, .
References
 Aarseth, S., Lin, D., & Papaloizou, J. 1988, ApJ, 324, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Aguilar, L., & Merritt, D. 1990, ApJ, 354, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Antonov, V. A. 1961, Sov. Astron., 4, 859 [NASA ADS] [Google Scholar]
 Barnes, E. I., Lanzel, P. A., & Williams, L. L. R. 2009, ApJ, 704, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Bellazzini, M., Bragaglia, A., Carretta, E., et al. 2012, A&A, 538, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benhaiem, D., Sylos Labini, F. 2015, MNRAS, 448, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Bett, P., Eke, V., Frenk, C. S., et al. 2007, MNRAS, 376, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Boily, C., Athanassoula, E., & Kroupa, P. 2002, MNRAS, 332, 971 [NASA ADS] [CrossRef] [Google Scholar]
 Boily, C. M., & Athanassoula, E. 2006, MNRAS, 369, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., Dekel, A., & Kolatt, T. S. 2001, ApJ, 555, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Carucci, I. P., Sparre, M., Hansen, S. H., & Joyce, M. 2014, J. Cosmol. Astropart. Phys., 6, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Fridman, A. M., & Polyachenko, V. L. 1984, Physics of gravitating systems. I. Equilibrium and stability (Springer Verlag) [Google Scholar]
 HénaultBrunet, V., Gieles, M., Evans, C. J., et al. 2012, A&A, 545, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henon, M. 1973, Ann. Astrophys., 24, 229 [NASA ADS] [Google Scholar]
 Hernandez, X., Park, C., CervantesSodi, B., & Choi, Y.Y. 2007, MNRAS, 375, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., & Sylos Labini, F. 2013, MNRAS, 429, 1088 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., Marcos, B., & Sylos Labini, F. 2009, MNRAS, 397, 775 [NASA ADS] [CrossRef] [Google Scholar]
 Kacharov, N., Bianchi, P., Koch, A., et al. 2014, A&A, 567, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Knebe, A., & Power, C. 2008, ApJ, 678, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, C. C., Mestel, L., & Shu, F. H. 1965, ApJ, 142, 1431 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D., & Aguilar, L. A. 1985, MNRAS, 217, 787 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P. J. E. 1969, ApJ, 155, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Romanowsky, A. J., & Fall, S. M. 2012, ApJS, 203, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J. 2015, ApJ, 799, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Sylos Labini, F. 2012, MNRAS, 423, 1610 [NASA ADS] [CrossRef] [Google Scholar]
 Sylos Labini, F. 2013a, MNRAS, 429, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Sylos Labini, F. 2013b, A&A, 552, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sylos Labini, F., Benhaiem, D., & Joyce, M. 2015, MNRAS, 449, 4458 [NASA ADS] [CrossRef] [Google Scholar]
 Theis, C., & Spurzem, R. 1999, A&A, 341, 361 [NASA ADS] [Google Scholar]
 van Albada, T. 1982, MNRAS, 201, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Vitvitska, M., Klypin, A. A., Kravtsov, A. V., et al. 2002, ApJ, 581, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Worrakitpoonpon, T. 2015, MNRAS, 466, 1335 [NASA ADS] [CrossRef] [Google Scholar]
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. 

Open with DEXTER  
In the text 
Fig. 2 Evolution from spherical initial condition with α = 1, and N = 10^{5}: top left: the fraction of particles with positive energy f_{+} (full black line) and the gravitational radius R_{g} (dotted red line) in units of the initial sphere radius R_{0}; top right: the flatness ratio ι_{80} for the 80% most bounded particles; bottom left: the total angular momentum L_{T} (full black line) with respect to the initial centre of mass, the angular momentum of the centre of mass of the bound particles (red dashed line) and the angular momentum of the bound particles (blue dasheddotted line) with respect to their centre of mass (all in units of L_{0} defined in Eq. (6)); bottom right: the spin parameter λ. 

Open with DEXTER  
In the text 
Fig. 3 Average in radial shells of the square of the transverse component of particle velocities (in units in which is unity), as a function of radius, at each of the indicated times, for a simulation with N = 10^{5} 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. 

Open with DEXTER  
In the text 
Fig. 4 Distribution of the absolute values of the particle angular momentum ℓ (in units of L_{0}) with respect to the centre of mass of bound particles P(  ℓ ), for α = 1 and N = 10^{5} at different times before and after the collapse. 

Open with DEXTER  
In the text 
Fig. 5 For the case α = 1 and N = 10^{5}, 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. 

Open with DEXTER  
In the text 
Fig. 6 Final values of f_{+}(f_{p}), f_{E}, ι_{80} and λ for the spherical initial conditions as a function of α, for N = 10^{4}. The points are the average values for N_{s} = 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 ). The lower panel also shows (red dashed points) the angular momentum of the bound mass normalised to L_{0}. 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. 

Open with DEXTER  
In the text 
Fig. 7 Final values of f_{+}(f_{p}), f_{E}, ι_{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 L_{0}). 

Open with DEXTER  
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 = 10^{4}. μ is normalized in units in which GM^{3}/R_{0} = 1. The dashed line is λ ∝ μ. 

Open with DEXTER  
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 = 10^{5}. The distance is normalised to the structure size R_{max}. 

Open with DEXTER  
In the text 
Fig. 10 Histograms of the cosine of the angle well after the collapse between 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. 

Open with DEXTER  
In the text 
Fig. 11 Pdf of the cosine of the angle well after the collapse between 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. 

Open with DEXTER  
In the text 
Fig. 12 Logarithm of the final angular momentum of bound virialised mass (in units of L_{0}) as a function of particle number N, for spherical initial conditions with α = 1. The point for the cases N = 10^{5} 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 = 10^{5}. 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}. 

Open with DEXTER  
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 twobody relaxation, which make the system evolve in time towards a more spherical quasiequilibrium. 

Open with DEXTER  
In the text 