# Stability of rotating spherical stellar systems

J.-M. Alimi, Jérôme Perez, A. Serna

► To cite this version:

J.-M. Alimi, Jérôme Perez, A. Serna. Stability of rotating spherical stellar systems. Monthly Notices of the Royal Astronomical Society, 1999, 305 (4), pp.859-865. 10.1046/j.1365-8711.1999.02475.x. hal-01010749

HAL Id: hal-01010749

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

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 rotating spherical stellar systems

Jean-Michel Alimi,¹ Jérôme Perez²* and Arturo Serna³

¹Laboratoire d'Astrophysique Extragalactique et de Cosmologie, CNRS UMR 8631, Observatoire de Meudon, 5 Place Jules Jansen, 92195 Meudon, France

²Ecole Nationale Supérieure des Techniques Avancées, SMP, CNRS URA D853, 32 Boulevard Victor, 75996 Paris, France

³Universidad Miguel Hernández, Area de Física Aplicada, Edif. La Galia, 03291-Elche, Spain

Accepted 1999 January 19. Received 1998 November 9; in original form 1998 January 23

## ABSTRACT

We study the stability of rotating collisionless self-gravitating spherical systems by using high-resolution N-body experiments on a Connection Machine CM-5.

We added rotation to Ossipkov–Merritt (OM) anisotropic spherical systems by using two methods. The first method conserves the anisotropy of the distribution function defined in the OM algorithm. The second method distorts the systems in velocity-space. We then show that the stability of systems depends both on their anisotropy and on the value of the ratio of the total kinetic energy to the rotational kinetic energy. We also test the relevance of the stability parameters introduced by Perez et al. for the case of rotating systems.

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

## 1 INTRODUCTION

Various analytical and numerical studies (Antonov 1973; Barnes, Goodman & Hut 1986; Palmer & Papaloizou 1986; Perez & Aly 1996; Perez et al. 1996 and references therein) have shown that spherical, collisionless, self-gravitating anisotropic systems with components moving mainly on radial orbits are unstable. However, all these works considered non-rotating systems, and it is well known that rotation can play an important role in the dynamical evolution of systems and modify their stability properties. It has been shown that rotation can be the cause of the deformation of systems like globular clusters or weakly elliptical galaxies (Staneva, Spassova & Golev 1996; Goodwin 1997).

The stability of rotating stellar systems is a very complex problem. Much work has been concerned with barred galaxies (which are rapidly rotating stellar systems) (for a review see, Buta & Crocker & Elmegreen 1996), but in a more general context, few studies have been devoted in the literature to this topic. For example, Papaloizou, Palmer & Allen (1991) have performed a series of numerical simulations to analyse the stability of systems where rotation was introduced by using the technique proposed by Lynden-Bell (1962). All their simulations produced end-states in which a triaxial bar appears. These important results cannot be considered general, since they were obtained for systems dominated by particles evolving on radial orbits, and put in rotation by a specific procedure. In order to analyse the influence of rotation on the (in)stability of a given system, it is necessary to consider not only spherical systems with different kinds of anisotropy but also different methods for introducing the rotation. This paper develops such an analysis. We are also interested in testing the relevance of the stability parameters (Perez et al. 1996) on rotating systems.

Perez et al. have shown that the stability of spherical self-gravitating non-rotating systems can be deduced from the ‘anisotropic’ component of the linear variation of the distribution function (see Section 2.2). Such stability parameters can be computed from rotating systems. We show that they are still relevant for anisotropic systems as long as the rotational kinetic energy is not too large.

The paper is arranged as follows. We describe in Section 2 the method that we use to obtain the initial non-rotating systems as well as the parameters describing the (in)stability of such systems. In Section 3 we detail the techniques used to introduce a parametrizable rotation to the initial conditions presented in Section 2. In Section 4 we show our numerical results on the (in)stability of various rotating systems generated with different procedures. Finally, the discussion and physical interpretation of our results are presented in Section 5.

## 2 STABILITY AND INSTABILITY OF NON-ROTATING SYSTEMS

### 2.1 Non-rotating initial conditions

In a previous paper (Perez et al. 1996) we used the Ossipkov–Merritt (OM) algorithm (Ossipkov 1979; Merritt 1985a,b; Binney & Tremaine 1987) to generate anisotropic self-gravitating spherical systems with various physical properties. This algorithm starts from a density given by $\rho_{iso}(r) \propto \psi_{iso}^n$, where $\psi_{iso}(r)$ is a known gravitational potential satisfying the Poisson equation, while $n$ is the polytropic index ($1/2 < n \le 5$). This density profile $\rho_{iso}(r)$ is then deformed according to

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

