# Dissipationless collapse of a set of N massive particles

Jérôme Perez, Fabrice Roy

► To cite this version:

Jérôme Perez, Fabrice Roy. Dissipationless collapse of a set of N massive particles. Monthly Notices of the Royal Astronomical Society, 2004, 348 (1), pp.62-72. 10.1111/j.1365-2966.2004.07294.x. hal-00988993

HAL Id: hal-00988993

https://ensta-paris.hal.science/hal-00988993

Submitted on 20 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dissipationless collapse of a set of N massive particles

Fabrice Roy<fnref n="1" />* and Jérôme Perez<fnref n="1" />,²*

¹École Nationale Supérieure de Techniques Avancées, Unité de Mathématiques Appliquées, 32 Bd Victor, 75015 Paris, France

²Laboratoire de l'Univers et de ses THéories, Observatoire de Paris-Meudon, 5 place Jules Janssen, 92350 Meudon, France

Accepted 2003 October 13. Received 2003 October 10; in original form 2003 January 13

## ABSTRACT

The formation of self-gravitating systems is studied by simulating the collapse of a set of N particles which are generated from several distribution functions. We first establish that the results of such simulations depend on N for small values of N. We complete a previous work by Aguilar & Merritt concerning the morphological segregation between spherical and elliptical equilibria. We find and interpret two new segregations: one concerns the equilibrium core size and the other the equilibrium temperature. All these features are used to explain some of the global properties of self-gravitating objects: origin of globular clusters and central black hole or shape of elliptical galaxies.

**Key words:** methods: numerical – methods: N-body simulations – globular clusters: general – galaxies: formation.

## 1 INTRODUCTION

It is intuitive that the gravitational collapse of a set of N masses is directly related to the formation of astrophysical structures such as globular clusters or elliptical galaxies (the presence of gas may complicate the pure gravitational N-body problem for spiral galaxies). From an analytical point of view, this problem is very difficult. When N is much larger than 2, direct approach is intractable, and according to Poincaré results of non-analyticity, exact solutions may be unobtainable. In the context of statistical physics, the situation is more favorable and, in a dissipationless approximation,<fnref n="1" /> leads to the Collisionless Boltzmann Equation (hereafter denoted by CBE)

$$ \frac{\partial f}{\partial t} + \mathbf{p} \cdot \frac{\partial f}{\partial \mathbf{r}} + m \frac{\partial \psi}{\partial \mathbf{r}} \cdot \frac{\partial f}{\partial \mathbf{p}} = 0 \quad (1) $$

where $f=f(r, \mathbf{p}, t)$ and $\psi = \psi(r, t)$ are, respectively, the distribution function of the system with respect to the canonically conjugated $(r, \mathbf{p})$ phase-space variables and the mean field gravitational potential. As noted initially by Hénon (1960), this formalism holds for such systems if and only if we consider N identical point masses equal to m. This problem splits naturally into two related parts: the time-dependent regime and the stationary state. We can reasonably think that these two problems are not completely understood. The transient time-dependent regime was investigated mainly considering self-similar solutions (Lynden-Bell & Eggleton 1980; Henriksen

& Widrow 1995; Blottiau, Bouquet & Chièze 1988; Lancellotti & Kiessling 2001). These studies conclude that power-law solutions can exist for the spatial dependence of the gravitational potential (with various powers). Nevertheless, there is no study which indicates clearly that the time-dependence of the solutions disappears in a few dynamical times, giving a well defined equilibrium-like state. On the other hand, applying Jeans theorem (e.g. Binney & Tremaine 1987, hereafter BT87, p. 220), it is quite easy to find a stationary solution. For example, every positive and integrable function of the mean field energy per mass unit E is a potential equilibrium distribution function for a spherical isotropic system. Several approaches are possible to choose the equilibrium distribution function. Thermodynamics (Violent Relaxation paradigm: Lynden-Bell 1967; Chavanis 2002; Nakamura 2000) indicate that isothermal spheres or polytropic systems are good candidates. Stability analyses can be split into two categories. In the CBE context (see Perez & Aly 1996, for a review), it is well known that spherical systems (with decreasing spatial density) are generally stable except in the case where a large radial anisotropy is present in the velocity space. This is the Radial Orbit Instability, hereafter denoted by ROI (see Perez & Aly 1996; Perez et al. 1996, for a detailed analytic and numeric study of these phenomena) which leads to a bar-like equilibrium state in a few dynamical times. In the context of thermodynamics of self-gravitating systems, in a pioneering work by Antonov (1984), it was shown that an important density contrast leads to the collapse of the core of system (see Chavanis 2003, for details).

In all these studies there is no definitive conclusion, and the choice of the equilibrium distribution remains unclear. Introducing observations and taking into account analytical constraints, several models are possible: chronologically, we can cite (see for example BT87, pp. 223–239) the Plummer model (or other polytropic models), de Vaucouleurs $r^{1/4}$ law, King and isochrone Hénon models or, more recently, the very simple but interesting Hernquist model

<footnote>*E-mail: roy@ensta.fr (FR); perez@ensta.fr (JP)</footnote>

<footnote n="1">¹ The dissipationless hypothesis is widely accepted in the context of gravitational N-body problem because the ratio of the two-body relaxation time over the dynamical time is of the order of N. For a system composed of more than ~10⁴ massive particles, a study during a few hundreds dynamical times can really be considered as dissipationless, the unique source of dissipation being two-body encounters.</footnote>
(Hernquist 1990) for spherical isotropic systems. In the anisotropic case, Ossipkov–Merritt or generalized polytropes can be considered. Finally, for non-spherical systems, there also exists some models reviewed in BT87 (pp. 245–266). Considering this wide variety of possibilities, one can try to make accurate numerical simulations to clarify the situation. Surprisingly, such a program has not been completely carried on. In a pioneering work, van Albada (1982) remarked that the dissipationless collapse of a clumpy cloud of N equal masses could lead to a final stationary state that is quite similar to elliptical galaxies. This kind of study was reconsidered in an important work by Aguilar & Merritt (1990), with more details and a crucial remark concerning the correlation between the final shape (spherical or oblate) and the virial ratio of the initial state. These authors explain this feature invoking ROI. Some more recent studies (Cannizzo & Hollister 1992; Boily, Clarke & Murray 1999; Theis & Spurzem 1999) concentrate on some particularities of the preceding works. Finally, two works (Dantas et al. 2002; Carpintero & Muzzio 1995) develop new ideas considering the influence of the Hubble flow on the collapse. However, the problem is only partially depicted.

The aim of this paper is to analyse the dissipationless collapse of a large set of N-body systems with a very wide variety of ‘realistic’ initial conditions. As we will see, the small number of particles involved, the numerical technique or the specificity of the previous works did not allow their authors to reach a sufficiently precise conclusion. The layout of this paper is as follows. In Section 2 we describe in detail the numerical procedures used in our experiments. Section 3 describes the results we have obtained. These results are then interpreted in Section 4, where some conclusions and perspectives are also proposed.

# 2 NUMERICAL PROCEDURES

## 2.1 Dynamics

The TREECODE used to perform our simulations is a modified version of the Barnes & Hut (1986) TREECODE, parallelized by D. Pfenniger using the MPI library. We implemented some computations of observables and adapted the code to suit our specific problems. The main features of this code are a hierarchical O[N log(N)] algorithm for the gravitational force calculation and a leap-frog algorithm for the time-integration. We introduced an adaptative time-step, based on a very simple physical consideration. The time-step is equal to a fraction $n_{\text{ts}}$ of the instantaneous dynamical time $T_d$, i.e. $\Delta t = T_d/n_{\text{ts}}$. The simulations were run on a Beowulf cluster (25 dual CPU processors, the speeds of which range from 400 MHz to 1 GHz).

