# Stability of spherical stellar systems II: Numerical results

Jérôme Perez, J.-M. Alimi, J.-J. Aly, H. Scholl

► To cite this version:

Jérôme Perez, J.-M. Alimi, J.-J. Aly, H. Scholl. Stability of spherical stellar systems II: Numerical results. Monthly Notices of the Royal Astronomical Society, 1996, 280 (3), pp.700-710. 10.1093/mn-ras/280.3.700. hal-01010753

HAL Id: hal-01010753

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

Submitted on 21 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
# Stability of spherical stellar systems – II. Numerical results

Jérôme Perez,¹,² Jean-Michel Alimi,³ Jean-Jacques Aly¹ and Hans Scholl⁴

¹CEA, DSM/DAPNIA, Service d'Astrophysique (CNRS URA 2052), CE Saclay, 91191 Gif-sur-Yvette Cedex, France

²ETCA/CREA, 16 bis avenue Prieur de la côte d'or, 94114 Arcueil Cedex, France

³Laboratoire d'Astronomie Extragalactique et de Cosmologie, Observatoire de Meudon, 5, Place Jules Janssen, 92195 Meudon, France

⁴Département Cassini, Observatoire de la Côte d'azur, BP 229, 06304 Nice Cedex 4, France

Accepted 1995 November 22. Received 1995 October 30; in original form 1995 July 12

## ABSTRACT

We have performed a series of high-resolution *N*-body experiments on a connection machine CM-5 in order to study the stability of collisionless self-gravitating spherical systems. We interpret our results in the framework of symplectic mechanics, which provides the definition of a new class of particular perturbations: the preserving perturbations, which are a generalization of the radial ones. Using models defined by the Ossipkov–Merritt algorithm, we show that the stability of a spherical anisotropic system is directly related to the preserving or non-preserving nature of the perturbations acting on the system. We then generalize our results to all spherical systems.

Since the ‘isotropic component’ of the linear variation of the distribution function cannot be used to predict the stability or instability of a spherical system, we propose a more useful stability parameter which is derived from the ‘anisotropic’ component of the linear variation.

**Key words:** instabilities – celestial mechanics, stellar dynamics.

## 1 INTRODUCTION

A steady state of a collisionless self-gravitating system – in which the distribution function (DF) is determined only by the isolating integrals of motion (Jeans’ theorem) – may be subject to some form of dynamical instability. After being slightly perturbed, it thus suffers a rapid evolution, which may lead to significant changes after a few crossing times. It is well known, for instance, that spherical anisotropic systems with components moving mainly on radial orbits are unstable. The physical origin of this instability has been extensively investigated analytically (Antonov 1973; Palmer & Papaloizou 1986) and numerically (Henon 1973; Merritt & Aguilar 1985; Barnes et al. 1986; de Zeeuw & Franx 1991; Dejonghe & Merritt 1993). In particular, Fridman & Polyachenko (1984) reached the conclusion that a spherical system is unstable once a parameter $\xi$ – actually equal to the ratio between the radial and transversal kinetic energies of the system – exceeds some critical value ($\xi > 1.7 \pm 0.2$). However, this criterion is only an empirical suggestion and, unfortunately, subsequent simulations (Palmer 1994a, 1994b) have shown its invalidity in many cases. It thus appears at present that there are no simple and general

stability criteria in which a single parameter intervenes depending on the state of the gravitational system.

In this paper, we revisit this problem numerically, and propose a stability parameter for all spherical anisotropic systems. Our work rests on some analytical results (Perez 1995, Perez & Aly 1996a, hereafter Paper I) which have been recently obtained by working in the framework of the symplectic formulation of Vlasov-Poisson (VP) equations (for an introduction to the symplectic formalism, see e.g., Arnold 1978). In this approach – first introduced in the gravitational context by Bartholomew and by Kandrup (Bartholomew 1971; Kandrup 1990, 1991a), and thus developed by Perez and by Aly (Perez 1995; Perez & Aly 1996b, in preparation) – the Vlasov-Poisson system is rewritten in algebraic form, from which one can obtain a formal expression giving the value at any time of an arbitrary function of the DF (whose initial value is assumed to be given). For a steady-state system submitted to a ‘symplectic’ perturbation induced by a ‘generator’ $g$, it is then possible in particular to write a Taylor expansion of its energy in terms of $g$, and to deduce linear stability criteria based on the consideration of the first non-vanishing term in that expansion. Using this general result, some form of
stability of the anisotropic spheres against the so-called 'preserving perturbations' (whose class contains as a small subclass all the radial perturbations) was demonstrated in Paper I, thus providing a generalization of an earlier result (see below). Clearly, this is not sufficient to give a complete description of the stability properties of a spherical system (except in the two limiting cases of a stable isotropic system and of an unstable one with purely radial orbits), and it was conjectured in Paper I that the stability in the intermediate cases is directly related to the preserving or non-preserving nature of the perturbations. The numerical simulations presented here seem to confirm this guess.

It may be worth recalling that two other analytical approaches have been developed in the past 30 years for studying the stability of collisionless self-gravitating systems. The first rests on an angle-action description and on a decomposition of the gravitational potential into normal modes, and it applies to a large class of systems (Fridman & Polyachenko 1984; Palmer & Papaloizou 1986; Goodman 1988; Polyachenko 1992). However, the calculations and the derived stability criteria are often tedious, and difficult to apply. The second method, inspired from plasma physics (Laval, Mercier & Pellat 1965), is related to the variational energy principle, and it has been adopted by Kulsrud & Mark (1970) to gravitational systems. This method has yielded interesting criteria in two situations. On one hand, the stability of radially perturbed isotropic or anisotropic spheres has been obtained – the most general result here being that found by Doremus and his collaborators (Doremus, Feix & Baumann 1971), whose analytical arguments were thus considerably simplified by Sygnet et al. (1984) and Kandrup & Sygnet (1985). On the other hand, the stability of isotropic spheres against non-radial perturbations has also been established. The proof here appeals to the well-known Antonov-Lebovitz theorem (Binney & Tremaine 1987) on the stability of gaseous spheres (Lynden-Bell & Sannit 1969; Antonov 1973), a simplified derivation of which has been reported most recently by Aly & Perez (1992). A strong limitation of the energy principle arises from the necessity to have the 'operator of the small motions' being Hermitian – which is not the case, e.g., for a non-radially perturbed anisotropic sphere.