where the anisotropic radius $r_a$ is a parameter which controls the deformation.

<footnote>*E-mail: perez@enstay.ensta.fr</footnote>
Using the Abel inversion technique, this procedure allows one to define an anisotropic distribution function (DF) which depends both on $E$ and $L^2$ through the variable $Q := E + L^2/2r_a^2$,

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

Once this DF has been computed, the initial conditions of our N-body numerical simulations are generated by choosing at random, from the above DF, the positions and velocities for the N particles. The density profile $\rho_{ani}(r)$ defined by equation (1) is the probability density from which the positions are generated. The velocities are generated from the velocity probability density deduced from the equation (2) (see Appendix A).

It must be noted that there is a fundamental limitation in the OM models: any given value of the polytropic index $n$ implies a critical value of $r_a$ below which the DF becomes negative and unphysical in some region of phase-space. Merritt (1985a) interprets this limitation as a simple illustration of the well-known fact that radial orbits cannot always reproduce an arbitrary spherical mass distribution. In these cases, in order to extend the OM algorithm to highly radially anisotropic ($r_a \approx 1$) systems, we have arbitrarily set the DF equal to zero in this region. Such a procedure on DF affects only particles with a large value of $Q$. This procedure is applied for systems with a small value of $r_a$ which contain mainly particles with a small value of $Q$. Such a procedure then affects a very small number of particles (less than 0.1 per cent of the total number of particles). The systems with a modified DF are not strictly OM systems. However, they conserve the properties which are for the present work: the density profiles deduced from the modified DF are indistinguishable from the density profiles given by equation (1) with the same value of $r_a$, the Lindblad diagrams are very peaked around a small value of $Q$ (Perez et al. 1996), they well correspond to highly radially anisotropic systems, and finally the radial dependence of the velocity anisotropy

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

is preserved.

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

## 2.2 Stability analysis

The equilibrium DF of a spherical self-gravitating system depends only on the one-particle energy $E$ and the squared total angular momentum $L^2$. If $g_1$ denotes the perturbation generator, the linear variation of the DF can be written as

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

If DF is a monotonic decreasing function,<fnref n="1" /> the stability is then related to the Poisson brackets $\{g_1, E\}$ and $\{g_1, L^2\}$ (Perez & Aly 1996; Perez et al. 1996). In our N-body simulations, these quantities appear as two random variables, $\epsilon$ and $\lambda$ respectively, defined for

each particle $i$ (Perez et al. 1996). The stability of the system can be predicted from the probability $P_\epsilon$ for $\epsilon$ to be negative, and the statistical Pearson index (Calot 1973) $P_\lambda$ of the variable $\lambda$. All anisotropic, collisionless, self-gravitating, non-rotating spherical systems with

$$P_\lambda \le 2.5 \text{ and } P_\epsilon \ge 20 \text{ per cent} \quad (5)$$

are unstable, while those with

$$P_\lambda \ge 2.5 \text{ and } P_\epsilon \le 20 \text{ per cent} \quad (6)$$

are stable. The two other regions of the $(P_\lambda P_\epsilon)$ plane correspond to a transition between a stable system and an unstable system. In the particular case of OM models, these regions correspond to an anisotropy radius $r_a$ close to 2 (Perez et al. 1996).

# 3 GENERATION OF ROTATING SYSTEMS

## 3.1 Definition

In order to generate virial-relaxed rotating spherical systems, we modify the non-rotating systems defined in the previous section by using techniques derived from the Lynden-Bell method (1962). Since this method preserves the position and the norm of the velocity of each particle, the systems are put in rotation without modifying their total potential and kinetic energy. In practice, we apply the following transformations to the velocity components $\{v_r, v_\theta, v_\phi\}$ of the particles:

Method 1 Method 2

$$\begin{align}
v_r &\rightarrow v_r & v_r &\rightarrow 0 \nonumber \\
v_\theta &\rightarrow v_\theta & v_\theta &\rightarrow v_\theta \nonumber \\
v_\phi &\rightarrow |v_\phi| & v_\phi &\rightarrow \sqrt{1 + \frac{v_r^2}{v_\phi^2}} |v_\phi|. \tag{7}
\end{align}$$