## 2.2 Initial conditions

The initial virial ratio is an important parameter in our simulations. The following method was adopted to set the virial ratio to the value $V_{\text{initial}}$. Positions $\mathbf{r}_i$ and velocities $v_i$ are generated. We can then

$$T_d = \frac{\sum_{i=1}^{N} \sqrt{x_i^2 + y_i^2 + z_i^2}}{\sum_{i=1}^{N} \sqrt{vx_i^2 + vy_i^2 + vz_i^2}}.$$

compute

$$V_p = \frac{2K}{U}, \qquad (2)$$

where

$$K = \sum_{i=1}^{N} \frac{1}{2} m_i v_i^2 \qquad (3)$$

and

$$U = -\frac{G}{2} \sum_{i \neq j}^{N} \frac{m_i m_j}{\{\max[(\mathbf{r}_i - \mathbf{r}_j)^2, \epsilon^2]\}^{1/2}} \qquad (4)$$

In this relation $\epsilon$ is a softening parameter whose value is discussed in Section 2.3.2. As the potential energy depends only on the positions, we obtain a system with a virial ratio equal to $V_{\text{initial}}$ simply by multiplying all the particle velocities by the factor $(V_{\text{initial}}/V_p)^{1/2}$. For convenience we define

$$\eta = |V_{\text{initial}}| \times 10^2. \qquad (5)$$

### 2.2.1 Homogeneous density distribution ($C_{\eta}^{n}$)

As we study large N-body systems, we can produce a homogeneous density by generating positions randomly. These systems are also isotropic. We produce the isotropic velocity distribution by generating velocities randomly.

### 2.2.2 Clumpy density distribution ($C_{\eta}^{u}$)

A type of inhomogeneous systems is made of systems with a clumpy density distribution. We first generate n small homogeneous spherical systems with radius $R_g$. The centres of these subsystems are uniformly distributed in the system. The empty space is then filled using a homogeneous density distribution. In the initial state, each clump contains about 1 per cent of the total mass of the system and has a radius which represents 5 per cent of the initial radius of the whole system. These systems are isotropic.

### 2.2.3 Power-law $r^{-\alpha}$ density distribution ($P_{\eta}^{\alpha}$)

We first generate the $\varphi$ and $z$ cylindrical coordinates using two uniform random numbers, $u_1$ and $u_2$:

$$(z, \varphi) = (2u_1 - 1, 2\pi u_2). \qquad (6)$$

Using the inverse transformation method, if

$$r = RF^{-1}(u), \quad \text{with} \quad F(r) = \frac{1}{S} \int_{1}^{r} x^{2-\alpha} \, dx, \qquad (7)$$

where R is the radius of the system, u is a uniform random number, $\iota \ll 1$ and

$$S = \int_{0}^{1} x^{2-\alpha} \, dx, \qquad (8)$$

then the probability density of r is proportional to $r^{2-\alpha}$, and the mass density $\rho$ is proportional to $r^{-\alpha}$. Finally, one gets

$$\mathbf{r} = \begin{bmatrix} r \sqrt{1 - z^2} \cos \varphi \\ r \sqrt{1 - z^2} \sin \varphi \\ rz \end{bmatrix}. \qquad (9)$$

These systems are isotropic.

<footnote n="2">² The fraction $n_{\text{ts}}$ is adapted to the virial parameter $\eta$ and ranges roughly from $n_{\text{ts}} = 300$ when $\eta = 90$ to $n_{\text{ts}} = 5000$ when $\eta = 08$. The dynamical time we used is given by</footnote>
### 2.2.4 Gaussian velocity distribution ($G_\eta^\sigma$)

Most of the systems we use have a uniform velocity distribution. But we have also performed simulations with systems presenting a Gaussian initial velocity distribution. These systems are isotropic, but the x-, y- and z-components of the velocity are generated following a Gaussian distribution. Using a standard method we generate two uniform random numbers, $u_1$ and $u_2$, and set

$$v_i = \sqrt{-2\sigma^2 \ln u_1} \cos(2\pi\sigma^2 u_2) \quad i = x, y, z, \quad (10)$$

where $\sigma$ is the Gaussian standard deviation.

### 2.2.5 Global rotation ($R_\eta^f$)

Some of our initial systems are homogeneous systems with a global rotation around the z-axis. The method we choose to generate such initial conditions is the following. We create a homogeneous and isotropic system (an H-type system). We then compute the average velocity of the particles.

$$\bar{v} = \frac{1}{N} \sum_{i=1}^{N} \|v_i\|. \quad (11)$$

We project the velocities on a spherical referential, and add a fraction of $\bar{v}$ to $v_\phi$ with regard to the position of the particle. We set

$$v_{i,\phi} = v_{i,\phi} + f \frac{\rho_i \bar{v}}{R}, \quad (12)$$

where $f$ is a parameter of the initial condition, $\rho_i$ is the distance from the particle to the z-axis and $R$ the radius of the system. The amount of rotation induced by this method can be evaluated through the ratio

$$\mu = K_{\text{rot}}/K, \quad (13)$$

where $K$ is the total kinetic energy defined above, whereas $K_{\text{rot}}$ is the rotation kinetic energy defined by Navarro & White (1993):

$$K_{\text{rot}} = \frac{1}{2} \sum_{i=1}^{N} m_i \frac{(\mathbf{J}_i \cdot \hat{\mathbf{J}}_{\text{tot}})^2}{[\mathbf{r}_i^2 - (\mathbf{r}_i \cdot \hat{\mathbf{J}}_{\text{tot}})^2]} \quad (14)$$

Above, $\mathbf{J}_i$ is the specific angular momentum of particle i, and $\hat{\mathbf{J}}_{\text{tot}}$ is a unit vector in the direction of the total angular momentum of the system. In order to exclude counter-rotating particles, the sum in equation (14) is actually carried out only over those particles which verify the condition $(\mathbf{J}_i \cdot \hat{\mathbf{J}}_{\text{tot}}) > 0$.

### 2.2.6 Power-law initial mass function ($M_\eta^k$)

Almost all the simulations we made assume particles with equal masses. However, we have created some initial systems with a power-law mass function, such as

$$n(M) = \alpha M^{\beta} \quad (15)$$

The number of particles of mass $M \le m \le M + dM$ is $n(M)$ dM. In some models, the value of $\alpha$ and $\beta$ depends on the range of mass that is considered. We have used several types of mass functions, among them the initial mass function given by Kroupa (2001) ($k = I$), the one given by Salpeter (1955) ($k = II$) and an $M^{-1}$ mass function ($k = III$). In order to generate masses following these functions, we first calculate $\alpha_k$ to produce a continuous function. We can then calculate the number of particles whose mass is between $M$ and $M + dM$. We generate $n(M)$ masses

$$m_i = M + u \, dM \quad 1 \le i \le n(M), \quad (16)$$

where $0 \le u \le 1$ is a uniform random number. In the initial state, these systems have a homogeneous number density, a quasi-homogeneous mass density, and they are isotropic.

### 2.2.7 Nomenclature

We indicate below the whole set of our non-rotating initial conditions.

(i) Homogeneous $H_\eta$ models: $H_{88}$, $H_{79}$, $H_{60}$, $H_{50}$, $H_{40}$, $H_{30}$, $H_{20}$, $H_{15}$ and $H_{10}$.

(ii) Clumpy $C_\eta^\eta$ models: $C_{67}^{20}$, $C_{65}^{20}$, $C_{61}^{20}$, $C_{48}^{20}$, $C_{39}^{20}$, $C_{29}^{20}$, $C_{14}^{20}$, $C_{10}^{20}$, $C_{07}^{20}$ and $C_{10}^{03}$.