The paper is arranged as follows. In Section 2, we outline the main analytical results obtained in the framework of symplectic mechanics concerning the stability of collisionless self-gravitating spherical systems. In Section 3 we present the numerical methods for simulating the initial conditions and the dynamical evolution of such systems. In Section 4, we give a physical interpretation of their (in)stability in terms of the 'symplectic quantities' introduced in Section 2. Finally (Section 5), we deduce from the numerical results two laws concerning the stability of any anisotropic collisionless self-gravitating spherical systems.

## 2 THEORETICAL CONTEXT AND ANALYTICAL RESULTS

If $f$ is the DF of the system in the $\{\xi := (x, v)\}$ phase space, $\psi(x)$ the gravitational potential and $G$ the gravitational constant, the evolution of the system (during a period compatible with the collisionless assumption) is given by the Vlasov-Poisson system

$$ 
\left\{ 
\begin{aligned}
\frac{\partial f}{\partial t} + v \cdot \nabla_x f - \nabla_x \psi \cdot \nabla_v f &= 0, \\
\psi(x,t) &= -4\pi G \int \frac{f(x', v', t)}{|x - x'|} d\xi',
\end{aligned}
\right. 
\quad (1) 
$$

where $d\xi := dx dv$ denotes an infinitesimal volume of phase space.

It is well known that Vlasov's equation (1) can be written in the standard Poisson Bracket form

$$ \frac{\partial f}{\partial t} = - \nabla_v E \cdot \nabla_x f + \nabla_x E \cdot \nabla_v f := \{E, f\}, \quad (2) $$

where $E = v^2/2 + \psi$ is the one-particle energy of the system. $E$ can be viewed as the generator of the canonical transformation describing the motion, and its conjugated quantity is the time $t$. More generally, the evolution equation of $f$ in any transformation defined by a generator $g_1$ and its conjugated quantity $\lambda$ can be written

$$ \frac{\partial f}{\partial \lambda} = \{g_1, f\}. \quad (3) $$

The previous equation can also be generalized to yield an evolution equation for any functional of $f$. The time evolution is then generated by the total physical energy $H$ associated with $f$

$$ H[f] = \int \frac{|v|^2}{2} f(x, v, t) \, d\xi - \frac{G}{2} \int d\xi \\ \qquad \qquad \qquad \qquad \int d\xi' \frac{f(x, v, t) f(x', v', t)}{|x - x'|} . \quad (4) $$

$H$ is a functional of $f$ and its functional derivative is

$$ \frac{\delta H[f]}{\delta f} = E. \quad (5) $$

For any $F[f]$ we have

$$ \dot{F}[f] := \frac{\partial}{\partial t} F[f] = \int \frac{\delta F}{\delta f} f \, d\xi. \quad (6) $$

Inserting the bracket form of Vlasov equations (2) and (5) into equation (6), we obtain, after integration by parts, the functional form of Vlasov's equation describing the time evolution of $F[f]$ (Morrisson 1980)

$$ \dot{F}[f] = \int f \left\{ \frac{\delta F}{\delta f}, \frac{\delta H}{\delta f} \right\} d\xi := [F, H][f]. \quad (7) $$

$[\cdot, \cdot]$ is a Lie bracket.

More generally, the functional evolution equation of $F[f]$ for any transformation defined by a functional generator $G_1(f) := \int g_1 f d\xi$ can be written

$$ \forall F \quad \frac{\partial}{\partial \lambda} F[f] = \int f \left\{ \frac{\delta F}{\delta f}, \frac{\delta G_1}{\delta f} \right\} d\xi = [F, G_1][f]. \quad (8) $$
The solution to (8) can be written in exponential form as

$$ \forall F \quad F[f] = (e^{\lambda_1 f} F)[f_{\lambda_0}], \quad (9) $$

or more explicitly as a Taylor expansion;

$$ \begin{aligned} \forall F \quad F[f] = &F[f_{\lambda_0}] - [G_1, F][f_{\lambda_0}] (\lambda - \lambda_0) \\ &+ \frac{[G_1, [G_1, F]] [f_{\lambda_0}]}{2!} (\lambda - \lambda_0)^2 \\ &- \frac{[G_1, [G_1, [G_1, F]]][f_{\lambda_0}]}{3!} (\lambda - \lambda_0)^3 + \dots, \end{aligned} \quad (10) $$

where $f_{\lambda_0}$ denotes the value of the DF at the point $\lambda = \lambda_0$.

In the particular case $F=H$, equation (10) provides a development of the total energy of the system. When $f_{\lambda_0}$ is an equilibrium DF $f_0$, the first non-vanishing term of this development gives a linear stability criterion for the system, independently of its geometry.

$$ \begin{aligned} H^{(1)}[f_0] = &-[G_1, H][f_0] = -\int f_0 \left\{ \frac{\delta G_1}{\delta f}, \frac{\delta H}{\delta f} \right\} d\xi \\ &= \int g_1 \{E, f_0\} d\xi \end{aligned} \quad (11) $$

is clearly vanishing ($\{E, f_0\} = \dot{f}_0 = 0$), while

$$ H^{(2)}[f_0] = \frac{[G_1, [G_1, F]] [f_0]}{2!} \quad (12) $$

can be written after some straightforward algebra (Perez 1995; Paper I)

$$ \begin{aligned} H^{(2)}[f_0] = &-\frac{1}{2} \int \{g_1, E\} \{g_1, f_0\} d\xi \\ &-\frac{1}{2} \int \frac{\{g_1, f_0\} \cdot \{g'_1, f'_0\}}{|x-x'|} d\xi d\xi'. \end{aligned} \quad (13) $$

The positiveness of $H^{(2)}[f_0]$ determines the stability of the system against perturbations generated by some $g_1$. Using a less general method, this result had been previously obtained by Bartholomew (1971) and Kandrup (1990, 1991a). It is important to notice that this criteria on the positiveness of $H^{(2)}[f_0]$ is only a stability criterion; it tells us nothing about the linear instability of the system. Kandrup (1991b) showed that if $H^{(2)}[f_0]$ is negative the system develops a secular instability in the presence of dissipation.

In the particular case of spherical self-gravitating systems, where the equilibrium DF depends only on the one-particle energy $E$ and squared total angular momentum $L^2$, the linear variation of the DF, $\delta f: = \{g_1, f_0\}$, which appears in the second-order energy variation (equation 13), can be written as

$$ \delta f = \frac{\partial f_0}{\partial E} \{g_1, E\} + \frac{\partial f_0}{\partial L^2} \{g_1, L^2\}. \quad (14) $$

When one considers a preserving perturbation defined by $\{g_1, L^2\} = 0$, an anisotropic spherical system behaves as a stable isotropic spherical system,<fnref n="1" /> with $H^2[f_0]$ positive (Perez 1995; Paper I). On the contrary, systems with components evolving only on radial orbits are unstable. Such systems can be described by an equilibrium distribution function $f_0(E, L^2) = v(E)\delta(L^2)$, where $v(E)$ is an arbitrary monotonic function of $E$, and $\delta(L^2)$ denotes a Dirac distribution for $L^2$. In this case, it is then easy to show that all preserving physical perturbations are vanishing. The two previous opposite cases (stable isotropic spherical systems and unstable radial orbit anisotropic spherical systems) suggest that the stability of a spherical self-gravitating system is directly related to the preserving or non-preserving nature of the perturbations acting on the system. In order to confirm this conjecture in all the intermediate cases, which cannot be studied analytically, we have performed numerical simulations. Moreover, under the light of previous analytical work, we will propose also a stability parameter.

# 3 NUMERICAL METHODS

## 3.1 Initial conditions

We use the Ossipkov-Merritt algorithm (Ossipkov 1979; Merritt 1985a,b; Binney & Tremaine 1987) which is well adapted for generating anisotropic self-gravitating spherical systems with various physical properties (see Appendix).

This algorithm has the effect of deforming the isotropic density $\rho_{iso}(r)$ deduced from a given isotropic gravitational potential $\psi_{iso}(r)$ in the following way ($r = |x|$):

$$ \rho_{\text{ani}}(r) := \left(1 + \frac{r^2}{r_a^2}\right) \rho_{\text{iso}}(r). \quad (15) $$

The anisotropic radius $r_a$ controls the deformation. Using the Abel inversion technique, this procedure allows us to generate an anisotropic DF that depends on both $E$ and $L^2$ through the variable $Q := E + L^2/2r_a^2$:

$$ f_0(Q) = \frac{\sqrt{2}}{4\pi^2} \frac{d}{dQ} \int_Q^0 \frac{d\psi_{iso}}{\sqrt{\psi_{iso} - Q}} \frac{d\rho_{\text{ani}}}{d\psi_{iso}}. \quad (16) $$

The velocity anisotropy at any radius $r$ is given by

$$ \frac{\sigma_r^2}{\sigma_t^2} := \frac{\langle v_r^2 \rangle}{1/2 \langle v_t^2 \rangle} = 1 + \frac{r^2}{r_a^2}. \quad (17) $$

the model is isotropic near the centre and becomes more and more anisotropic outwards.

For the isotropic potential, we have chosen here the polytropic model. Thus $\psi_{iso}$ is a solution to the Lame-Emden differential equation

$$ \frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{d\psi}{dr} \right) = (-1)^n 4\pi G c_n \psi^n, \quad c_n = \frac{(2\pi)^{3/2} \Gamma(n-1/2)}{\Gamma(n+1)}, \quad (18) $$

<footnote n="1">¹We consider all along this paper only systems with a DF which admits a monotonous decreasing dependence with respect to all the isolating integrals of motion ($\rho f_0/\partial E < 0$ and $\partial f_0/\partial L^2 < 0$).</footnote>
with $\psi(0)$ being taken as a free parameter. We have set $\psi(0)$ equal to $-1$. Because of the spherical symmetry, we have $d\psi/dr=0$ at $r=0$.

The system admits a finite density and a finite mass provided the polytropic index falls in the $1/2 < n \le 5$ (Chan-

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Phase-space DF for $n=4.5$ with $r_a=0.75$ (dotted line), 1, 2, 5, 10 and 100 (solid lines), respectively. As $r_a$ decreases, $f$ becomes progressively smaller at low energy ($Q \to 1$).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Relative gravitational potential for $n=4.5$ with $r_a=0.75$ (dotted line), 1, 2, 5, 10 and 100 (solid lines), respectively. The radius $R$ at which $\psi$ vanishes increases with the anisotropy ($r_a \to 0$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Radial density profile for a system with $n=4$ and $r_a=2$. The solid line is derived from equations (16) and (18), the symbols correspond to the radial density profile deduced from particles.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The radial dependence of the ratio between the radial velocity dispersion and tangential velocity dispersion for $n=4$ and $r_a=2$. The plain line corresponds to the equation (18), and the symbols correspond to the same quantities deduced from the particles.</figcaption></figure>

drasekar 1957). The total mass of the system is fixed equal to 1 in our numerical simulations. Tuning $n$ and $r_a$, we can then modify the only free physical parameters of the system: the size and the dynamical time. The size increases as $n$ for a given $r_a$, while the dynamical time increases as the anisotropy of the system for a given $n$. All the models are virial-relaxed ($\eta=2E_{\text{kinetic}}/E_{\text{potential}} = -1$). One family of DFs is plotted on Fig. 1. The DFs generated from such a procedure do not present any fixed point as in fig. 2(a) of Merritt (1985a). As a matter of fact, since we have not constrained the size of the system, the gravitational potential depends in our case on $r_a$ for a given $n$ (Fig. 2). Consequently, there is no fixed point.

In order to set up the initial conditions of our N-body numerical simulations, we now randomly choose the positions and the velocities for N particles from the DF. We plot in Fig. 3 the radial density profile and in Fig. 4 the velocity anisotropy deduced from one simulation of the system with $n=4$ and $r_a=2$. The agreement between our initial numerical conditions and the theoretical models are fully satisfactory.

The Ossipkov-Merrit models admit a fundamental limitation. As a matter of fact, for a given polytropic index $n$, there does exist a critical value of $r_a$ for which the DF becomes negative and unphysical in some region of the phase space (Fig. 1). Merrit (1985a) interprets this limitation as a simple illustration of the well-known fact that an arbitrary spherical mass distribution cannot always be reproduced by radial orbits. We have nevertheless considered models with $r_a=1$ or $r_a=0.75$ for example, which admit a negative DF in a region of phase space (Fig. 1). In these cases, however, we have arbitrarily set the DF equal to zero in this region, the positions and the velocities of the particles being thus chosen only after this truncation. As we shall see hereafter, our conclusions concerning the stability or the instability of such systems with a strong anisotropy have been obtained in the same way as for the realistic Ossipkov-Merrit systems with $r_a \ge 2$.

Finally, since each particle is initialized independently, the equilibrium DF $f_0(E,L^2)$ of the system is in fact slightly perturbed. The perturbation is caused by local Poissonian
fluctuations of the density. The dynamical evolution of the system represents then the response of an anisotropic self-gravitating spherical equilibrium system subject to such a perturbation.

## 3.2 Numerical integrations

Particles of mass $m_i$ and $m_j$ interact through the softened potential

$$ \phi_{ij} = Gm_i m_j / (r_{ij}^2 + \epsilon^2)^{1/2}, \quad (19) $$

where $r_{ij}$ is their separation. The softening parameter $\epsilon$ is essentially a particle radius. In our simulations, $\epsilon=0.05$. The most obvious algorithm to integrate equations (19) is direct summation (Aarseth 1972, 1985). Unfortunately, the number of operations in this algorithm grows as $N^2$. On classical sequential computers, it is at present not possible to consider a large number of particles ($N \le 10000$) and to perform many numerical simulations. On the other hand, the direct summation algorithm is particularly well adapted to massively parallel computers like a Connection Machine. In this paper, numerical simulations have been performed on a Connection Machine CM-5 by using a ‘Digital Orrery’ algorithm, first developed on a Connection Machine CM-2 (Hillis & Barnes 1987; Alimi & Scholl 1993; Serna, Alimi & Scholl 1994). In this case, one physical or virtual processor is assigned to each particle. One can then imagine two rings of processors both containing the same set of N-body coordinates. One ring starts to rotate stepwise according to the CM instructions. At each step, forces between all adjacent couples of bodies placed on the two different rings are calculated in parallel. A complete force calculation between all particles requires $N-1$ such steps. To illustrate the computational requirements, for a simulation with 8192 particles, the main part, computing all interacting forces between all particles, took less than 10 s per time-step on a 32-node CM-5. The accuracy of the integration can be illustrated by the variation of the total energy throughout a complete calculation, typically 30 dynamical times resolved by 10000 time-steps. This variation was smaller than 0.01 per cent. Hereafter, the set of numerical simulations performed have been made with 8192 particles. In order to get statistics to put error bars on our results, each model (a couple $n, r_a$) has had six independent analyses. Some experiments have been performed using more particles (65536); these show no significant changes in comparison with the work presented here.

# 4 NUMERICAL RESULTS

## 4.1 Morphological (in)stability

The physical mechanism of radial orbit instability for collisionless self-gravitating systems is well known. It is described in detail by several authors (Antonov 1962; Palmer 1994b). The morphological deformation resulting from such an instability is mainly caused by the trapping of particles with a low angular momentum in a bounded area of space. We can see such trapping in our two opposite stable ($n=4, r_a=100$) and unstable ($n=4, r_a=0.75$) cases in Fig. 5. To evaluate this deformation of an initial spherical system into an ellipsoid, it is convenient to use the axial ratio

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left diagram: the orbits of the two particles in a stable case $n=4, r_a=100$; we can see the regular precession owing to the angular momentum and the sphericity of the potential. Right diagram: by opposition, in an unstable $n=4, r_a=0.75$ case the orbit of a particle with low angular momentum can be trapped.</figcaption></figure>

defined from the moment of inertia tensor $I$ (Allen, Palmer & Papaloizou 1990). From the three real eigenvalues of $I$, $\lambda_1 \le \lambda_2 \le \lambda_3$, we compute the axial ratios $a_1 = \lambda_1/\lambda_2$ and $a_2 = \lambda_3/\lambda_2$. These two quantities, which can always be defined because an eigenvalue never vanishes, satisfy $a_1 \le 1 \le a_2$.

In Figs 6, 7 and 8 we plot the evolution of both the axial ratios and the virial ratio $\eta$ for three classes of systems defined by three different polytropic indices $n$. In each class, we have considered six different models which range from a purely radial model with $r_a=0.75$ to a quasi-isotropic one with $r_a=100$. A large range of systems with various physical properties (size, dynamical time...) (see Appendix), are thus taken into account in the set of these simulations.

All systems are initially spherical and virially relaxed. $a_1$ and $a_2$ are equal to 1, and $\eta$ is equal to -1 (see Section 3.1). The temporal average value for $\eta$ along the evolution is also equal to -1. During the first steps, however, at the same time as when the systems possibly deform, a significant fluctuation of $\eta$ appears clearly. It disappears later, $\eta$ staying equal to -1. Independently of the value of the polytropic index $n$, weakly anisotropic systems ($r_a > 2$) keep their initial spherical geometry during the evolution. Their axial ratios never differ significantly from unity ($a_1 \approx a_2 \approx 1$). They are morphologically stable. On the contrary – but still independently of $n$ – systems with a low initial anisotropic radius $r_a \le 2$ inevitably deform, reach a new non-spherical configuration, and stabilize in a new ellipsoidal configuration ($a_1$ or $a_2$ are significantly different from unity). This is exactly the effect of the radial orbit instability, which flattens the system (see de Zeeuw & Franx 1991 for a review). The collisionless hypothesis is fundamental for interpreting our results. Consequently, we have not continued our numerical simulations beyond a few hundred dynamical times in order to avoid a later evolution where two-body relaxation arises. However, all of our models reach a new relaxed state before $30 T_d$. We thus present our results for this interval. Moreover, it is important to notice on Figs 6, 7 and 8 the evolution of the virial ratio. This parameter, which stabilizes also at the -1 value before $30 T_d$, is a good indicator of the dynamical state of the system, which presents all the guarantees of no further dynamical evolution. Previous studies, which make no use of this parameter, must illustrate their results on much larger time-scales. In the next section, we interpret all of these numerical results with the help of the symplectic formalism described in Section 2. We also propose a stability parameter.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The axial (bold curves, left ordinates) and the virial (dashed curve, right ordinates) ratio versus dynamical time, for models with a polytropic index *n* = 3.5 and an anisotropy radius *r*<sub>a</sub> = 0.75 1, 2, 5, 10 and 100.</figcaption></figure>

## 4.2 Physical interpretation in terms of {$g_1, E$} and {$g_1, L^2$}

The analytical developments presented in Section 2 concern the dynamics and the properties in the mean field approximation of the six-variable DF $f$. In an *N*-body simulation we have access to the 6*N*-variable exact DF, $f^{(N)}(x_1, v_1, \dots, x_N, v_N)$ which is a solution of the Liouville equation (for a detailed justification of the approximations required for reconstructing $f$ from $f^{(N)}$, see Saslaw 1985)

$$ \frac{\partial f^{(N)}}{\partial t} = \sum_{\alpha=1}^{N} \left[ \nabla_x \cdot (f^{(N)} \dot{x}^\alpha) + \nabla_v \cdot (f^{(N)} \dot{v}^\alpha) \right] = 0. \quad (20) $$

Consequently, the Poisson brackets {$g_1, E$} and {$g_1, L^2$} appearing in equation (14), defining the linear perturbation, are estimated in *N*-body simulations from statistical properties of the random variables $\epsilon_i$ and $\lambda_i$ defined for each particle *i*.

$$ \epsilon_i := \{g_1, E\}_{x=x_i}^{v=v_i} = (\nabla_x g_1 \cdot \mathbf{v})_i - (\nabla_v g_1 \cdot \nabla_x \mathbf{v})_i = \left(\frac{\partial E}{\partial \mathbf{v}} \cdot \mathrm{d}\mathbf{v}\right)_i \\ + \left(\frac{\partial E}{\partial \mathbf{x}} \cdot \mathrm{d}\mathbf{x}\right)_i = dE_i, \quad (21) $$

and

$$ \lambda_i := \{g_1, L^2\}_{x=x_i}^{v=v_i} = (\nabla_x g_1 \cdot \nabla_v L^2)_i - (\nabla_v g_1 \cdot \nabla_x L^2)_i = \left(\frac{\partial L^2}{\partial \mathbf{v}} \cdot \mathrm{d}\mathbf{v}\right)_i \\ + \left(\frac{\partial L^2}{\partial \mathbf{x}} \cdot \mathrm{d}\mathbf{x}\right)_i = dL_i^2. \quad (22) $$

$\epsilon_i$ and $\lambda_i$ are computed from positions and velocities of each particle *i* at the initial time $t_0$ and at a time $t_1 = t_0 + \delta t$. In order to compare different models (defined by different *n* and $r_a$) we tune $\delta t$ such that $\delta f/f = 0.01$. In all cases, we have $\delta t < T_d/100$.

The variable $\epsilon_i$ is related to the energy perturbation seen by the system (equation 21) and thus it is related to the stability of the system. In Fig. 9, we plot, for different models, the fraction of $\epsilon_i$ in the numerical simulations at initial time with a negative value, i.e., the probability $P_{nb}(\epsilon)$ that $\epsilon$ is negative.

The variable $\lambda_i$ is related to the anisotropic component of the linear variation of the DF of the system (equation 22). We have shown (Section 2, Perez 1995; Paper I) that the systems have particular stability properties against all preserving perturbations, i.e., the perturbations generated by $g_1$ such that {$g_1, L^2$} = 0. Except for the purely isotropic stable case, however, there exist no physical systems that are sub-
<figure><img src="image_8.png" /></figure>

**Figure 7.** The axial (bold curves, left ordinates) and the virial (dashed curve, right ordinates) ratio versus dynamical time, for models with a polytropic index *n* = 4 and an anisotropy radius *r*<sub>a</sub> = 0.75 1, 2, 5, 10 and 100.

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

**Figure 8.** The axial (bold curves, left ordinates) and the virial (dashed curve, right ordinates) ratio versus dynamical time, for models with a polytropic index *n* = 4.5 and an anisotropy radius *r*<sub>a</sub> = 0.75 1, 2, 5, 10 and 100.

© 1996 RAS, MNRAS 280, 700–710
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Probability $P_{\text{nb}}(\epsilon)$ for $\epsilon$ to be negative (dashed curves, left ordinate) and Pearson index of the random variable $\lambda$ (plain curves, right ordinates) versus anisotropy radius, for models with a polytropic index $n=3.5, 4$ and $4.5$, respectively.</figcaption></figure>

mitted only to preserving perturbations. In the opposite case, a fully anisotropic system with all components evolving only on radial orbits, is unstable. Consequently, in order to characterize the preserving nature of the perturbations acting on the system, we suggest the use of the distribution of $\lambda_i$ around their vanishing mean value. In order to discriminate between the preserving and non-preserving nature of perturbations we must take into account not only how weakly the $\lambda_i$ are scattered around their mean, but also how highly they are peaked at their mean. We thus compute the statistical Pearson index of the random variable $\lambda$ (hereafter $P_{\lambda}$). This parameter is both a shape and a scattering parameter. A flattened distribution (platikurtic) of a random variable is characterized by a small or negative Pearson index. A distribution which is both concentrated around its mean value and which is highly peaked (leptokurtic) has a large Pearson index (Calot 1973). However, we are interested only in the distribution of $\lambda_i$ around their mean value and, in order to eliminate some aberrant contributions, we define a truncated variable

$$ \bar{\lambda}_i = \begin{cases} \lambda_i & \text{when } |\lambda_i| \le k\sigma_\lambda \\ \bar{\lambda} & \text{when } |\lambda_i| > k\sigma_\lambda \end{cases} $$

where $k$ is such that $K$ per cent of particles are taken into account, with

$$ K \text{ per cent} = \min \left( 95 \text{ per cent}, \frac{\bar{v}_T}{3\bar{v}_R} \times 100 \right), $$

$\sigma_\lambda$ is the standard deviation of the random variable $\lambda_i$ and $\bar{\lambda}$ is the mean of $\bar{\lambda}_i$. We compute the statistical Pearson index

of the variable $\bar{\lambda}_i$.

$$ P_{\bar{\lambda}} = \frac{\sum_{i=1}^{N} (\bar{\lambda}_i - \bar{\lambda})^4}{\left[ \sum_{i=1}^{N} (\bar{\lambda}_i - \bar{\lambda})^2 \right]^2} - 3, \quad (23) $$

with $\bar{\lambda}$ being the mean $\bar{\lambda}_i$.

In Fig. 9 we plot $P_{\bar{\lambda}}$ for different models. The error bars both for $P_{\bar{\lambda}}$ and $P_{\text{nb}}(\epsilon)$ are calculated from six numerical simulations of a given model. The bars represent $\pm \sigma$ variations (i.e., 95 per cent of the values for these non-Gaussian variables).

For any polytropic index $n$, the respective behaviours of $P_{\text{nb}}(\epsilon)$ and $P_{\bar{\lambda}}$ are similar. As a matter of fact, the strongly anisotropic initial systems with $r_a \lesssim 2$ which are morphologically unstable (Section 4.1, Figs 5, 6 and 7), are characterized by a large $P_{\text{nb}}(\epsilon)$,

$$ P_{\text{nb}}(\epsilon) \gtrsim 20 \text{ per cent}, \quad (24) $$

and by a small $P_{\bar{\lambda}}$

$$ P_{\bar{\lambda}} \lesssim 2.5. \quad (25) $$

On the contrary, the median and weakly anisotropic initial systems ($r_a \gtrsim 2$), which are morphologically stable (Section 4.1, Figs 5, 6 and 7), are characterized by a small $P_{\text{nb}}(\epsilon)$,

$$ P_{\text{nb}}(\epsilon) \lesssim 20 \text{ per cent}, \quad (26) $$

and a large $P_{\bar{\lambda}}$

$$ P_{\bar{\lambda}} \gtrsim 2.5. \quad (27) $$
<table><caption>**Table 1.** Physical characteristic of used general anisotropic spherical systems.</caption><thead><tr><th></th><th>R<sub>1/2</sub></th><th>T<sub>d</sub></th><th>P<sub>nb</sub>(ε)</th><th>P<sub>λ&#773;</sub></th></tr></thead><tbody><tr><td>Model 1</td><td>2.69</td><td>4,42</td><td>18,75</td><td>8,26</td></tr><tr><td>Model 2</td><td>9.94</td><td>31,35</td><td>22,82</td><td>2,12</td></tr><tr><td>Model 3</td><td>3.27</td><td>5,94</td><td>17,46</td><td>2,93</td></tr><tr><td>Model 4</td><td>3.18</td><td>5,68</td><td>18,04</td><td>3,42</td></tr></tbody></table>

Moreover, these two classes of systems are clearly separated with respect to the logarithmic scale used for $r_a$.

## 5 UNIVERSALITY OF OUR PREDICTIONS AND CONCLUSIONS

In order to generalize to all collisionless self-gravitating spherical systems the stability criterion – defined in terms of the parameters $P_{nb}(\epsilon)$ and $P_{\lambda\lambda'}$ – which has been obtained for

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Lindblad diagram for polytropic Ossipkov–Merritt models with $n=4.5$, $r_a=1, 2, 10$ and $100$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Lindblad diagram for models 1, 2, 3 and 4.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Axial ratio (top diagram) and virial ratio (bottom diagram) versus dynamical time for models 1, 2, 3 and 4.</figcaption></figure>

the Ossipkov–Merritt models (equations 24–27), we have generated four independent spherical models. These models do not follow an Ossipkov–Merritt distribution (Section 3.1). They are virialized and present an anisotropy in the velocity space. We show in Table 1 the physical characteristics (size and dynamical time) of these different configurations. Since all these models are spherical and anisotropic, their DF depends both on *E* and *L*². We plot the isocontour of this function in the *E*–*L*² plane, thus forming the so-called Lindblad diagram of the system (see Merritt 1985a and Lindblad 1933). For comparison, we have plotted in Fig. 10 a set of such diagrams for the Ossipkov–Merritt models with *n* = 4.5, *r*<sub>a</sub> = 1, 2, 10 and 100. We show in Fig. 11 the Lindblad diagram of the systems described in Table 1.

We have considered four models in order to illustrate the three classes of foreseeable spherical systems: quasi-isotropic models (M1 and M3); circular anisotropic models (M4); and a non-Ossipkov–Merritt radial anisotropic system (M2).

Models 1 and 3 are weakly anisotropic. They are similar to the Ossipkov–Merritt systems (respectively with *n* = 4.5, *r*<sub>a</sub> = 10 and *n* = 4.5, *r*<sub>a</sub> = 100). Model 2 represents a system with particles mainly on radial orbits. Its DF is highly peaked around *L*²=0. An analysis in velocity space of model 4 shows that this model presents a strong circular orbit anisotropy.

In Fig. 12, we present the evolution of the axial and virial ratios for models 1 to 4. Only model 2 is morphologically

unstable, contrary to model 1 which is very stable. We have computed for these models the parameters *P*<sub>nh</sub>(*ε*) and *P*<sub>z</sub>. The stability criteria (equations 24–27) are confirmed. For the unstable model 2, parameter *P*<sub>nh</sub>(*ε*) is larger than 20 per cent (Table 1) and parameter *P*<sub>z</sub> is smaller than 2.5 (Table 1). The very stable model 1 is characterized by a low *P*<sub>nh</sub>(*ε*) (=18.75 per cent < 20 per cent) and by a very large *P*<sub>z</sub> (=8.26 > 2.5) (Table 1). We are thus able to generalize our stability criteria from the previous section to all collisionless self-gravitating spherical systems.

(i) All anisotropic collisionless self-gravitating spherical systems with parameter *P*<sub>nh</sub>(*ε*) smaller than 20 per cent are stable. As we do not have analytical suggestions for the opposite case, we do not propose any conclusions on the instability of such systems from parameter *P*<sub>nh</sub>(*ε*).

(ii) The dynamical evolution of all anisotropic collisionless self-gravitating spherical systems with $P_{\bar{z}} < 2.5$ is dominated by non-preserving perturbations (Section 2). Such systems are unstable.

## REFERENCES

Aarseth S. J., 1972, in Lecar M., ed., IAU Symp. 10, Gravitational N-Body Problem. Reidel, Dordrecht

Aarseth S. J., 1985, in Brackbill J. U., Cohen B. H., eds, Multiple Time Scale. Academic Press, New York

Alimi J. M., Scholl H., 1993, Int. J. Mod. Phys. C, 4, 197

Allen A. J., Palmer P. L., Papaloizou J., 1990, MNRAS, 242, 576

Aly J. J., Perez J., 1992, MNRAS, 259, 95

Antonov V.A., 1962, Vest. Leningrad Univ., 7, 135

Antonov V A., 1973, in Omarov G. B., ed., Nauka, Alma Ata, The Dynamics of Galaxies and Star Clusters, p. 139

Arnold V. I., 1978, Mathematical Methods of Classical Mechanics. Springer-Verlag, New York

Barnes J., 1986, in Goodman J., Hut P., eds, IAU Symp. 113, Dynamics of Stellar Systems. Reidel, Dordrecht, p. 297

Barnes J., Goodman J., Hut P., 1986, ApJ, 300, 112

Bartholomew P., 1971, MNRAS, 151, 333

Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton Univ. Press, Princeton, NJ

Calot G., 1973, Cours de Statistique Descriptives. Dunod, Paris

Chandrasekhar S., 1957, An Introduction to the Study of the Stellar Structure. Dover, New York

Dejonghe H., Merritt D., 1993, ApJ, 328, 93

de Zeeuw P. T., Franx M., 1991, ARA&A, 29, 239

Doremus J. P., Feix M. R., Baumann G., 1971, Phys. Rev. Lett., 26, 725

Fridman A. M., Polyachenko V. L., 1984, Physics of Gravitating Systems. Springer-Verlag, New York

Goodman J., 1988, ApJ, 329, 612

Hénon M., 1973, A&A, 24, 229

Hillis W. D., Barnes J., 1987, Nat, 326, 27

Kandrup H. E., 1990, ApJ, 351, 104

Kandrup H. E., 1991a, ApJ, 370, 312

Kandrup H. E., 1991b, ApJ, 380, 511

Kandrup H. E., Sygnet J. F., 1985, ApJ, 298, 27

Kulsrud R. M., Mark J. W. K., 1970, ApJ, 160, 471

Larsson J., 1991, Phys. Rev. Lett., 66, 1466

Laval G., Mercier C., Pellat R., 1965, Nuclear Fusion, 5, 156

Lindblad B., 1933, Handbuch des Astrophysik, Vol. 5. Springer-Verlag, Berlin

Lynden-Bell D., Sannit N., 1968, MNRAS, 143, 167

Merritt D., 1985a, AJ, 86, 1027

Merritt D., 1985b, MNRAS, 214, 25P
Merritt D., Aguilar L. A., 1985, MNRAS, 217, 787

Morrisson P. J., 1980, Phys. Lett. A, 80, 383

Ossipkov L. P., 1979, Pis'ma Astron. Zh., 5, 77

Palmer P. L., 1994a, in Contopoulos, ed, Proc. EADN Astrophys. School VI, Galactic Dynamics and N-Body Simulations. Lecture Notes Phys. 433, Springer-Verlag, New York

Palmer P. L., 1994b, Stability of Collisionless Stellar Systems. Kluwer, Dordrecht

Palmer P. L., Papaloizou J., 1986, MNRAS, 224, 1043

Palmer P. L., Papaloizou J., 1988, MNRAS, 231, 935

Perez J., 1995, PhD Thesis, University Paris VII, Denis Diderot

Perez J., Aly J. J., 1996, MNRAS, 280, 689 (Paper I, this issue)

Polyachenko V. L., 1992, Sov. Phys. JETP, 74, 755

Saslaw W. C., 1985, Gravitational Physics of Stellar and Galactic Systems. Cambridge Univ. Press, Cambridge

Serna A., Alimi J.-M., Scholl H., 1994, ApJ, 427, 574

Sygnet J. F., Des Forets G., Lachieze-Rey M., Pellat R., 1984, ApJ, 276, 737

## APPENDIX A: STATISTICS OF SIMULATIONS

We have randomly chosen the position and the velocity for 8192 particles from the DF. Six independent simulations have been carried out for each value of *n* and *r*<sub>a</sub>. We present in the following table (Table A1) the physical properties and the stability parameters calculated for a large class of generalized polytropic Ossipkov-Merritt models. The dynamical time is the ratio between the typical size of the system and the modulus of median velocity of a particle. *R*<sub>1/2</sub> represents the radius of the sphere containing half of the system total mass. *P*<sub>nb</sub>(*ε*) represents the probability for the {g<sub>1</sub>, E} calculated for each particle to be negative (expressed in per cent). Finally, *P*(*λ*) represents the Pearson index of the distribution of {g<sub>1</sub>, L<sup>2</sup>}.

The number indicates the mean value of the quantity over the six simulations, and the error represents an interval of two standard deviations of the quantity. Units are such that the gravitational constant, the total mass of the system and the initial value of the relative gravitational potential are unity.

<table><caption>**Table A1.** Physical characteristics for models with *n* = 3.5, 4, 4.5 and *r*<sub>a</sub> = 0.75, 1, 2, 5, 10 and 100.</caption><thead><tr><th>n = 3.5</th><th>T<sub>d</sub></th><th>R<sub>1/2</sub></th><th>P<sub>nb</sub>(&epsilon;)</th><th>P(&lambda;)</th></tr></thead><tbody><tr><td>r<sub>a</sub> = 0.75</td><td>5.455 &plusmn; 0.062</td><td>3.099 &plusmn; 0.023</td><td>23.04 &plusmn; 0.72</td><td>0.953 &plusmn; 0.053</td></tr><tr><td>r<sub>a</sub> = 1</td><td>5.127 &plusmn; 0.057</td><td>2.973 &plusmn; 0.022</td><td>21.03 &plusmn; 0.73</td><td>1.205 &plusmn; 0.064</td></tr><tr><td>r<sub>a</sub> = 2</td><td>4.157 &plusmn; 0.039</td><td>2.585 &plusmn; 0.016</td><td>16.38 &plusmn; 0.61</td><td>2.088 &plusmn; 0.076</td></tr><tr><td>r<sub>a</sub> = 5</td><td>3.238 &plusmn; 0.022</td><td>2.189 &plusmn; 0.010</td><td>12.87 &plusmn; 0.55</td><td>2.85 &plusmn; 0.13</td></tr><tr><td>r<sub>a</sub> = 10</td><td>3.043 &plusmn; 0.020</td><td>2.0997 &plusmn; 0.0095</td><td>12.38 &plusmn; 0.46</td><td>3.06 &plusmn; 0.13</td></tr><tr><td>r<sub>a</sub> = 100</td><td>2.973 &plusmn; 0.020</td><td>2.0676 &plusmn; 0.0097</td><td>12.10 &plusmn; 0.48</td><td>3.082 &plusmn; 0.13</td></tr></tbody><thead><tr><th>n = 4</th><th>T<sub>d</sub></th><th>R<sub>1/2</sub></th><th>P<sub>nb</sub>(&epsilon;)</th><th>P(&lambda;)</th></tr></thead><tbody><tr><td>r<sub>a</sub> = 0.75</td><td>8.720 &plusmn; 0.096</td><td>4.237 &plusmn; 0.031</td><td>25.63 &plusmn; 0.60</td><td>0.935 &plusmn; 0.082</td></tr><tr><td>r<sub>a</sub> = 1</td><td>8.257 &plusmn; 0.087</td><td>4.085 &plusmn; 0.029</td><td>23.45 &plusmn; 0.72</td><td>1.079 &plusmn; 0.044</td></tr><tr><td>r<sub>a</sub> = 2</td><td>6.473 &plusmn; 0.037</td><td>3.472 &plusmn; 0.013</td><td>18.74 &plusmn; 0.57</td><td>1.962 &plusmn; 0.065</td></tr><tr><td>r<sub>a</sub> = 5</td><td>4.503 &plusmn; 0.037</td><td>2.727 &plusmn; 0.015</td><td>14.77 &plusmn; 0.64</td><td>3.04 &plusmn; 0.15</td></tr><tr><td>r<sub>a</sub> = 10</td><td>4.015 &plusmn; 0.033</td><td>2.526 &plusmn; 0.014</td><td>13.98 &plusmn; 0.60</td><td>3.40 &plusmn; 0.12</td></tr><tr><td>r<sub>a</sub> = 100</td><td>3.835 &plusmn; 0.031</td><td>2.450 &plusmn; 0.013</td><td>13.59 &plusmn; 0.55</td><td>3.53 &plusmn; 0.16</td></tr></tbody><thead><tr><th>n = 4.5</th><th>T<sub>d</sub></th><th>R<sub>1/2</sub></th><th>P<sub>nb</sub>(&epsilon;)</th><th>P(&lambda;)</th></tr></thead><tbody><tr><td>r<sub>a</sub> = 0.75</td><td>17.09 &plusmn; 0.21</td><td>6.634 &plusmn; 0.053</td><td>26.67 &plusmn; 0.38</td><td>1.161 &plusmn; 0.070</td></tr><tr><td>r<sub>a</sub> = 1</td><td>16.26 &plusmn; 0.19</td><td>6.418 &plusmn; 0.050</td><td>24.79 &plusmn; 0.56</td><td>1.241 &plusmn; 0.042</td></tr><tr><td>r<sub>a</sub> = 2</td><td>12.56 &plusmn; 0.12</td><td>5.403 &plusmn; 0.035</td><td>19.78 &plusmn; 0.60</td><td>1.898 &plusmn; 0.094</td></tr><tr><td>r<sub>a</sub> = 5</td><td>7.206 &plusmn; 0.081</td><td>3.730 &plusmn; 0.028</td><td>16.12 &plusmn; 0.61</td><td>3.16 &plusmn; 0.13</td></tr><tr><td>r<sub>a</sub> = 10</td><td>5.655 &plusmn; 0.062</td><td>3.174 &plusmn; 0.023</td><td>15.14 &plusmn; 0.57</td><td>3.85 &plusmn; 0.13</td></tr><tr><td>r<sub>a</sub> = 100</td><td>5.064 &plusmn; 0.051</td><td>2.949 &plusmn; 0.020</td><td>14.73 &plusmn; 0.57</td><td>4.18 &plusmn; 0.14</td></tr></tbody></table>