The first method then conserves the radial anisotropy defined in the OM algorithm, while the second method distorts the system in velocity-space.

The amount of rotation introduced by these methods can be evaluated through the ratio

$$\mu = K_{\text{rot}}/K_{\text{tot}}, \quad (8)$$

where $K_{\text{tot}}$ is the total kinetic energy, and $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{(L_i \cdot \hat{L}_{\text{tot}})^2}{[r_i^2 - (r_i \cdot \hat{L}_{\text{tot}})^2]}. \quad (9)$$

Here, $L_i$ is the specific angular momentum of particle $i$, and $\hat{L}_{\text{tot}}$ is a unit vector in the direction of the total angular momentum of the system; when the system does not rotate at all, the $\hat{L}_{\text{tot}}$ vector is the null vector, and $N$ is the total number of particles. In order to exclude counter-rotating particles, the sum in equation (9) is actually carried out only over those particles satisfying the condition $(L_i \cdot \hat{L}_{\text{tot}}) > 0$.

In order to have systems with different strengths of homogeneous rotation (HR), we have applied either Method 1 or Method 2 to a fraction $\tau$ of the total number of particles. This fraction has been constructed by choosing the particles at random in the overall system. When $\tau \to 0$, the system does not turn, while, when $\tau \to 1$, the system rotation reaches its maximum value.

In principle, there is no reason to consider only the case of homogeneous rotation. Moreover, in order to roughly model the presence of a rotating massive object like those sometimes considered in the centre of some elliptical galaxies, we also consider inhomogeneous rotation (IR). In this case, we apply the above

<footnote n="1">¹We consider throughout this paper only systems with a DF which admits a monotonous decreasing dependence with respect to all the isolating integrals of motion ($\partial f_o/\partial E < 0$, and $\partial f_o/\partial L^2 < 0$).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** $\mu$ versus $\tau$ for $n = 4$; $HR_1$, $r_a = 1.5$ (thin solid line); $HR_1$, $r_a = 100$ (thin dashed line) and $HR_2$, $r_a = 1.5$ (bold solid line); $HR_2$, $r_a = 100$ (bold dashed line).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** $\mu$ versus $k$ for $n = 4$, $HR_1$, $r_a = 1.5$ (thin solid line); $HR_1$, $r_a = 100$ (thin dashed line) and $HR_2$, $r_a = 1.5$ (bold solid line); $HR_2$, $r_a = 100$ (bold dashed line).</figcaption></figure>

velocity transformations only to those particles placed at a radial distance smaller than $k \times R_1$, where $k$ is a positive parameter and $R_1$ is the radius containing half of the system mass. If $k = 0$, the system does not turn, while, if $k \to +\infty$, the rotation has its maximum value. We have then four possible procedures for introducing a rotation motion on our initial conditions. The first two possibilities introduce a homogeneous rotation by choosing particles at random in the whole system and modifying their velocities according either to the method 1 or 2, which defines the $HR_1$ and $HR_2$ procedures respectively. The other two possibilities introduce an inhomogeneous rotation by applying either method 1 or method 2 to modify the velocities of those particles placed within a given radial distance, which defines the procedures $HR_1$ and $HR_2$ respectively.

Figs 1 and 2 show the values of $\mu$ (equation 8) obtained from the four possible procedures. As we can see from such figures, only the $HR_2$ and $HR_2$ procedures lead to large fractions ($\mu \ge 10$ per cent) of kinetic rotation energy. We also note from these figures that the dependence of $\mu$ on $r_a$ depends on whether velocities have been modified according to method 1 or method 2. As a matter of fact, for the $HR_1$ and $HR_1$ procedures, the amount of rotation obtained ($\mu$) is greater for large $r_a$ values than for small ones. On the contrary, for $HR_2$ and $HR_2$, $\mu$ is larger for small $r_a$ values than for large ones.

## 4 INFLUENCE OF THE ROTATION ON THE STABILITY

Using the N-body code described in Alimi & Scholl (1993), we have performed on Connection-Machine 5 a series of numerical simulations<fnref n="2" /> of the evolution of the systems defined in the previous section. As the collisionless hypothesis is fundamental for interpreting our results, we have not continued our simulations beyond a few hundred dynamical times in order to avoid the later evolution where two-body relaxation arises. However, all our models reach a steady state before about 50 $T_d$ (where the initial dynamical time is estimated by the following formula $T_d = \sqrt{\sum v_i^2 / \sum v_i^2}$, and the summations on initial positions and velocities are done on all the particles). We will then present our results for this interval.