(iii) Power-law $P_\eta^\alpha$ models: $P_{50}^{2.0}$, $P_{09}^{2.0}$, $P_{50}^{1.0}$, $P_{50}^{0.5}$, $P_{10}^{1.0}$, $P_{08}^{1.5}$ and $P_{40}^{1.5}$.

(iv) Gaussian velocity profiles $G_\eta^\sigma$ models: $G_{50}^1$, $G_{50}^2$, $G_{50}^3$, $G_{50}^4$, $G_{12}^5$.

(v) Mass spectra $M_\eta^k$ models: $M_{50}^I$, $M_{50}^{II}$, $M_{51}^{III}$, $M_{35}^I$, $M_{25}^{II}$, $M_{15}^{III}$ and $M_{07}^I$.

For all these models we ran the numerical simulations with 30 000 particles (see Section 3.1).

## 2.2 Observables

### 2.3.1 Units

Our units are not the commonly used ones (see Heggie & Mathieu 1986). We did not set the total energy $E$ of the system to $-0.25$ because we wanted to prescribe instead the initial virial ratio $V_{\text{initial}}$, the size $R$ of the system and its mass $M$. We thus have $M = 1$, $R = 10$ and $G = 1$, and values of $V_{\text{initial}}$ and $E$ depending on the simulation. We can link the units we have used with more standard ones. We have chosen the following relationships between our units of length and mass and common astrophysical ones:

$$M = 10^6 M_{\odot} \quad \text{and} \quad R = 10 \, \text{pc}. \quad (17)$$

Our unit of time $u_t$ is given by:

$$1 u_t = \sqrt{\frac{R_c^3 G_s M_s}{R_s^3 G_c M_c}} \approx 4.72 \, 10^{11} \, \text{s} = 1.50 \, 10^4 \, \text{yr}, \quad (18)$$

where the variables $X_s$ are expressed in our simulation units and variables $X_c$ in standard units.

### 2.3.2 Potential softening and energy conservation

The non-conservation of the energy during the numerical evolution has three main sources.

The softening parameter $\epsilon$ introduced in the potential calculus (cf. equation 4) is an obvious one. This parameter introduces a lower cut-off $\Lambda$ in the resolution of length in the simulations. Following Barnes & Hut (1989), structural details up to scale $\Lambda \lesssim 10\epsilon$ are sensitive to the value of $\epsilon$. Moreover, in order to be compatible with the collisionless hypothesis, the softening parameter must be greater than the scale where important collisions can occur. Still following Barnes & Hut (1989), this causes

$$\frac{\epsilon}{10} \gtrsim \frac{G \langle m \rangle}{\langle v^2 \rangle}. \quad (19)$$

In our collapse simulation with $3 \times 10^4$ particles, this results in $\epsilon\eta \gtrsim 2/3$. The discretization of time-integration introduces inevitably another source of energy non-conservation, particularly during the collapse. The force computation also generates errors. The choice
of the opening angle $\Theta$, which governs the accuracy of the force calculation of the TREECODE, is a compromise between speed and accuracy. For all these reasons, we have adopted $\epsilon = 0.1$. This choice imposes $\eta \gtrsim 6$ (for $30 \times 10^4$ particles). This trade-off allowed to perform simulations with less than 1 per cent energy variation without requiring too much computing time. For each of our experiments, the total CPU time ranges between 3 to 24 h for 3000 $u_i$ and $3 \times 10^4$ particles. The total aggregated CPU time of all our collapse experiments is approximately 6 months.

We have tested two other values of the softening parameter ($\epsilon = 0.03$ and $\epsilon = 0.3$) for several typical simulations. These tests did not reveal significant variations of the computed observables.

### 2.3.3 Spatial indicators

As indicators of the geometry of the system, we computed axial ratios, radii containing 10 ($R_{10}$), 50 ($R_{50}$) and 90 per cent ($R_{90}$) of the mass, the density profile $\rho(r)$ and the equilibrium core radius. The axial ratios are computed with the eigenvalues $\lambda_1, \lambda_2$ and $\lambda_3$ of the $(3 \times 3)$ inertia matrix $I$, where $\lambda_3 \le \lambda_2 \le \lambda_1$ and, if the position of the particle $i$ is $r_i = (x_{1,i}; x_{2,i}; x_{3,i})$,