The physical mechanism of the radial-orbit instability for collisionless self-gravitating systems is well known. It has been described in detail by several authors (see, e.g., Palmer 1994). The morphological deformation of the initial gravitational system resulting from this instability is mainly due to the trapping of particles with a low angular momentum in a bounded area of space. This trapping favours a deformation of the initial spherical system to an ellipsoidal or even a bar-like structure. To evaluate such a deformation, it is convenient to use the axial ratio defined from the moment of inertia tensor $I$ (Allen, Palmer & Papaloizou 1990). From the three real eigenvalues of $I$, $\lambda_1 \ge \lambda_2 \ge \lambda_3$, we compute the axial ratios $a_1 = \lambda_2/\lambda_1$ and $a_2 = \lambda_2/\lambda_3$. These two quantities, which can always be defined because these eigenvalues never vanish, satisfy $a_1 \le 1 \le a_2$. In order to discriminate clearly between a bar-like structure, a quasi-sphere and a disc-like structure, we define the quantity $f$ from $a_1$ and $a_2$:

$$f = \frac{1 - a_1}{a_2 - 1} \qquad (10)$$

A bar-like structure is characterized by $a_1 < 1$ and $a_2 \approx 1$, which implies an $f$-value significantly larger than 1. A disc-like structure is characterized by $a_1 \approx 1$, $a_2 > 1$ and $0 < f < 1$. Any system with a $f$-value of order unity has a quasi-spherical structure.

### 4.1 Rotating systems according to method 1

This type of rotation preserves the anisotropy of the non-rotating OM systems. The distribution function of the rotating system depends only on the variable $Q = E + L^2/2r_a^2$ as in the case of the non-rotating systems. We see in Table 1 that, in the case $n = 4$, the stability parameters defined in Section 2.2 are very weakly modified whatever the $\tau$ and $k$ parameters values are, that is to say, whatever the rotational kinetic energy is (low with this method). According to the conditions given by equations (5) and (6), we expect the (in)stability of the original non-rotating systems not to be modified when they are put in rotation. Our numerical simulations confirm this. In Fig. 3 we see that the evolution of axial ratios is similar for the rotating (dashed and dotted lines) and non-rotating (solid lines) systems. This results holds for the homogeneous and inhomogeneous rotations whatever the $n$ parameter value is.

### 4.2 Rotating systems according to method 2

The situation is now more complicated because the rotation procedure modifies the system's anisotropy. In the second

<footnote n="2">²The set of numerical simulations performed have been made with 16384 particles. Some experiments have been performed using more particles (65536); no significant change in comparison with the work presented here have been obtained.</footnote>
<table><caption>**Table 1.** Stability parameter evolution for the first method for rotating systems with $n = 4$.</caption><thead><tr><th colspan="5">HR<sub>1</sub></th><th colspan="5">IR<sub>1</sub></th></tr><tr><th rowspan="2">τ</th><th colspan="3">r<sub>a</sub> = 1.5</th><th colspan="3">r<sub>a</sub> = 100</th><th colspan="3">r<sub>a</sub> = 1.5</th><th colspan="3">r<sub>a</sub> = 100</th></tr><tr><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>0.0</td><td>22.42</td><td>1.39</td><td>13.84</td><td>4.07</td><td>0.4</td><td>24.54</td><td>1.31</td><td>15.25</td><td>3.89</td><td></td><td></td><td></td></tr><tr><td>0.2</td><td>22.83</td><td>1.41</td><td>14.28</td><td>4.12</td><td>0.6</td><td>24.04</td><td>1.31</td><td>14.60</td><td>3.84</td><td></td><td></td><td></td></tr><tr><td>0.4</td><td>22.76</td><td>1.39</td><td>14.05</td><td>4.23</td><td>0.8</td><td>24.02</td><td>1.33</td><td>14.44</td><td>3.75</td><td></td><td></td><td></td></tr><tr><td>0.6</td><td>22.22</td><td>1.41</td><td>13.68</td><td>4.22</td><td>1.0</td><td>23.57</td><td>1.27</td><td>17.18</td><td>3.77</td><td></td><td></td><td></td></tr><tr><td>0.8</td><td>22.04</td><td>1.42</td><td>13.34</td><td>4.08</td><td>1.2</td><td>22.04</td><td>1.25</td><td>14.09</td><td>3.79</td><td></td><td></td><td></td></tr><tr><td>1.0</td><td>21.93</td><td>1.44</td><td>12.98</td><td>4.18</td><td>1.4</td><td>23.59</td><td>1.26</td><td>13.99</td><td>3.79</td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The axial ratio versus dynamical time for $n = 4$; $r_a = 100$-HR$_1$ (top-left panel), $r_a = 100$-IR$_1$ (top-right panel), $r_a = 1.5$-HR$_1$ (bottom-left panel) and $r_a = 1.5$-IR$_1$ (bottom-right panel). The amount of rotation is represented by using different kinds of lines. Solid lines: $τ = 0$ (HR) or $k = 0$ (IR), dotted lines: $τ = 0.3$ (HR) or $k = 0.6$ (IR), and dashed lines: $τ = 0.8$ (HR) or $k = 1.4$ (IR).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The evolution of the virial ratio for the initial systems defined respectively by $n = 4$, $r_a = 1.5$, and $n = 4$, $r_a = 100$ which have acquired their rotation according to the procedure HR$_2$ with $τ = 30$ per cent (thin line) and IR$_2$ with $k = 1.4$ (bold line). A similar evolution for the virial ratio (remaining very close to 1) is obtained for all runs.</figcaption></figure>

method, the stationary OM distribution function is modified by a positive definite and time-independent transformation. The resulting DF is then always stationary and positive definite. Moreover, as the modified systems are always spherical (no modification on positions have been performed), the new DF depends only on isolating integrals of motion, the energy and the squared angular momentum (Perez & Aly 1996). It is therefore a stationary solution of the collisionless Boltzmann-Poisson system. We have also verified that the resulting systems are always virialized, as confirmed by Fig. 4.

Let us first consider rotating systems with high values of $r_a$. We recall that non-rotating OM systems with the same anisotropy parameter are stable (Perez et al. 1996). The dynamical evolution obtained for such systems in our present simulations (see Fig. 5, top panels) allows us to distinguish two classes of behaviour. When $r_a$ is large and the rate of rotation stays modest (typically $\mu < 10$ per cent), we find that systems remain stable and spherical ($a_1 \approx a_2$, $f \approx 1$). However, when $r_a$ is large and the rate of rotation becomes
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The axial ratios versus dynamical time for *n* = 4; *r*<sub>a</sub> = 100-HR<sub>2</sub> (top-left panel), *r*<sub>a</sub> = 100-IR<sub>2</sub> (top-right panel), *r*<sub>a</sub> = 1.5-HR<sub>2</sub> (bottom-left panel) and *r*<sub>a</sub> = 1.5-IR<sub>2</sub> (bottom-right panel). The amount of rotation is represented by using different kinds of lines: solid lines: τ = 0(HR) or *k* = 0 (IR), dotted lines: τ = 0.3(HR) or *k* = 0.6 (IR), and dashed lines: τ = 0.8(HR) or *k* = 1.4 (IR).</figcaption></figure>

<table><caption>**Table 2.** Evolution of the stability parameter for systems with *n* = 4 and put in rotation by using the second method.</caption><thead><tr><th colspan="4">HR<sub>2</sub><br>r<sub>a</sub> = 100</th><th colspan="4">IR<sub>2</sub><br>r<sub>a</sub> = 1.5</th><th colspan="4">IR<sub>2</sub><br>r<sub>a</sub> = 100</th></tr><tr><th>τ</th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>k</th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th><th>k</th><th>P<sub>ε</sub></th><th>P<sub>λ</sub></th></tr></thead><tbody><tr><td>0.0</td><td>22.42</td><td>1.39</td><td>13.84</td><td>4.07</td><td>0.4</td><td>23.50</td><td>1.81</td><td>17.93</td><td>3.51</td><td>0.0</td><td>22.42</td><td>1.39</td></tr><tr><td>0.2</td><td>21.31</td><td>1.94</td><td>13.43</td><td>4.08</td><td>0.6</td><td>21.68</td><td>2.17</td><td>13.34</td><td>3.57</td><td>0.2</td><td>21.31</td><td>1.94</td></tr><tr><td>0.4</td><td>19.21</td><td>2.92</td><td>12.32</td><td>4.10</td><td>0.8</td><td>18.43</td><td>2.81</td><td>11.67</td><td>3.67</td><td>0.4</td><td>19.21</td><td>2.92</td></tr><tr><td>0.6</td><td>16.96</td><td>3.55</td><td>10.46</td><td>4.17</td><td>1.0</td><td>16.56</td><td>3.17</td><td>10.55</td><td>3.75</td><td>0.6</td><td>16.96</td><td>3.55</td></tr><tr><td>0.8</td><td>14.08</td><td>4.04</td><td>8.97</td><td>4.22</td><td>1.2</td><td>15.23</td><td>3.62</td><td>9.39</td><td>3.80</td><td>0.8</td><td>14.08</td><td>4.04</td></tr><tr><td>1.0</td><td>11.24</td><td>4.58</td><td>7.40</td><td>4.46</td><td>1.4</td><td>13.85</td><td>3.97</td><td>8.64</td><td>3.96</td><td>1.0</td><td>11.24</td><td>4.58</td></tr></tbody></table>