$$ \left\{ \begin{array}{ll} I_{\mu\nu} = -\sum_{i=1}^{N} m_i x_{\mu,i} x_{\nu,i} & \text{for } \mu \neq \nu = 1, 2, 3 \\ I_{\mu\mu} = \sum_{i=1}^{N} m_i (r_i^2 - x_{\mu,i}^2) & \text{for } \mu = 1, 2, 3 \end{array} \right. \qquad (20) $$

The axial ratios $a_1$ and $a_2$ are given by $a_2 = \lambda_1/\lambda_2 \ge 1.0$ and $a_1 = \lambda_3/\lambda_2 \le 1.0$.

The density profile $\rho$, which depends only on the radius $r$, together with the $R_\delta$ ($\Delta= 10, 50, 90$) have a physical meaning only for spherical or nearly spherical systems. For all the spatial indicator computations we have only considered particles whose distance to the centre of mass of the system is less than $6 \times R_{50}$ of the system. This assumption excludes particles which are inevitably 'ejected' during the collapse.<fnref n="3" />

After the collapse a core halo structure forms in the system. In order to measure the radius of the core, we have computed the density radius as defined by Casertano and Hut (see Casertano & Hut 1985). The density radius is a good estimator of the theoretical and observational core radius.

We have also computed the radial density of the system. The density is computed by dividing the system into spherical bins and by calculating the total mass in each bin.

### 2.3.4 Statistical indicators

When the system has reached an equilibrium state, we compute the temperature of the system

$$ T = \frac{2 \langle K \rangle}{3Nk_B}, \qquad (21) $$

where $K$ is the kinetic energy of the system, $k_B$ is the Boltzmann constant (which is set to 1) and the notation $\langle A \rangle$ denotes the mean

value of the observable $A$, defined by

$$ \langle A \rangle = \frac{1}{N} \sum_{i=1}^{N} A_i. \qquad (22) $$

In order to characterize the system in the velocity space we have computed the function

$$ \kappa(r) = \frac{2 \langle v_{i, \text{rad}}^2 \rangle_{r \le r_i < r+dr}}{\langle v_{i, \text{tan}}^2 \rangle_{r \le r_i < r+dr}}, \qquad (23) $$

where $v_{i,\text{rad}}$ is the radial velocity of the $i$th particle, and $v_{i,\text{tan}}$ its tangential velocity. For spherical and isotropic systems ($a_1 \simeq a_2 \simeq 1$ and $\kappa(r) \simeq 1$), we have fitted the density by

(i) a polytropic law

$$ \rho = \rho_0 \psi^\gamma \qquad (24) $$

which corresponds to a distribution function

$$ f(E) \propto E^{\gamma-3/2}; \qquad (25) $$

(ii) an isothermal sphere law

$$ \rho = \rho_1 e^{\psi/s^2} \qquad (26) $$

which corresponds to a distribution function

$$ f(E) = \frac{\rho_1}{(2\pi s^2)^{3/2}} e^{E/s^2} \qquad (27) $$

Using the least-squares method in the $\ln(\rho)-\ln(\psi)$ plane we get [$\gamma$, $\ln(\rho_0)$] and [$s^2$, $\ln(\rho_1)$].

## 3 DESCRIPTION OF THE RESULTS

We have only studied systems with an initial virial ratio corresponding to $\eta \in [7, 88]$. In such systems, the initial velocity dispersion cannot balance the gravitational field. This produces a collapse. After this collapse, in all our simulations the system reaches an equilibrium state characterized by a temporal mean value of the virial ratio equal to $-1$, i.e. $\eta = 100$, and stationary physical observables. These quantities (defined in Section 2) are presented in a table of results in the Appendix of this paper. The following results will be discussed and interpreted in Section 4.

### 3.1 Relevant number of particles

In all previous works on this subject (van Albada 1982; Aguilar & Merritt 1990; Cannizzo & Hollister 1992; Boily et al. 1999), the authors did not really consider the influence of the number of particles on their results. In the first two and more general works, this number is rather small (not more than a few thousands in the largest simulations). The two other studies are more specific and use typically $10^4$ and, in a few reference cases, $2 \times 10^4$ particles. In order to test the influence of the number of particles on the final results, we have computed several physical observables of some collapsing systems with various numbers of particles. The results are presented in Fig. 1. In order to check the influence of $N$ in the whole phase space, we have studied positions and velocities related observables: $a_1, a_2, R_{10}, R_{50}$ and $R_{90}$ and parameters of isothermal and polytropic fit models, namely $\gamma$ and $s^2$. Moreover, in order to be model-independent, we have studied three representative initial conditions: $H_{80}, H_{50}$ and $H_{10}$, i.e. initially hot, warm and cold systems, respectively. The number of particles used in each case ranges from $10^2$ to $10^5$. We can see in Fig. 1 that some observables

<footnote n="3">³ The number of excluded particles ranges from 0 to 30 per cent of the total number of particles, depending mostly on $\eta$. For example, the number of excluded particles is 0 per cent for $H_{80}$, 3 per cent for $C_{67}^{20}$, 5 per cent for $H_{50}$, 22 per cent for $C_{10}^{20}$ and 31 per cent for $H_{10}$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Influence of the number of particles on the physical observables of a collapsing system. Axial ratios are on the top panel, radius containing 50 per cent of the total mass and density radius are on the middle panel and the best $s^2$ and $\gamma$ fit for respectively isothermal and polytropic distribution function are on the bottom panel. All cases are initially homogeneous with $\eta = 10$ (solid line), $\eta = 50$ (dotted line) and $\eta = 80$ (dashed line). The number $N$ of particles used is in units of $10^3$.</figcaption></figure>

are $N$-dependent when $N < 3 \times 10^4$. In particular, $R_{50}$ and $R_d$ present a monotonic variation larger than 50 per cent when $N$ varies from $10^2$ to $10^5$ and the ellipticity of the final state is overestimated for small values of $N$. As a conclusion of this preliminary study, we claim that the relevant number of particles for collapse simulations is $N \ge 3 \times 10^4$. As all simulations have been completed with a total energy loss smaller than 1 per cent, we state moreover that this result is independent of the numerical scheme used (TREECODE or Direct N-Body). As a consequence, the simulations presented hereafter have been performed using $N = 3 \times 10^4$ particles.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Axial ratios of equilibrium states reached from Homogeneous, Clumpy, Gaussian velocity dispersion, Power law and Mass spectrum initial conditions.</figcaption></figure>

## 3.2 Morphological segregation

An important study by Aguilar & Merritt (1990) shows that, in the case of an initial density profile $\rho \propto r^{-1}$, the shape of the virialized state depends on $\eta$: a very small $\eta$ leads to a flattened equilibrium state, when a more quiet collapse produces a spherical one. Our investigations concern a wide range of different initial conditions (homogeneous, clumpy, ...) and show that the influence of $\eta$ depends on the properties of the initial system.<fnref n="4" /> Fig. 2 shows the axial ratios of the equilibrium state reached by our simulations. In fact, only a few simulations produced a final state with an ellipticity greater than E1. Every homogeneous initial condition (i.e. $H_\eta$, $G_\eta^\sigma$ and $M_\eta^k$) resulted in a spherical equilibrium state independently of the values of $\eta$ we tested. Cold clumpy systems have a weakly flattened equilibrium state. The only final systems with an ellipticity significantly greater than E1 are those produced by the collapse of cold $P_\eta^\alpha$.

Previous studies invoked ROI to explain this morphological segregation. However, it seems that ROI requires inhomogeneities near the centre to be triggered.

## 3.3 Characteristic size segregation

In addition to the morphological segregation, presented in the previous section, we discovered a finer phenomenon.

In Figs 3, 4, 5, 6 and 7, we have plotted the mass density of all equilibrium states produced by the collapse of our initial conditions as a function of the ratio $r/R_{50}$. These plots represent the density at the end of our simulations (after about 100 crossing times). These functions do not significantly evolve after the collapse except for $M_{07}^k$. For this special case, a comparative plot is the subject of Fig. 8. All equilibrium states we obtain clearly fall into two categories.

(i) Flat Core Systems. All these systems present a core halo structure, i.e. a large central region with a constant density and a steep envelope. These systems are typically such that $\ln(R_{50}/R_d) < 0.5$ and $\ln(R_{10}/R_d) < -0.05$.

(ii) Small Core Systems. For such systems, the central density is two orders of magnitude larger than for Flat Core systems. There

<footnote n="4">⁴ The particular case of rotating initial conditions is discussed in a special section.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Density profile for $H_\eta$ models plotted in units of $R_{50}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Density profile for $C_\eta^n$ models plotted in units of $R_{50}$. The dashed line corresponds to the $C_{10}^{03}$ model.</figcaption></figure>

is no central plateau and the density falls down regularly outward. These systems are typically such that $\ln(R_{50}/R_d) > 0.7$ and $\ln(R_{10}/R_d) > 0.1$.

The diagram $\ln(R_{10}/R_d)$ versus $\ln(R_{50}/R_d)$ is the subject of Fig. 9. One can see in this figure that each equilibrium state belongs to one or the other family except in a few particular cases. In the Flat Core family we found all $H_\eta$, $G_\eta^\sigma$ and $M_\eta^k$ systems except $M'_\eta$, and two $P_\eta^\alpha$ systems, namely $P_{50}^{0.5}$ and $P_{50}^1$. These systems are all initially homogeneous or slightly inhomogeneous (e.g. $P_{50}^{0.5}$ and $P_{50}^1$ systems). In the Small Core family, we found all $C_\eta^{20}$ and the $P_{50}^{20}$ and $P_{09}^{20}$ systems. These systems are all initially rather very inhomogeneous. Finally, there are five systems in between the two categories: $C_{10}^{03}$, $P_{40}^{1.5}$, $P_{08}^{1.5}$, $P_{10}^1$ and $M'_{07}$. This last model is the only one which migrates from Flat Core set (when $t \simeq 10T_d$) to the edge of the Small Core region (when $t \simeq 100T_d$).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Density profile for $P_\eta^\sigma$ models plotted in units of $R_{50}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Density profile for $G_\eta^\sigma$ models plotted in units of $R_{50}$.</figcaption></figure>

## 3.4 Equilibrium distribution function

In order to compare systems in the whole phase space, we fitted the equilibrium state reached by each system with two distinct isotropic models, e.g. polytropic and isothermal (see equations 24–27 or BT87, pp. 223–232). Fig. 10 shows these two fits for the $P_{50}^{0.5}$ simulation. The technique used for the fit is described in Section 2 of this paper. The result obtained for this special study is as follows. The equilibrium states reached by our initial conditions can be fitted by the two models with a good level of accuracy. As long as $\eta < 70$, the polytropic fit gives a mean value $\gamma = 4.77$ with a standard deviation of $2.48 \ 10^{-1}$. This deviation represents 5.1 per cent of the mean value. This value corresponds typically to the well known Plummer model for which $\gamma = 5$ (see BT87, p. 224 for details). When the collapse is very quiet (typically $\eta > 70$) the polytropic fit is always very good but the value of the index is much larger than Plummer model, e.g. $\gamma = 6.86$ for $H_{79}$ and $\gamma = 7.37$ for $H_{88}$. The corresponding plot is the subject of the Fig. 11. All the data can be found in the Appendix. As we can see on the example plotted in Fig. 10, the isothermal fit is generally not as good as the polytropic
<figure><img src="image_8.png" /><figcaption>Figure 7. Density profile for $M_{\eta}^k$ models plotted in units of $R_{50}$.</figcaption></figure>

one. On the whole set of equilibria, isothermal fits give a mean value $s^2 = 2.5 \times 10^{-2}$ with a standard deviation of $1.6 \times 10^{-2}$ (60 per cent). The corresponding plot is the subject of Fig. 12. All the data can be found in the Appendix.

In fact both isothermal and polytropic fits are reasonable: as long as the model is able to reproduce a core halo structure, the fit is correct. The success of the Plummer model, the density of which is given by

$$ \rho(r) = \frac{3}{4\pi b^3} \left[ 1 + \left( \frac{r}{b} \right)^2 \right]^{-5/2}, $$

can be explained by its ability to fit a wide range of models with various ratio of the core size over the half-mass radius. The adjustment of this ratio is made possible by varying the free parameter $b$. We expect that other core halo models, such as King or Hernquist models, work as well as the Plummer model. As a conclusion of this section, let us say that, as predicted by theory, there is not a single universal model to describe the equilibrium state of isotropic spherical self-gravitating system.

## 3.5 Influence of rotation

We saw in Section 3.2 a source of flattening for self-gravitating equilibrium. Let us now show the influence of initial rotation, which is a natural candidate to produce flattening. The way we have added a global rotation and the significance of our rotation parameters $f$ and $\mu$ are explained in Section 2.2.5. The set of simulations performed for this study contains 31 different elements. The initial virial ratio ranges from $\eta = 10$ to $\eta = 50$, and the rotation parameter from $f = 0$ (i.e. $\mu = 0$) to $f = 20$ (i.e. $\mu = 0.16$ when $\eta = 50$). As a matter of fact, equilibrium states always preserve a rather important part of the initial rotation<fnref n="5" /> and observed elliptical gravitating systems generally possess very small amount of rotation (see e.g. Combes et al. 1995). We thus exclude large values of $f$.

<figure><img src="image_9.png" /><figcaption>Figure 8. Comparison between the evolution of the mass density with respect to time for $C_{10}^{20}$ (top-panel) and $M_{07}^I$. For each case, plotted times are 10,20,30,40,50,75 and 100 $T_{d}$.</figcaption></figure>

Our experiments exhibit two main features (see Fig. 13): on the one hand, rotation produces a flattened equilibrium state only when $f$ exceeds a triggering value (typically $f = f_o \simeq 4$). On the other hand, we have found that for a given value of $\eta$, the flatness of the equilibrium is roughly $f$-independent, provided that $f > f_o$.

## 3.6 Thermodynamical segregation

As we study isolated systems, the total energy $E$ contained in the system is constant during the considered dynamical evolution. This property remains true as long as we consider collisionless evolutions. For gravitational systems, this means that we cannot carry out any simulation of duration larger than a few hundred dynamical times. We have obviously taken this constraint into account in our experiments. All systems which experience a violent relaxation reach an equilibrium state which is stationary in the whole phase space. Spatial behaviour such as morphological segregation produced by ROI was confirmed and further detailed thanks to our study. A new size segregation was found in Section 3.3. Now let us consider another new segregation appearing in the velocity space. Each equilibrium state is associated with a constant temperature $T$, calculated using equation (21). More precisely, we have calculated

<footnote n="5">⁵ We observed that $\mu$ is always smaller in the equilibrium state than in the initial one, typically each rotating systems conserves 65 per cent of the initial $\mu$.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The core size segregation: ln(*R*<sub>10</sub>/*r*<sub>c</sub>) versus ln(*R*<sub>50</sub>/*r*<sub>c</sub>) is plotted for all non-rotating systems.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Polytropic and isothermal fit for the $P_{50}^{0.5}$ simulation.</figcaption></figure>

the temporal mean value<fnref n="6" /> of the temperatures, evaluated every one hundred time units. As we can see in Fig. 14, after the collapse and whatever the nature of the initial system is, the temperature is a very stable parameter.

Fig. 15 shows the *E*–*T* diagram of the set of all non-rotating simulations. It reveals a very interesting feature of post-collapse self-gravitating systems.

On the one hand, the set of systems with a total energy $E > -0.054$ can be linearly fitted in the *E*–*T* plane. We call this set Low Branch 1 (hereafter denoted by LB1, see Fig. 15). On the other hand, the set of the systems with a total energy $E < -0.054$ splits into two families. The first is an exact continuation of LB1. Hence we named it Low Branch 2 (hereafter denoted by LB2). The second can also be linearly fitted, but with a much greater slope (one order of magnitude). We label this family High Branch (hereafter denoted by HB).

In LB1 or LB2, we find every H, G and M system with $\eta > 25$, and every P and C with $n > 10$. In HB, we find $C_{10}^{0.3}$ and every H, G and M system with $\eta < 25$. This segregation thus affects violent

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Best fit of the $\gamma$ parameter for a polytropic model for all non-rotating systems studied. The error bars correspond to the least-square difference between the data and the model.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Best fit of the $s^2$ parameter for an isothermal model for all non-rotating systems studied. The error bar corresponds to the least-square difference between the data and the model.</figcaption></figure>

collapses (cold initial data): on the one hand, when $\eta > 25$ all systems are on LB1, on the other hand for $\eta < 25$, initially homogeneous or quasi-homogeneous (e.g. $C_{10}^{0.3}$) systems reach HB when initially inhomogeneous systems stay on LB2 instead.

## 4 INTERPRETATIONS, CONCLUSIONS AND PERSPECTIVES

Let us now recapitulate the results we have obtained and propose an interpretation:

(i) The equilibrium state produced by the collapse of a set of N gravitating particles is N-independent provided that $N > 3.0 \times 10^4$.

(ii) Without any rotation, the dissipationless collapse of a set of gravitating particles can produce two relatively distinct equilibrium states: (1) if the initial set is homogeneous, the equilibrium has a large core and a steep envelope; or (2) if the initial set contains significant inhomogeneities ($n > 10$ for clumpy systems or $\alpha > 1$

<footnote n="6">⁶ The temporal mean value is computed from the time when the equilibrium is reached until the end of the simulation.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Axial ratios as for different values of $\eta$ as a function of the initial solid rotation parameter $f$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Evolution of the temperature as a function of time.</figcaption></figure>

for power law systems), the equilibrium state has only a Small Core around which the density falls down regularly.

The explanation of this core size segregation is clear: it is associated to the Antonov core-collapse instability occurring when the density contrast between the central and outward regions of a gravitating system is very big. As a matter of fact, if the initial set contains inhomogeneities, these collapse much more rapidly than the whole system<fnref n="7" /> and quickly fall into the central regions. The density contrast then becomes very large and the Antonov instability triggers producing a core-collapse phenomenon. The rest of the system then smoothly collapses around this collapsed core. If there are no inhomogeneities in the initial set, the system collapses as a whole, and the central density grows slowly without reaching the triggering value of the Antonov instability. A large core then forms. Later evolution can also produce core collapse: this is what occurs for our $M_{07}^\alpha$ system (see Fig. 8). This is an initially homogeneous system with a Kroupa mass spectrum which suffers a very strong collapse. As the mass spectrum is not sufficient to bring a lot of mass into the centre of the system quickly enough, Antonov instability does not trigger and a large core forms. As the collapse is very violent,

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Energy–temperature diagram.</figcaption></figure>

an increasing significant part of particles are progressively ‘ejected’ and the core collapse takes progressively place. This is the same phenomenon which is generally invoked to explain the collapsed core of some old globular clusters (e.g. Djorgowski et al. 1986): during its dynamical evolution in the galaxy, some stars are tidally extracted from a globular cluster, and to compensate for this loss the cluster concentrates its core, increasing the density contrast, which sooner or later triggers the Antonov instability.

(iii) Without any rotation, the collapse (violent or quiet) of an homogeneous set of gravitating particles produces an E0 (i.e. spherical) isotropic equilibrium state. There are two possible ways to obtain a flattened equilibrium: (1) introduce a large amount of inhomogeneity near the centre in the initial state, and make a violent collapse ($\eta < 25$); or (2) introduce a sufficient amount ($f > 4$) of rotation in the initial state.

These two ways do not have the same origin and do not produce the same equilibrium state.

In the first case, one can reasonably invoke the ROI: as a matter of fact, as it is explained in a lot of works (see Perez et al. 1996, for example), two features are associated to this phenomenon. First of all, it is an instability which needs an equilibrium state from which it grows. Secondly, it triggers only when a sufficient amount of radial orbits are present. The only non-rotating flattened systems we observed just combine this two conditions: sufficient amount of radial orbits because the collapse is violent and something from which ROI can grow because we have seen in the previous point that inhomogeneities collapse first and quickly join the centre. The fact that cold $P_\eta^\alpha$ systems are more flattened than $C_\eta^\alpha$ ones is in complete accordance with our interpretation: as a matter of fact, by construction, power law systems have an initial central overdensity, whereas clumpy systems create (quickly but not instantaneously) this overdensity bringing the collapsed clumps near the centre. The ROI flattening is oblate ($a_2 \simeq 1$ and $a_1 < 1$).

The rotational flattening is more natural and occurs when the centrifugal force overcomes the gravitational pressure. The rotational flattening is prolate ($a_2 > 1$ and $a_1 \simeq 1$). We notice that initial rotation must be invoked with parsimony to explain the ellipticity of some globular clusters or elliptical galaxies. As a matter of fact, these objects are very weakly rotating systems and our study has shown that the amount of rotation is almost constant during the collapse.

(iv) Spherical equilibria can be suitably fitted by both isothermal and polytropic laws with various indexes. It suggests that any

<footnote n="7">⁷ Because their Jeans length is much more smaller than the one of the whole system.</footnote>
distribution function of the energy exhibiting an adaptable core halo structure (Polytrope, Isothermal, King, Hernquist, ...) can suitably fit the equilibrium produced by the collapse of our initial conditions.

(v) There exists a temperature segregation between equilibrium states. It concerns only initially cold systems (i.e. systems which will suffer a violent collapse): for such systems when $\eta$ decreases, the equilibrium temperature $T$ increases much more for initially homogeneous systems than for initially inhomogeneous systems. On the other hand, whatever their initial homogeneity, quiet collapses are rather all equivalent from the point of view of the equilibrium temperature: $T$ increases in the same way for all systems as $\eta$ decreases. This feature may be the result of the larger influence of the dynamical friction induced by the primordial core on the rapid particles in a violent collapse.

All these properties may be directly confronted to physical data from globular clusters (see the Harris catalogue; Harris 1996) or galaxies observations.

As a matter of fact, in the standard 'bottom-up' scenario of the hierarchical growth of structures, galaxies naturally form from very inhomogeneous medium. Our study then suggests for the equilibrium state of such objects a potential flattening and a collapsed core. This is in very good accordance with the E0 to E7 observed flatness of elliptical galaxies and may be a good explanation for the presence of massive black hole in the centre of galaxies (see Schodel et al. 2002).

On the other hand, globular cluster observations show that these are spherical objects (the few low, flattened clusters all possess a low amount of rotation), and that their core is generally not collapsed (the collapsed core of almost 10 per cent of the galactic globular clusters can be explained by their dynamical evolution through the galaxy). Our study then expect that globular clusters form from homogeneous media.

These conclusions can be tested using the $E-T$ plane. As a matter of fact, we expect that an $E-T$ plane build from galactic data would not present any High Branch whereas the same plane build from globular clusters data would.

## ACKNOWLEDGMENTS

We thank the referee for the relevance of his remarks and suggestions. We thank Joshua E. Barnes, who wrote the original TREECODE. We also particularly thank Daniel Pfenniger for the use of the parallel TREECODE. The simulations were done on the Beowulf cluster at the Laboratoire de Mathématiques Appliquées from the École Nationale Supérieure de Techniques Avancées.

<table><caption>Table A1. Homogeneous initial conditions: $H_\eta$.</caption><thead><tr><th>η</th><th>ΔE<br>(per cent)</th><th>η<br>(end)</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>R<sub>10</sub></th><th>R<sub>50</sub></th><th>R<sub>90</sub></th><th>R<sub>d</sub></th><th>γ</th><th>Σ<sup>2</sup><sub>γ</sub><br>(×10<sup>2</sup>)</th><th>s<sup>2</sup><br>(×10<sup>2</sup>)</th><th>Σ<sup>2</sup><sub>s<sup>2</sup></sub><br>(×10<sup>3</sup>)</th><th>T<br>(×10<sup>11</sup>)</th><th>-E<br>(×10<sup>2</sup>)</th></tr></thead><tbody><tr><td>88</td><td>0.0</td><td>98</td><td>1.02</td><td>0.99</td><td>3.39</td><td>6.53</td><td>10.6</td><td>5.02</td><td>7.37</td><td>-6</td><td>1.12</td><td>-3</td><td>2.30</td><td>3.30</td></tr><tr><td>79</td><td>0.0</td><td>99</td><td>1.00</td><td>1.00</td><td>3.11</td><td>6.01</td><td>9.8</td><td>4.59</td><td>6.86</td><td>-4</td><td>1.37</td><td>-3</td><td>2.20</td><td>3.60</td></tr><tr><td>60</td><td>0.0</td><td>96</td><td>1.01</td><td>0.98</td><td>2.41</td><td>4.73</td><td>11.5</td><td>3.41</td><td>5.05</td><td>-1</td><td>2.04</td><td>-2</td><td>2.89</td><td>4.20</td></tr><tr><td>50</td><td>0.0</td><td>96</td><td>1.01</td><td>0.98</td><td>2.04</td><td>4.09</td><td>15.1</td><td>2.81</td><td>4.73</td><td>-2</td><td>2.49</td><td>-2</td><td>3.04</td><td>4.50</td></tr><tr><td>40</td><td>0.1</td><td>96</td><td>1.01</td><td>0.99</td><td>1.72</td><td>3.51</td><td>25.7</td><td>2.30</td><td>4.72</td><td>-2</td><td>2.79</td><td>-2</td><td>3.39</td><td>4.80</td></tr><tr><td>30</td><td>0.0</td><td>96</td><td>1.02</td><td>0.99</td><td>1.36</td><td>2.88</td><td>253.2</td><td>1.75</td><td>4.66</td><td>-2</td><td>3.27</td><td>-3</td><td>3.74</td><td>5.10</td></tr><tr><td>20</td><td>0.0</td><td>101</td><td>1.01</td><td>0.99</td><td>0.95</td><td>2.22</td><td>874.1</td><td>1.17</td><td>4.68</td><td>-1</td><td>4.03</td><td>-5</td><td>4.78</td><td>5.39</td></tr><tr><td>15</td><td>0.0</td><td>108</td><td>1.01</td><td>0.99</td><td>0.74</td><td>1.89</td><td>1143.0</td><td>0.88</td><td>4.66</td><td>-1</td><td>4.73</td><td>-6</td><td>5.83</td><td>5.53</td></tr><tr><td>10</td><td>1.4</td><td>120</td><td>1.02</td><td>0.98</td><td>0.52</td><td>1.59</td><td>1448.0</td><td>0.60</td><td>4.59</td><td>-1</td><td>5.76</td><td>-9</td><td>7.71</td><td>5.66</td></tr></tbody></table>

## REFERENCES

Aguilar L., Merritt D., 1990, ApJ, 354, 33

Antonov V. A., 1984, in Goodman J., Hut P., eds, Proc. IAU Symp. 113, Dynamics of Star Clusters. Reidel, Dordrecht, p. 525 (Original Russian version: 1962, Vest. Len. Univ., 7, 135)

Barnes J., Hut P., 1986, Nat, 324, 446

Barnes J., Hut P., 1989, ApJS, 70, 389

Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton University Press Princeton

Blottiau P., Bouquet S., Chièze J.-P., 1988, A&A, 207, 24

Boily C. M., Clarke C. J., Murray S. D., 1999, MNRAS, 302, 399

Cannizzo J. K., Hollister T. C., 1992, ApJ, 400, 58

Carpintero D. D., Muzzio J. C., 1995, ApJ, 440, 5

Casertano S., Hut P., 1985, ApJ, 298, 80

Chavanis P.-H., 2002, Phys. Rev. E, 66, 036105

Chavanis P.-H., 2003, A&A, 401, 15

Combes F., Boissé P., Mazure A., Blanchard A., 1995, Galaxies and Cosmology. Springer-Verlag Heidelberg

Dantas C. C., Capelato H. V., de Carvalho R. R., Ribiero A. L. B., 2002, A&A, 384, 772

Djorgovski S., King I. R., Vuosalo C., Oren A., Penner H., 1986, in Hearnshaw J. B., Cottrell P.L., eds, Proc. IAU Symp. 118, Instrumentation and Research Programs for Small Telescopes. Reidel, Dordrecht, p. 281

Harris W. E., 1996, AJ, 112, 1487

Heggie D. C., Mathieu R. D., 1986, in Hut P., McMillan S. L. W., eds, The Use of Supercomputers in Stellar Dynamics. Springer-Verlag, Heidelberg

Hénon- M., 1960, An. Astro, 23, 474

Henriksen R. N., Widrow L. M., 1995, MNRAS, 276, 679

Hernquist L., 1990, ApJ, 356, 359

Kroupa P., 2001, MNRAS, 322, 231

Lancellotti C., Kiessling M., 2001, ApJ, 549, L93

Lynden-Bell D., 1967, MNRAS, 136, 601

Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483

Nakamura T.K., 2000, ApJ, 531, 739

Navarro J. E., White S. D. M., 1993, MNRAS, 265, 271

Perez J., Aly J.-J., 1996, MNRAS, 280, 689

Perez J., Alimi J.-M., Aly J.-J., Scholl H., 1996, MNRAS, 280, 700

Salpeter E. E., 1955, ApJ, 121, 161

Schodel R. et al., 2002, Nat, 419, 694

Theis C., Spurzem R., 1999, A&A, 341, 361

van Albada T. S., 1982, MNRAS, 201, 939

## APPENDIX A: TABLES OF RESULTS

The numerical values of the observables we have computed for all our equilibrium states are presented in the following tables. These observables are defined in Section 2, except $\sum_{\gamma}^2$ and $\sum_{s^2}$ which are the residuals of the least-squares method used to determine the $\gamma$ and $s^2$ observables.
<table><caption>**Table A2.** Clumpy Initial Condition: $C_{\eta}^{n}$.</caption><thead><tr><th>η</th><th>n</th><th>ΔE<br>(per cent)</th><th>η<br>(end)</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>R<sub>10</sub></th><th>R<sub>50</sub></th><th>R<sub>90</sub></th><th>R<sub>d</sub></th><th>γ</th><th>Σ<sup>2</sup><sub>γ</sub><br>(×10<sup>2</sup>)</th><th>s<sup>2</sup><br>(×10<sup>2</sup>)</th><th>Σ<sub>s2</sub><br>(×10<sup>3</sup>)</th><th>T<br>(×10<sup>11</sup>)</th><th>-E<br>(×10<sup>2</sup>)</th></tr></thead><tbody><tr><td>10</td><td>3</td><td>0.1</td><td>97</td><td>1.03</td><td>0.96</td><td>0.55</td><td>1.85</td><td>1241.0</td><td>0.58</td><td>4.61</td><td>-1</td><td>5.21</td><td>-11</td><td>6.99</td><td>5.82</td></tr><tr><td>67</td><td>20</td><td>0.8</td><td>95</td><td>1.01</td><td>0.94</td><td>0.93</td><td>6.14</td><td>16.1</td><td>0.66</td><td>5.00</td><td>-1</td><td>1.98</td><td>-10</td><td>2.86</td><td>3.96</td></tr><tr><td>65</td><td>20</td><td>0.9</td><td>95</td><td>1.02</td><td>0.96</td><td>0.93</td><td>5.63</td><td>13.5</td><td>0.63</td><td>5.27</td><td>-2</td><td>2.00</td><td>-13</td><td>3.07</td><td>4.29</td></tr><tr><td>61</td><td>20</td><td>0.6</td><td>95</td><td>1.05</td><td>0.98</td><td>0.86</td><td>5.15</td><td>12.3</td><td>0.58</td><td>5.43</td><td>-3</td><td>2.09</td><td>-11</td><td>3.32</td><td>4.64</td></tr><tr><td>48</td><td>20</td><td>0.6</td><td>94</td><td>1.03</td><td>0.99</td><td>0.81</td><td>4.24</td><td>11.8</td><td>0.59</td><td>5.15</td><td>-3</td><td>2.61</td><td>-13</td><td>3.79</td><td>5.31</td></tr><tr><td>39</td><td>20</td><td>0.5</td><td>94</td><td>1.04</td><td>0.99</td><td>0.76</td><td>3.84</td><td>12.8</td><td>0.40</td><td>4.72</td><td>-1</td><td>3.20</td><td>-10</td><td>4.06</td><td>5.65</td></tr><tr><td>29</td><td>20</td><td>1.2</td><td>94</td><td>1.04</td><td>0.99</td><td>0.72</td><td>3.42</td><td>15.4</td><td>0.56</td><td>4.42</td><td>-1</td><td>3.81</td><td>-6</td><td>4.24</td><td>5.97</td></tr><tr><td>14</td><td>20</td><td>0.3</td><td>93</td><td>1.09</td><td>0.98</td><td>0.64</td><td>2.72</td><td>39.5</td><td>0.48</td><td>4.56</td><td>-1</td><td>4.27</td><td>-7</td><td>4.66</td><td>6.50</td></tr><tr><td>10</td><td>20</td><td>1.4</td><td>97</td><td>1.13</td><td>0.99</td><td>0.61</td><td>2.54</td><td>345.5</td><td>0.54</td><td>4.70</td><td>-9</td><td>4.25</td><td>-15</td><td>4.98</td><td>6.66</td></tr><tr><td>7</td><td>20</td><td>0.3</td><td>94</td><td>1.14</td><td>0.92</td><td>0.57</td><td>2.44</td><td>224.6</td><td>0.45</td><td>4.74</td><td>-1</td><td>4.41</td><td>-12</td><td>5.36</td><td>6.76</td></tr></tbody></table>

<table><caption>**Table A3.** Power-law initial conditions: $P_{\eta}^{\alpha}$.</caption><thead><tr><th>η</th><th>α</th><th>ΔE<br>(per cent)</th><th>η<br>(end)</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>R<sub>10</sub></th><th>R<sub>50</sub></th><th>R<sub>90</sub></th><th>R<sub>d</sub></th><th>γ</th><th>Σ<sup>2</sup><sub>γ</sub><br>(×10<sup>2</sup>)</th><th>s<sup>2</sup><br>(×10<sup>2</sup>)</th><th>Σ<sub>s2</sub><br>(×10<sup>3</sup>)</th><th>T<br>(×10<sup>11</sup>)</th><th>-E<br>(×10<sup>2</sup>)</th></tr></thead><tbody><tr><td>50</td><td>0.5</td><td>0.0</td><td>95</td><td>1.01</td><td>0.99</td><td>1.84</td><td>3.92</td><td>13.5</td><td>2.53</td><td>4.66</td><td>-1</td><td>2.65</td><td>-2</td><td>3.54</td><td>4.69</td></tr><tr><td>50</td><td>1</td><td>0.0</td><td>94</td><td>1.01</td><td>0.99</td><td>1.56</td><td>3.77</td><td>12.1</td><td>2.01</td><td>4.77</td><td>-6</td><td>2.78</td><td>-3</td><td>3.69</td><td>5.00</td></tr><tr><td>10</td><td>1</td><td>0.1</td><td>96</td><td>1.00</td><td>0.80</td><td>0.69</td><td>2.71</td><td>382.2</td><td>0.70</td><td>4.61</td><td>-8</td><td>4.05</td><td>-8</td><td>4.85</td><td>6.32</td></tr><tr><td>8</td><td>1.5</td><td>0.1</td><td>96</td><td>1.01</td><td>0.71</td><td>0.62</td><td>2.63</td><td>25.1</td><td>0.61</td><td>4.63</td><td>-7</td><td>4.42</td><td>-9</td><td>5.52</td><td>7.18</td></tr><tr><td>50</td><td>2</td><td>1.7</td><td>93</td><td>1.02</td><td>0.99</td><td>0.53</td><td>3.20</td><td>9.2</td><td>0.34</td><td>5.30</td><td>-6</td><td>3.35</td><td>-9</td><td>5.30</td><td>7.32</td></tr><tr><td>40</td><td>1.5</td><td>0.1</td><td>96</td><td>1.00</td><td>0.99</td><td>0.97</td><td>3.31</td><td>11.0</td><td>1.03</td><td>4.71</td><td>-8</td><td>3.44</td><td>-6</td><td>4.42</td><td>5.99</td></tr><tr><td>9</td><td>2</td><td>1.6</td><td>96</td><td>1.01</td><td>0.78</td><td>0.38</td><td>2.51</td><td>10.6</td><td>0.18</td><td>4.68</td><td>-10</td><td>5.21</td><td>-20</td><td>6.73</td><td>9.31</td></tr></tbody></table>

<table><caption>**Table A4.** Mass spectrum initial conditions: $M_{\eta}^{k}$.</caption><thead><tr><th>η</th><th>k</th><th>ΔE<br>(per cent)</th><th>η<br>(end)</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>R<sub>10</sub></th><th>R<sub>50</sub></th><th>R<sub>90</sub></th><th>R<sub>d</sub></th><th>γ</th><th>Σ<sup>2</sup><sub>γ</sub><br>(×10<sup>2</sup>)</th><th>s<sup>2</sup><br>(×10<sup>2</sup>)</th><th>Σ<sub>s2</sub><br>(×10<sup>3</sup>)</th><th>T<br>(×10<sup>11</sup>)</th><th>-E<br>(×10<sup>2</sup>)</th></tr></thead><tbody><tr><td>7</td><td>Krou</td><td>5.0</td><td>132</td><td>1.02</td><td>0.99</td><td>0.25</td><td>2.04</td><td>1721.0</td><td>0.37</td><td>4.26</td><td>-1</td><td>5.16</td><td>-15</td><td>9.97</td><td>5.62</td></tr><tr><td>15</td><td>1/M</td><td>0.6</td><td>101</td><td>1.01</td><td>0.98</td><td>0.68</td><td>2.04</td><td>366.8</td><td>0.97</td><td>4.65</td><td>-1</td><td>4.35</td><td>-10</td><td>5.72</td><td>5.53</td></tr><tr><td>25</td><td>Salp</td><td>0.4</td><td>99</td><td>1.01</td><td>0.98</td><td>1.18</td><td>2.51</td><td>225.5</td><td>1.54</td><td>4.66</td><td>-2</td><td>3.59</td><td>-3</td><td>4.08</td><td>5.23</td></tr><tr><td>35</td><td>Krou</td><td>0.2</td><td>98</td><td>1.01</td><td>0.99</td><td>1.55</td><td>3.18</td><td>39.8</td><td>2.09</td><td>4.66</td><td>-2</td><td>2.97</td><td>-3</td><td>3.51</td><td>4.95</td></tr><tr><td>51</td><td>1/M</td><td>0.2</td><td>95</td><td>1.01</td><td>0.98</td><td>1.79</td><td>4.19</td><td>14.2</td><td>2.83</td><td>4.67</td><td>-9</td><td>2.39</td><td>-4</td><td>3.15</td><td>4.48</td></tr><tr><td>50</td><td>Krou</td><td>0.1</td><td>96</td><td>1.02</td><td>0.98</td><td>1.93</td><td>4.09</td><td>15.1</td><td>2.87</td><td>4.60</td><td>-1</td><td>2.45</td><td>-3</td><td>3.19</td><td>4.50</td></tr><tr><td>50</td><td>Salp</td><td>0.1</td><td>96</td><td>1.01</td><td>0.98</td><td>2.03</td><td>4.13</td><td>15.1</td><td>2.87</td><td>4.70</td><td>-2</td><td>2.45</td><td>-2</td><td>3.08</td><td>4.49</td></tr></tbody></table>

<table><caption>**Table A5.** Gaussian velocity dispersion initial conditions: $G_{\eta}^{\alpha}$.</caption><thead><tr><th>η</th><th>σ</th><th>ΔE<br>(per cent)</th><th>η<br>(end)</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>R<sub>10</sub></th><th>R<sub>50</sub></th><th>R<sub>90</sub></th><th>R<sub>d</sub></th><th>γ</th><th>Σ<sup>2</sup><sub>γ</sub><br>(×10<sup>2</sup>)</th><th>s<sup>2</sup><br>(×10<sup>2</sup>)</th><th>Σ<sub>s2</sub><br>(×10<sup>3</sup>)</th><th>T<br>(×10<sup>11</sup>)</th><th>-E<br>(×10<sup>2</sup>)</th></tr></thead><tbody><tr><td>48</td><td>G1</td><td>0.0</td><td>95</td><td>1.00</td><td>1.00</td><td>1.72</td><td>3.98</td><td>14.6</td><td>2.23</td><td>4.66</td><td>-5</td><td>2.64</td><td>-3</td><td>3.13</td><td>4.50</td></tr><tr><td>49</td><td>G2</td><td>0.0</td><td>95</td><td>1.02</td><td>0.99</td><td>1.74</td><td>4.02</td><td>14.4</td><td>2.28</td><td>4.65</td><td>-6</td><td>2.61</td><td>-3</td><td>3.13</td><td>4.50</td></tr><tr><td>50</td><td>G3</td><td>0.0</td><td>95</td><td>1.00</td><td>1.00</td><td>1.90</td><td>4.08</td><td>14.1</td><td>2.56</td><td>4.72</td><td>-8</td><td>2.52</td><td>-2</td><td>3.09</td><td>4.50</td></tr><tr><td>12</td><td>G4</td><td>0.4</td><td>118</td><td>1.03</td><td>0.98</td><td>0.56</td><td>1.77</td><td>1312.0</td><td>0.64</td><td>4.56</td><td>-1</td><td>5.33</td><td>-10</td><td>6.08</td><td>5.63</td></tr><tr><td>50</td><td>G5</td><td>0.0</td><td>96</td><td>1.01</td><td>0.99</td><td>1.98</td><td>4.11</td><td>14.8</td><td>2.72</td><td>4.71</td><td>-1</td><td>2.50</td><td>-2</td><td>3.19</td><td>4.50</td></tr></tbody></table>

This paper has been typeset from a TeX/Latex file prepared by the author.