important (typically $\mu > 10$ per cent), we find that initially spherical systems develop a very soft bar-like instability ($a_2 \approx 1$, $a_1 \approx 0.85$, $f \approx 1.3$).

Systems with small $r_a$-values have a different behaviour, which depends on whether the rotational motion has been introduced by using a homogeneous or an inhomogeneous procedure. In the first case (HR<sub>2</sub> procedure), Fig. 5 (bottom-left panel) shows that systems which are radial orbit-unstable without rotation (e.g., $a_1 \approx 0.55$, $a_2 \approx 1$, $f \gg 1$) (solid line), become quasi-spheroidal ($a_1 \approx 0.85$, $a_2 \approx 1.1$, $f \approx 1.3$) when they have a modest rotation motion ($\mu < 10$ per cent) (dotted line). However, when rotation is important (typically $\mu > 10$ per cent), such systems develop a disc-like instability ($a_1 \approx 1$, $a_2 \approx 1.25$, $f \approx 0$) (dashed line). In the second case (IR<sub>2</sub> procedure)(bottom-right panel), the fact that rotation has been introduced only in a central region of the system prevents one from obtaining quasi-spherical systems, and therefore the radial-orbit instability persists ($a_1 \approx 0.65$, $a_2 \approx 1.1$, $f \approx 3.5$) for systems with a modest amount of rotation ($n = 4$; IR<sub>2</sub>;

$r_a = 1.5$; $k = 0.6$) (dotted line). When the amount of rotation is high enough ($\mu > 10$ per cent) a disc-like structure appears ($a_1 \approx 1$, $a_2 \approx 1.25$, $f \approx 0$), as in the HR<sub>2</sub> procedure. The evolution of the axial ratio for the system ($n = 4$; IR<sub>2</sub>; $r_a = 1.5$; $k = 1.4$) is represented by dashed line. These results hold whatever the *n* parameter value is. In practice, we have performed numerical simulations for three values of *n* ($n = 3.5, 4$ and $4.5$).

Are the stability parameters $P_\epsilon$ and $P_\lambda$ still discriminating for such rotating systems? We can see in Table 2 that non-rotating stable systems ($r_a = 100$) are predicted to remain stable when the second rotation method is applied ($P_\epsilon$ stays smaller than 20 and $P_\lambda$ stays larger than 2.5) whatever the $\tau$ and $k$ parameters values are. On the other hand, non-rotating radial-orbit unstable systems ($r_a = 1.5$) are predicted to become stable when sufficient method 2 rotation is applied, from $\tau = 0.4$ and $k = 0.8$ which correspond in Figs 1 and 2 to typically $\mu \approx 5$ per cent. Actually, in this case $P_\epsilon$ becomes less than 20, and $P_\lambda$ becomes greater than 2.5. Consequently, these stability parameters which have been constructed to
predict the stability of non-rotating spherical systems are not relevant as soon as the quantity of kinetic rotation energy becomes large, typically $\mu \ge 10$ per cent.

## 5 PHYSICAL INTERPRETATION AND CONCLUSIONS

The rotational properties of collapsed systems depend to a large extent on the amount of angular moment before the collapse. In order to study in a realistic way the importance of rotation for the dynamics of self-gravitating systems, it is necessary either to attempt an analytical approach, or to perform a complete numerical study modelling the collapse and relaxation phases prior to the two-body relaxation phase. However, although the collapse of a system can be studied by using the introduced amount of rotational kinetic energy as parameter, it is difficult to extract general conclusions from this kind of experiments. As a matter of fact, in this way the post-collapse physical features of the object cannot be completely controlled, and hence it can then be difficult to study with these methods the influence of the rotation on post-collapse systems.

This justifies the method that we have used in this paper. Starting from virialized systems with exactly known dynamical properties, we can study the influence of rotation by controlling its features. If the initial systems cover a wide variety of physical properties, and if the methods to introduce the rotation preserve certain fundamental features of these systems (invariance of mean energy, conservation or controlled modification of the distribution function), the numerical study will then be able to be used as a model to extract some general conclusions. As a matter of fact, our simulations start from a wide variety of initial conditions fully controllable through the dependence on the two parameters $n$ and $r_a$. On the other hand, the techniques used to introduce rotation to the systems preserve (as explained in Section 3) the properties that ensure that our models are spherical stationary solutions of the collisionless Boltzmann-Poisson system.

The main properties found in our study are the following.

(1) There do not exist spherical self-gravitating systems in 'fast' rotation. Our simulations show in fact that, when $\mu \ge 10$ per cent, the systems do not remain spherical but become lengthened along one or two axes depending on whether they are isotropic or anisotropic respectively, when they do not have a rotational motion.

(2) Rotation (in our case $HR_2$ and $IR_2$) can allow for a reorganization of systems in velocity-space able to modify their dynamical behaviour. We have, in fact, shown that a moderate rotation (typically $0 \le \mu \le 10$ per cent) can stabilize and confer a quasi-spherical structure to systems that, when they are not rotating, suffer a radial-orbit instability. Therefore there can exist rotating spherical self-gravitating systems. This is the case for our models with $r_a \ge 2$ and $\mu \le 10$ per cent.

(3) We have finally found that the stability parameters introduced in (Perez et al. 1996) remain valid as long as $\mu \le 10$ per cent.

## REFERENCES

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

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

* Antonov V. A., 1973, in Omarov G. B., ed., The Dynamics of Galaxies and Star Clusters. Alma Ata: Nauka 1973 (in Russian), and also in de Zeeuw T., ed., Proc. IAU Symp. 127, Structure and Dynamics of Elliptical Galaxies. Reidel, Dordrecht, p. 549

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

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

* Buta R., Crocker D. A., Elmegreen B. G., eds, 1996, IAU Colloq. 157, ASP Conf. Ser. Vol. 191, Barred Galaxies. Astron. Soc. Pac., San Francisco

* Calot G., 1973, Cours de Statistiques Descriptives. Dunod, Paris

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

* Goodwin S. P., 1997, MNRAS, 286, 139

* Lynden-Bell D., 1962, MNRAS, 124, 1

* Merritt D., 1985a, AJ, 90, 1027

* Merritt D., 1985b, MNRAS, 214, 25P

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

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

* Palmer P. L., Papaloizou J., 1987, MNRAS, 224, 1043

* Palmer P. L., 1994, Stability of Collisionless Stellar Systems. Kluwer Academic Publishers, Dordrecht

* Papaloizou J., Palmer P. L., Allen A. J., 1991, MNRAS, 253, 129

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

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

* Stavena A., Spassova N., Golev V. 1996, A&AS, 116, 447

## APPENDIX A: GENERATION OF INITIAL CONDITIONS

### A1 The initial positions for the particles

Let us consider a density $\rho_{iso}$ given by the polytropic relation
$\rho_{iso} = c_n \psi_{iso}^n$, where

$$c_n = \frac{(2\pi)^n \Gamma(n - 1/2)}{\Gamma(n + 1)}, \quad \Gamma(n + 1) := \int_0^\infty x^n e^{-x} dx,$$

and $\psi_{iso}$ is the solution of the Lame-Emden equation

$$\frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{d\psi}{dr} \right) + 4\pi G c_n \psi^n = 0.$$

This isotropic model is then deformed according to

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

where the anisotropic radius $r_a$ is a parameter which controls the deformation. The polytropic index $n$ is chosen in the range $[0.5, 5]$ in order to the system admits a finite mass $M(<r) = 4\pi \int_0^r r^2 \rho_{\text{ani}}(r) dr$. The total mass of the system is normalized to unity, and we then compute for a large set of particles ($1 \le i \le N$) from the inverse function of $M(x)$ and $\sin(x)$, the components

$$r_i = r_{\text{max}} M^{-1}(x),$$

$$\theta_i = 2 \arcsin(\sqrt{x}),$$

$$\phi_i = 2\pi x,$$

where $x$ is a uniform random variable in the range $[0, 1]$. The size of the system $r_{\text{max}}$ is chosen such that $\rho_{\text{ani}}(r > r_{\text{max}}) < 10^{-5}$.

### A2 The initial velocities for the particles

Let us consider the velocity components in spherical coordinates $(v_r, v_\theta, v_\phi)$; we compute $v_t = \sqrt{v_\theta^2 + v_\phi^2}$ and $\alpha = \arctan(v_\theta/v_\phi)$. The probability $\rho(\Gamma)$ of finding a particle in a volume $d\Gamma := dr d\theta d\phi dv_r dv_\theta d\alpha$ of the phase-space is defined by the
DF of the system

$$p(\Gamma)d\Gamma = \frac{1}{N}f(\Gamma)r^2dr \sin\theta d\theta d\phi v_t dv_t dv_t d\alpha. \quad (A1)$$

In the OM model, DF is a function only of the Q variable:

$$f = f(Q) \quad \text{with} \quad Q = \frac{1}{2}v_r^2 + \frac{1}{2}\left(1 + \frac{r^2}{r_a^2}\right)v_t^2 + \psi(r); \quad (A2)$$

equation (A1) then reduces to

$$\int p(\Gamma)d\Gamma = \int \frac{8\pi^2 r^2 v_t}{N} f(r, v_r, v_t) dr \, dv_t \, dv_r, \quad (A3)$$

where $r$, $v_r$ and $v_t$ are dependent random variables. In order to continue the integration of equation (A3), we introduce the variables $R$ and $\beta$, defined as follows

$$v_r = R \cos \beta,$$

$$v_t \sqrt{1 + \frac{r^2}{r_a^2}} = R \sin \beta.$$

We then get

$$f(r, v_r, v_t) \frac{8\pi r^2 v_t}{N} = f\left(\frac{R^2}{2} + \psi(r)\right) \frac{8\pi^2 r^2 R^2 \sin \beta}{N\left(1 + \frac{r^2}{r_a^2}\right)} dr dR d\beta.$$

We see from the previous expression that the random variable $\beta$ is independent of $r$ and $R$, and we have

$$p(\beta)d\beta = \frac{\sin \beta}{2} d\beta$$

$$p(r, R)drdR = \frac{16\pi^2 r^2 R^2}{1 + \frac{r^2}{r_a^2}} f\left[\frac{R^2}{2} + \psi(r)\right] drdR.$$

The conditional probability of finding a particle with a velocity defined by $R$ at a given distance $r_0$ is then

$$p(R|r_0) = \frac{p(r_0, R)}{p(r_0)} \quad \text{with} \quad p(r_0) = \frac{4\pi r^2 \rho(r_0)}{M(= 1)}, \quad (A4)$$

and finally

$$\begin{aligned} P(R|r=r_o) &= \frac{4\pi}{\rho(r_o)\left(1+\frac{r_o^2}{r_a^2}\right)} \int_0^R R'^2 f\left[\frac{R'^2}{2} + \psi(r_o)\right] dR' \\ &= \frac{2\pi}{\rho(r_o)\left(1+\frac{r_o^2}{r_a^2}\right)} \int_{\psi(r_o)}^{\psi(r_o)+\frac{r_o^2}{2}} \frac{f(Q)dQ}{\sqrt{2[Q-\psi(r_o)]}}. \end{aligned} \quad (A5)$$

We are now able to assign a velocity for each particle for which the position have been previously determined:

$$R_i = P^{-1}(x|r=r_o),$$

$$\beta_i = 2 \sin^{-1}(\sqrt{x}),$$

$$\alpha_i = 2\pi x$$

where $x$ is a uniform random variable in the range $[0, 1]$, and $P^{-1}$ is the inverse function of the probability defined by equation (A5). Finally,

$$v_{r i} = R_i \cos \beta_i,$$

$$v_{\theta i} = \sqrt{1 + \frac{r_o^2}{r_a^2} R_i \sin \beta_i \cos \alpha_i},$$

$$v_{\phi i} = \sqrt{1 + \frac{r_o^2}{r_a^2} R_i \sin \beta_i \sin \alpha_i}.$$

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