# Stability of spherical stellar systems I: Analytical results

Jérôme Perez, J.-J. Aly

► To cite this version:

Jérôme Perez, J.-J. Aly. Stability of spherical stellar systems I: Analytical results. Monthly Notices of the Royal Astronomical Society, 1996, 280 (3), pp.689-699. 10.1093/mnras/280.3.689. hal-01010751

HAL Id: hal-01010751

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

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
# Stability of spherical stellar systems – I. Analytical results

Jérôme Perez<fnref n="1" />,<fnref n="2" /> and Jean-Jacques Aly<fnref n="1" />

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

## ABSTRACT

The so-called 'symplectic method' is used for studying the linear stability of a self-gravitating collisionless stellar system, in which the particles are also subjected to an external potential. The system is steady and spherically symmetric, and its distribution function $f_0$ thus depends only on the energy $E$ and the squared angular momentum $L^2$ of a particle. Assuming that $\partial f_0 / \partial E < 0$, it is first shown that stability holds with respect to all the spherical perturbations – a statement which turns out also to be valid for a rotating spherical system. Thus it is proven that the energy of an arbitrary aspherical perturbation associated with a 'preserving generator' $\delta g_1$ [i.e., one satisfying $\partial f_0 / \partial L^2 \{\delta g_1, L^2\} = 0$] is always positive if $\partial f_0 / \partial L^2 \le 0$ and the external mass density is a decreasing function of the distance $r$ to the centre. This implies in particular (under the latter condition) the stability of an isotropic system with respect to all the perturbations.

Some new remarks on the relation between the symmetry of the system and the form of $f_0$ are also reported. It is argued, in particular, that a system with a distribution function of the form $f_0=f_0(E, L^2)$ is necessarily spherically symmetric.

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

## 1 INTRODUCTION

The analytical study of the linear stability of a steady collisionless stellar system is a quite formidable task which has great difficulties in its progression since the pioneering work of Antonov (Antonov 1961, 1973). The methods of investigation that have been used up to now fall into two categories: the normal mode approach and the energetic approach. The first of these is presented in a systematic way in the most recent monograph by Palmer (Palmer 1994, see also Fridman & Polyachenko 1984, and references in these two books). It consists of a derivation, from the linearized equations of motion (Vlasov-Poisson's system, hereafter VP) of a dispersion relation for the eigenfrequencies of the system, and an attempt to extract from it as much information as possible on the nature of these numbers. An instability is assumed to be present if at least one of the eigenfrequencies has a non-zero imaginary part, while stability holds if they are all real. The second method consists of constructing a quadratic functional $W$ over a set $\mathcal{L}$ of admissible test functions that satisfy some of the constraints actually fulfilled by any solution of VP. The system is found to be stable if $W$ keeps the same sign over $\mathcal{L}$. Actually there

are several types of energy method, which correspond to different choices for $W$ and $\mathcal{L}$. The simplest one is the 'thermodynamic method' (Lynden-Bell & Sannit 1969; Ipser 1974; Ipser & Horwitz 1979), in which only a very few constraints are taken into account. On the contrary, all of the VP constraints are retained in the more elaborate 'energy principle' (Antonov 1962; Kulsrud & Mark 1970; Kandrup & Sygnet 1985), as well as in the 'symplectic method' (Bartholomew 1971; Kandrup 1990, 1991; Perez 1995), the latter differing from the former by a restriction imposed on the perturbations, which are taken to be 'symplectic', i.e., they are generated by the infinitesimal canonical transforms acting on a steady state.

Our aim here is to reinvestigate the linear stability of a spherically symmetric system by fully working in the symplectic framework, with the hope of extending – or at least clarifying – some of the results which have been obtained thus far by the other energy methods. Our paper thus presents an interesting illustration of the general principles in Perez (1995). At the same time, it also provides the theoretical background for the following paper of this series (Perez et al. 1996, hereafter Paper II), in which the problem is tackled from a numerical point of view. It should be

<footnote n="1">¹CEA, DSM/DAPNIA, Service d'Astrophysique (CNRS URA 2052), CE Saclay, 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="2">²ETCA/CREA, 16 bis avenue Prieur de la côte d'or, 94114 Arcueil Cedex, France</footnote>
690 J. Perez and J.-J. Aly

stressed from the very beginning that it has not thus far been proved that the method used here will be able to reproduce all the results based on the normal mode analysis reported in Palmer (1994). The latter, however, are sometimes based on approximations which appear quite difficult to justify rigorously, while ours are derived without the help of any approximation.

Our plan is as follows. In Section 2, we state precisely our assumptions (which include the possible presence of an external potential, generally not considered by previous authors), and we develop some new considerations on the relationship between the symmetry of the steady states we are interested in and the form of their distribution functions $f_0$. In Section 3, we recall the definition of a linear symplectic perturbation, the expression of its energy (which we cast into a particularly useful form) and the relationship of the latter to the stability properties of the system. Thus we consider in detail stability with respect to spherical (Section 4) and aspherical (Section 5) perturbations, respectively, assuming that $f_0$ is a decreasing function of the energy. Some extensions of our results to spherical rotating systems (Lydden-Bell 1960), and to systems with non-monotonic $f_0$, are reported in Section 6.

## 2 EQUILIBRIUM

### 2.1 Notations and equations

We consider a system consisting of a large number of gravitating particles interacting together, and also subjected to an 'external' potential $\Phi_e$. We denote as $x$ and $v$, respectively, the position and the velocity of a particle with respect to a Galilean frame (O; $\hat{x}, \hat{y}, \hat{z}$), and as $\xi := (x, v)$ the corresponding point in the phase space $R^6$. Most often, we shall represent $\xi$ either by its standard spherical coordinates $(r, \theta, \phi, v_r, v_\theta, v_\phi)$ (with $\theta$ measured from the z-axis), or by its 'canonical' coordinates $(r, \theta, \phi, p_r, p_\theta, p_\phi)$, with $p_r=v_r, p_\theta=rv_\theta, p_\phi=r \sin \theta v_\phi$ being the 'conjugated' momenta of $(r, \theta, \phi)$ (e.g., Landau & Lifshitz 1966).

We assume the following.

(a) The external potential is spherically symmetric about O, i.e., $\Phi_e = \Phi_e(r)$. In most practical applications, $\Phi_e$ will represent the potential either of a central massive object or of a massive halo. For our theoretical purpose, however, we do not need to fix its form precisely, and we simply consider it as being created by some matter distributed with the given smooth density $\rho_e(r)$. The total respective masses $M_e$ of that distribution and M of the system are finite.

(b) The statistical state of the system is described at each time t by a one-particle distribution function $\mathcal{F}(\xi, t, m)$, with $\mathcal{F}(\xi, t, m) d\xi dm$ representing the number of particles of mass between m and $m + dm$ contained in the phase-space volume $d\xi$ around $\xi$.

(c) Collisions between the particles are negligible. It is then sufficient to consider the mass-averaged distribution function

$$f(\xi, t) := \int m \mathcal{F}(\xi, t, m) \, dm. \qquad (1)$$

f solves the VP system

$$\frac{\partial f}{\partial t} + \nu \cdot \nabla_x f - \nabla_x \Psi \cdot \nabla_x f = \frac{\partial f}{\partial t} + \{f, h\} = 0, \qquad (2)$$

$$\nabla_x^2 \Psi = 4\pi G (\rho + \rho_e) := 4\pi G \left( \int f d\nu + \rho_e \right), \qquad (3)$$

$$\Psi = \lim_{r \to \infty} O(r^{-1}) \quad \text{and} \quad \lim_{(r, v) \to \infty} f(x, v) = 0, \qquad (4)$$

with the decrease of f at infinity supposed to be sufficiently fast. Here, $\Psi := \Phi + \Phi_e$ is the total gravitational potential created by the particles and the external masses.

$$h(x, v, t) := \frac{v^2}{2} + \Psi(x, t) \qquad (5)$$

is the one-particle Hamiltonian, and {$, \cdot, \cdot$} denotes the Poisson bracket, defined by

$$\{f_1, f_2\} := \nabla_x f_1 \cdot \nabla_x f_2 - \nabla_x f_1 \cdot \nabla_x f_2. \qquad (6)$$

Note that there is no explicit dependence on the mass of the particles appearing in equations (2) and (3).

Instead of using both (2) and (3) we may use only equation (2), in which we have substituted for $\Phi$ the unique solution to equations (3) and (4). The latter can be written in integral form as $\Phi = K[f]$, with the operator K being quite generally defined by

$$K[f](x) := -G \int \frac{f(\xi')}{|x - x'|} d\xi'. \qquad (7)$$

### 2.2 Spherically symmetric steady states

We are interested here in a steady state which is spherically symmetric with respect to O and, in particular, occupies a sphere of radius $R (\le \infty)$ in the physical space. Hence $f_0$ – assumed from now on to be a sufficiently smooth function – and the associated potential $\Phi_0 = K[f]$ are left unchanged by any rotation of $R^3$ of centre O, which implies that they can depend only on the three geometrical invariants $r, v$ and $x \cdot v = rv$, i.e., they are necessarily of the form $f_0 = f_0(r, v, rv)$ and $\Phi_0 = \Phi_0(r)$, respectively.

Let us introduce the two quantities (where we set $v_r := v - v_r \hat{r}$)

$$E := h_0 := \frac{v^2}{2} + \Psi_0(r) := \frac{v_r^2}{2} + \frac{v_\theta^2}{2} + \Phi_0(r) + \Phi_e(r) \qquad (8)$$

and

$$L^2 := r^2 v_r^2 := r^2 (v_\theta^2 + v_\phi^2) = p_\theta^2 + \frac{p_\phi^2}{\sin^2 \theta}, \qquad (9)$$

which represent the energy of a particle in the potential $\Psi_0$ and the squared modulus of its angular momentum with respect to O. Clearly, we can express v in terms of $r, v_r$ and $L^2$, and thus we may write

$$f_0 = f_0(r, v_r, L^2). \qquad (10)$$
Inserting the expression into the time-independent version of equation (2), we obtain

$$v_r \frac{\partial f_0}{\partial r} + \left( \frac{L^2}{r^2} - \Psi'_0 \right) \frac{\partial f_0}{\partial v_r} = 0. \quad (11)$$

From the theory of first-order partial differential equations, we know that the general solution of equation (11) is a function of two independent integrals of its characteristic system

$$\frac{dr}{v_r} = \frac{dv_r}{L^2/r^2 - \Psi'_0}. \quad (12)$$

Two such integrals are obvious here: $L^2$ (as $dL^2$ does not appear) and $E(r, v_r, L^2)$. We can therefore conclude at once that

$$f_0 = f_0(E, L^2). \quad (13)$$

The particular case $f_0=f_0(E)$ corresponds to an isotropic equilibrium.

The general form of $f_0$ in equation (13) is well known, but it is usually derived (see, e.g., Dejonghe 1986, Binney & Tremaine 1987, Palmer 1994) by appealing to the so-called strong Jeans theorem, which deals with the isolating integrals of motion. As this theorem is far from being easy both to derive and to apply, we felt that it was a useful task to provide the reader with a completely straightforward proof. In particular, the proof presented here shows most clearly the basic reasons why $f_0$ depends on only two variables: symmetry requires that $f_0$ depends on only three variables, and the steady VP reduces this number by one unit.

Eventually, the potential $\Psi_0$ appearing in $E$ is a solution of

$$\nabla_x^2 \Psi_0 = \frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{d\Psi_0}{dr} \right) = 4\pi G (\rho_0 + \rho_c) = 4\pi G \times \left[ \int f_0(E, L^2) \, d\nu + \rho_c \right], \quad (14)$$

$$\Psi_0 = r \to \infty O(r^{-1}). \quad (15)$$

An immediate consequence (which will prove quite important hereafter) of these equations is that $\Psi_0(r)$ is a strictly increasing negative function (as $\rho_0 \ge 0$ and $\rho_c \ge 0$). $\Phi_0$ and $\Phi_c$ (if $\rho_c \ne 0$) obviously have the same property.

The characteristics of the steady-state VP describe the motion of the individual particles in the (self-consistent) potential $\Psi_0$. In particular, if we substitute $dr/dt$ for $v_r$, equation (12) and its integral $E$ determine the radial motion of a particle of energy $E$ and angular momentum modulus $L$, which can be described as a one-dimensional motion in the effective potential $\Psi_{\text{eff}} := \Psi_0 + L^2/2r^2$. When $r$ increases from 0 to $+\infty$, the latter decreases from $+\infty$ to a negative minimum, reached for $r_0(L^2)$, and then increases to 0.<fnref n="1" /> As the system is confined inside a sphere of radius $R$ and $v_r^2 \ge 0$, we have, necessarily,

$$E^{-}(L^2) := \Psi_{\text{eff}}[r_0(L^2)] = \Psi_0[r_0(L^2)] + \frac{L^2}{2[r_0(L^2)]^2} \le E \le \Psi_0(R) \\ + \frac{L^2}{2R^2} =: E^{+}(L^2) \le 0 \quad (16)$$

[with $E^+(L^2)=0$ if $R=\infty$; for an isotropic state, we can even take $E^+(L^2)=\Psi_0(R)$, because $v_r=0$ implies that $v=0$]: the particle thus bounces between the two solutions $r^-(E, L^2)$ and $r^+(E, L^2)$ [$r^- \le r_0 \le r^+ \le R$] to the equation $\Psi_{\text{eff}}(r)=E$. Equation (16) imposes $L^2$ to be bounded from above by the smallest solution $L_{\text{max}}^2$ of the equation $E^{-}(L^2)=E^{+}(L^2)$, or equivalently of $r_0(L_{\text{max}}^2)=R$ [the existence of that solution results immediately from the monotonic increase of $r_0(L^2)^{1/2}$]. We have eventually reached the conclusion that $f_0(E, L^2)$ must vanish outside the energy range $[E^{-}(L^2), E^{+}(L^2)]$, where $0 \le L^2 \le L_{\text{max}}^2$.

Hereafter, we shall denote as $\Omega_0$ the region of the phase space where $f_0 > 0$, and set

$$f_{0E} := \frac{\partial f_0}{\partial E} \quad \text{and} \quad f_{0L^2} := \frac{\partial f_0}{\partial L^2}. \quad (17)$$

Unless otherwise specified, it will be assumed that

$$f_{0E} < 0 \quad (18)$$

in $\Omega_0$.

## 2.3 Symmetry of $f_0$: converse theorems

In the previous subsection, we have shown that a steady spherically symmetric system has an $f_0$ of the form (13). It is then natural to address the converse question: if a steady state is characterized by a distribution function of the form $f_0(E, L^2)$, is it necessarily spherically symmetric?

Let us then assume that we have some steady state described by a sufficiently smooth distribution function $f_0(E, L^2)$ (condition 18 is not imposed here). Hence, $f_0$ does satisfy

$$\{f_0, E\} = f_{0E} \{E, E\} + f_{0L^2} \{L^2, E\} = -2f_{0L^2} v_r \cdot \nabla_x \Psi_0(\mathbf{x}) = 0, \quad (19)$$

while $u := -\Psi_0(\mathbf{x})$ is a positive solution (with no a priori prescribed symmetry) to the non-linear boundary value problem

<footnote n="1">¹This last statement is quoted very often in the literature, but without proof. To establish it, we remark that $\Psi'_{\text{eff}}=0$ implies $G(M, M_{cr}) = L^2/r$, where $(M, M_{cr})$ denotes the total mass in a sphere of radius $r$, and we have used the integral of equation (14) (Newton's theorem); as the left-hand side of that equation is a non-decreasing function of $r$, while its right-hand side is a decreasing one, it has one and only one solution $r_0(L^2)$. The conclusion thus follows by noting the obvious behaviour of $\Psi_{\text{eff}}$ when $r \to 0$ and when $r \to \infty$. It is also interesting to note that $r_0(L^2)$ increases monotonically with $L^2$ ($dr_0/dL^2 > 0$) from 0 to $\infty$, with $r_0 \propto L^{1/2}$ when $L^2 \to 0$, and $r_0 \propto L^2$ when $L^2 \to \infty$.</footnote>
692 J. Perez and J.-J. Aly

$$-\nabla_x^2 u(\mathbf{x}) = 4\pi G \left\{ \int f_0 \left[ \frac{\mathbf{v}^2}{2} - u(\mathbf{x}), r^2 v_t^2 \right] d\mathbf{v} + \rho_e(r) \right\} =: g[u(\mathbf{x}), r], \quad (20)$$

$$u = r \to \infty O(r^{-1}). \quad (21)$$

The function g defined in equation (20) is supposed to be sufficiently regular. Note that we have distinguished in the dependence of g on x a part which arises from the dependence of $f_0$ on E, and a part which is related to the dependence of both $f_0$ on $L^2$ and of $\rho_e$ on r (this separation being clearly possible without any ambiguity).

Consider first the case where $f_{0L^2} \neq 0$. Then, equation (19) implies that $\Psi_0(\mathbf{x}) = \Psi_0(r)$, (i.e., spherical symmetry holds) if, whenever $\rho_0(\mathbf{x}) > 0$, there is a value of $v$ such that $f_{0L^2}(\mathbf{x}, v) \neq 0$ [this situation is certainly met if $f_{0L^2} \neq 0$ in $\Omega_0$ – except maybe on a negligible set of points – which happens, for instance, if $f_0$ is of one of the standard forms $f_0(E \pm L^2/2r_\alpha^2)$, $f_1(E)L^{2\alpha}, \dots$ (Binney & Tremaine 1987)]. If the condition above is not satisfied, equation (19) only implies that $\Psi_0(\mathbf{x}) = \Psi_0(r)$ inside some spherical shells, but the general conclusion can still be reached by applying a continuation argument to equation (20) (note that both $\Psi_0$ and its normal derivative $\partial_r \Psi_0$ take constant values on the boundary of the shells).

The situation appears to be more intricate when $f_0$ does not depend explicitly on $L^2$ ($f_{0L^2} \equiv 0$). Fortunately, however, we may then use some recently obtained mathematical results which apply to the solutions to equations (20) and (21) under some technical assumptions (Gidas, Ni & Nirenberg 1981: the interested reader is referred to this paper for complete statements and proofs of the theorems we will use). Let us first assume that $\rho_e \equiv 0$ and $g(u) = O(u^2)$ near $u=0$, with $\alpha > 4$ [the last condition is satisfied in most cases of practical interest; for instance, it holds true when $f_0(E) = 0$ for $E^+ \le E$, with the constant $E^+ < 0$, in which case the system is contained in a bounded volume, and for Plummer's model, where $g(u) \propto u^5$ (Binney & Tremaine 1987)]. We are then within the conditions of theorem 1 of (Gidas et al. 1981), and we can assert at once that $u$ and $\Psi_0$ are necessarily spherically symmetric about the origin $O$ (possibly after a translation). The fact that $du/dr = -d\Psi_0/dr < 0$ for $r > 0$, also asserted by the theorem, results more directly here from g being non-negative (see Section 2.2). If $\rho_e \neq 0$, symmetry around $O$ is obtained as a consequence of theorem 1" of Gidas et al. (1981) (after some slight changes made possible by the particular form of our g); we need only assume that $\rho_e \le 0$, $\rho_0 \rho_e \neq 0$ at some point, and $g(u, r) \le cu^\alpha$ near $u=0$, with $c > 0$ and $\alpha > 4$.

We have thus eventually reached the conclusion that, under reasonable assumptions, a state characterized by $f_0(E, L^2)$ is spherically symmetric.

# 3 SYMPLECTIC ENERGY FUNCTIONAL AND LINEAR STABILITY

## 3.1 Symplectic perturbations

Let us now recall some of the basic facts about the symplectic approach to the stability of a steady state, referring the readers to Bartholomew (1971) and Kandrup (1990, 1991) for details. A linear 'symplectic' perturbation is

defined to be a first-order change of the distribution function which is of the form

$$f_1 = -\{f_0, g_1\}. \quad (22)$$

Here, $g_1$ is some arbitrary regular phase-space function – the so-called ‘generator’ – that it is convenient to consider formally as being complex valued. Clearly, a symplectic perturbation admits of an infinity of generators, two of them differing by some function $\tilde{g}$ commuting with $f_0$ ($\{f_0, \tilde{g}\} = 0$). We will say that $g_1$ is non-trivial if $f_1 \neq 0$.

To any $g_1$ are also associated a first-order variation of the potential,

$$\Phi_1 = K[f_1] = K[-\{f_0, g_1\}], \quad (23)$$

and a variation of the total energy which turns out to be second order in $g_1$, being given by

$$H^{(2)}[g_1] = -\frac{1}{2} \int \{f_0, g_1\} \{E, g_1^*\} d\xi - \frac{1}{8\pi G} \int |\nabla_x \Phi_1|^2 d\mathbf{x} \quad (24)$$

('*' denotes complex conjugation). It is worth noting the following.

(a) $H^{(2)}[g_1]$ is a conserved quantity if $g_1$ is taken to be a solution of the linearized equation of motion

$$\frac{\partial g_1}{\partial t} = \mathcal{L}[g_1] := \{E, g_1\} + K[\{f_0, g_1\}] \quad (25)$$

(in which case $f_1 = -\{f_0, g_1\}$ solves the linearized version of VP).

(b) If the system is translated as a whole by an infinitesimal vector $a$ ($g_1 = a \cdot v$ and $f_1 = -a \cdot \nabla_x f_0$), then $H^{(2)}[g_1] = 0$ in the absence of external potential (this is intuitively obvious, and easy to check formally). This type of 'neutral' displacement being of no interest, we shall eliminate it by imposing the admissible generators, whose set will be denoted by $\mathcal{G}$, to satisfy

$$\int f_1 \mathbf{x} \, d\xi = \int g_1 \nabla_x f_0 \, d\xi = 0 \quad (26)$$

in the case where $\Phi_e \equiv 0$, this condition impeding the centre of mass of the system to be moved by the perturbation.

## 3.2 Linear stability criteria

The stability properties of the steady state turn out to be strongly related to the values taken by $H^{(2)}$ on $\mathcal{G}$, and it is always an important first step when considering stability to classify as precisely as possible the generators $g_1$ according to the sign of that quantity. In fact, this study may even provide in some cases useful sufficient conditions of stability (Bartholomew 1971; Kandrup 1990, 1991).

(a) If

$$H^{(2)}[g_1] > 0 \quad (\text{or} < 0) \quad (27)$$

for any non-trivial $g_1$ of $\mathcal{G}$, then the system is linearly stable with respect to any perturbation. Unfortunately, this condi-
$$H^{(2)}[\bar{g}_1] = \frac{1}{2} \int (-f_{0E}) |\{E, \bar{g}_1\}|^2 d\xi - \frac{1}{8\pi G} \int |\nabla_x \bar{\Phi}_1|^2 dx, \quad (39)$$

$$H^{(2)}[\delta g_1] = -\frac{1}{2} \int \{f_0, \delta g_1\} \{E, \delta g_1^*\} d\xi \\ -\frac{1}{8\pi G} \int |\nabla_x \delta \Phi_1|^2 dx. \quad (40)$$

We now consider in turn each of the two terms appearing in the right-hand side of equation (38).

# 4 STABILITY WITH RESPECT TO SPHERICAL PERTURBATIONS

## 4.1 A further decomposition

For studying the term $H^{(2)}[\bar{g}_1]$, it is useful to effect one further decomposition by proceeding as follows. For any average $\bar{g}$, we denote as $\bar{g}^{+/-}$, respectively, its symmetric and antisymmetric parts (with respect to the transform $v_r \to -v_r$), defined by

$$\bar{g}(r, v, rv_r) = \frac{1}{2}[\bar{g}(r, v, rv_r) + \bar{g}(r, v, -rv_r)] + \frac{1}{2}[\bar{g}(r, v, rv_r) \\ - \bar{g}(r, v, -rv_r)] = :\bar{g}^{+}(r, v, rv_r) + \bar{g}^{-}(r, v, rv_r). \quad (41)$$

With

$$\bar{g}_1 = \bar{g}_1^+ + \bar{g}_1^-, \quad (42)$$

$$\bar{f}_1 = \bar{f}_1^+ + \bar{f}_1^-, \quad (43)$$

we thus have

$$\bar{f}_1^{\pm} = -\{f_0, \bar{g}_1^{-}\}, \quad (44)$$

$$\bar{\Phi}_1 = K[\{f_0, \bar{g}_1^{-}\}], \quad (45)$$

and, after some simple algebra, we obtain

$$H^{(2)}[\bar{g}_1] = H^{(2)}[\bar{g}_1^+] + H^{(2)}[\bar{g}_1^-], \quad (46)$$

with

$$H^{(2)}[\bar{g}_1^+] = \frac{1}{2} \int (-f_{0E}) |\{E, \bar{g}_1^+\}|^2 d\xi, \quad (47)$$

$$H^{(2)}[\bar{g}_1^-] = \frac{1}{2} \int (-f_{0E}) |\{E, \bar{g}_1^-\}|^2 d\xi - \frac{1}{8\pi G} \int |\nabla_x \bar{\Phi}_1|^2 dx. \quad (48)$$

The first part $H^{(2)}[\bar{g}_1^+]$ of $H^{(2)}[\bar{g}_1]$ is clearly non-negative, and thus we just need to consider in detail the second part $H^{(2)}[\bar{g}_1^-]$.

## 4.2 The sign of $H^{(2)}[\bar{g}_1^{-}]$

To study the sign of $H^{(2)}[\bar{g}_1^{-}]$, we closely follow the method of Sygnet et al. (1984), which originates in the work of Gillon, Doremus & Baumann (1976). However, in addition

tion may be satisfied only if $f_0=f_0(E)$ and $f_{0E} \le 0$ (actually, it may be fulfilled for a more general form of $f_0$ corresponding to a uniformly rotating system). Indeed, if $f_0$ depends explicitly on some other integral, both positive and negative energy perturbations always exist [Perez 1995: the main step in proving this result is the construction of two phase-space functions $\alpha$ and $\beta$ and a number $\epsilon$ such that the first term in the right-hand side of equation (24) be positive (respectively, negative) for the generator $g_1 := \alpha e^{i\beta\epsilon}$].

(b) Suppose that we cannot prove equation (27) for all the admissible $g_1$, but that we can identify in $\mathcal{G}$ a linear subspace $\mathcal{G}'$ that is 'closed' under evolution, which means that, if an arbitrary $g_1 \in \mathcal{G}'$ is taken as an initial condition for an evolution governed by equation (25), then the solution $g_1(t)$ stays in $\mathcal{G}'$ forever. Then the system is stable with respect to the corresponding restricted class of perturbations if equation (27) holds true for all the non-trivial elements of $\mathcal{G}'$.

## 3.3 Decomposition of $H^{(2)}$

For our spherically symmetric system, it will prove convenient to transform the expression for $H^{(2)}[g_1]$ into a different form. For that, we first introduce the following notations associated with the arbitrary phase function $g$ and the arbitrary rotation $\mathcal{R}$ of $\mathbb{R}^3$ about $O$:

$$g_{\mathcal{R}}(\xi) = g(\mathcal{R}x, \mathcal{R}v), \quad (28)$$

$$\bar{g}(\xi) = \frac{1}{4\pi} \int g_{\mathcal{R}}(\xi) d\mathcal{R}, \quad (29)$$

$$\delta g = g - \bar{g}. \quad (30)$$

Clearly, the average $\bar{g}$ is spherically symmetric, i.e., $\bar{g}(\mathcal{R}\xi) = \bar{g}(\xi)$ for any $\mathcal{R}$ and $\bar{g} = \bar{g}(r, v, rv_r)$. We also have the useful simple relations

$$\bar{\delta g} = \frac{1}{4\pi} \int (\delta g)_* d\mathcal{R} = 0, \quad (31)$$

and, for any $g$ and $g'$,

$$\{g_{\mathcal{R}}, g'_{\mathcal{R}}\}(\xi) = \{g, g'\}(\mathcal{R}\xi). \quad (32)$$

For a perturbation generated by $g_1$, we can thus write

$$g_1 = \bar{g}_1 + \delta g_1, \quad (33)$$

$$f_1 = \bar{f}_1 + \delta f_1, \quad (34)$$

$$\Phi_1 = \bar{\Phi}_1 + \delta \Phi_1. \quad (35)$$

Using the symmetry of $f_0$ and $E$, the rotational invariance of the Laplacian and the properties of the averaging process quoted above, we obtain immediately

$$\bar{f}_1 = -\{f_0, \bar{g}_1\} = -f_{0E}\{E, \bar{g}_1\} \quad \text{and} \quad \delta f_1 = -\{f_0, \delta g_1\}, \quad (36)$$

$$\bar{\Phi}_1 = K[\{f_0, \bar{g}_1\}] \quad \text{and} \quad \delta \Phi_1 = K[\{f_0, \delta g_1\}], \quad (37)$$

and after some straightforward algebra

$$H^{(2)}[g_1] = H^{(2)}[\bar{g}_1] + H^{(2)}[\delta g_1], \quad (38)$$

with, of course,
694 J. Perez and J.-J. Aly

to the fact that we have an external potential $\Phi_e$ and that we
do not use $(r, L^2, E)$ as independent variables, there are
some essential differences between their calculations and
ours, which we shall discuss below.

First, we transform the expression $H^{(2)}[\bar{g}_1^-]$ into a more
convenient form. We integrate once Poisson's equation

$$
\begin{align}
\nabla_x^2 \bar{\Phi}_1 &= \frac{1}{r^2} \frac{d}{dr} r^2 \frac{d\bar{\Phi}_1}{dr} = -4\pi G \int \{f_0, \bar{g}_1^-\} \, dv \nonumber \\
&= 4\pi G \frac{1}{r^2} \frac{d}{dr} \left( r^2 \int f_{0E} v_r \bar{g}_1^- \, dv \right), \tag{49}
\end{align}
$$

(where we have effected two integrations by part in the
second member to arrive at the right-hand side), which
gives

$$
\nabla \bar{\Phi}_1 = \frac{d\bar{\Phi}_1}{dr} \hat{r} = \left( 4\pi G \int f_{0E} v_r \bar{g}_1^- dv \right) \hat{r}. \quad (50)
$$

Injecting this result into the second term of the right-hand
side of equation (39) and using Schwartz's inequality, we
obtain

$$
\begin{align}
\frac{1}{8\pi G} \int |\nabla_x \bar{\Phi}_1|^2 \, dx &= 2\pi G \iint \left| \int f_{0E} v_r \bar{g}_1^- \, dv \right|^2 \, dx \nonumber \\
&\le 2\pi G \iint \left| \int (-f_{0E}) v_r^2 \, dv \right| \nonumber \\
&\qquad \times \left| \int (-f_{0E}) |\bar{g}_1^-|^2 \, dv \right| \, dx. \tag{51}
\end{align}
$$

Noticing that

$$
\rho_0 = \int f_0 \, dv = \int (-f_{0E}) v_r^2 \, dv \quad (52)
$$

(using $dv = 2\pi v_r dv$, $dv_t$ and integrating by parts with respect to $v_r$), we thus obtain

$$
H^{(2)}[\bar{g}_1^{-}] \ge \frac{1}{2} \int (-f_{0E}) (|\{E, \bar{g}_1^{-}\}|^2 - 4\pi G \rho_0 |\bar{g}_1^{-}|^2) \, d\xi. \quad (53)
$$

We now make the change of variable

$$
\bar{g}_1^- := rv_r \mu, \tag{54}
$$

with the new function $\mu$ still being regular for $rv_r=0$ as $\bar{g}_1^-$ ($r$, $v$, $rv_r=0$)=0 owing to the antisymmetry of $\bar{g}_1^-$ . Using the derivative property of the Poisson bracket, and the fact that the integral over the phase space of a Poisson bracket vanishes, we can rewrite equation (53) in the form

$$
\begin{equation}
\begin{split}
H^{(2)}[\bar{g}_1^{-}] \ge \frac{1}{2} \int & (-f_{0E}) [(rv_r)^2 |\{E, \mu\}|^2 + rv_r \{E, rv_r\} \{E, |\mu|^2\} \\
& + |\mu|^2 \{E, rv_r\}^2 - 4\pi G \rho_0 (rv_r |\mu|^2) ] \, d\xi
\end{split}
\tag{55}
\end{equation}
$$

$$
\begin{align}
&= \frac{1}{2} \int (-f_{0E}) [(rv_r)^2 |\{E, \mu\}|^2 + \{E, |\mu|^2 rv_r \{E, rv_r\}\} \nonumber \\
&\qquad - |\mu|^2 rv_r \{E, \{E, rv_r\}\} - |\mu|^2 \{E, rv_r\}^2 \nonumber \\
&\qquad + |\mu|^2 \{E, rv_r\}^2 - 4\pi G \rho_0 |rv_r, \mu|^2] \, d\xi \nonumber \\
&= \frac{1}{2} \int (-f_{0E}) [(rv_r)^2 |\{E, \mu\}|^2 - |\mu|^2 (rv_r \{E, \{E, rv_r\}\} \nonumber \\
&\qquad + 4\pi G \rho_0 (rv_r)^2)] \, d\xi. \tag{55}
\end{align}
$$

By a straightforward calculation, we obtain

$$
\begin{align}
\{E, \{E, rv_r\}\} &= -rv_r \left( \frac{d^2\Psi_0}{dr^2} + \frac{3}{r} \frac{d\Psi_0}{dr} \right) \nonumber \\
&= -rv_r \left[ 4\pi G (\rho_0 + \rho_e) + \frac{1}{r} \frac{d\Psi_0}{dr} \right], \tag{56}
\end{align}
$$

whence, injecting this expression into equation (55),

$$
H^{(2)}[\bar{g}_1^{-}] \ge \frac{1}{2} \int (-f_{0E})(rv_r)^2 \left[ |\{E, \mu\}|^2 + |\mu|^2 \right. \\
\left. \qquad \times \left( 4\pi G \rho_e + \frac{1}{r} \frac{d\Psi_0}{dr} \right) \right] d\xi \ge 0, \quad (57)
$$

the last inequality sign resulting from the negativeness of $f_{0E}$
and the positiveness of $\rho_e$ and $d\Psi_0/dr$.

Let us now assume that $H^{(2)}[\bar{g}_1^-] = 0$. Then $H^{(2)}[\bar{g}_1^+] = 0$, whence $\{f_0, \bar{g}_1^+\} = 0$; and $H^{(2)}[\bar{g}_1^-] = 0$, whence $f_{0E}\bar{g}_1^- = 0$, and $\{f_0, \bar{g}_1^-\} = 0$. Therefore $\bar{f}_1 = 0$, and we can eventually conclude that $H^{(2)}[\bar{g}_1^-] > 0$ for any $g_1$ such that $\bar{f}_1 = -\{f_0, \bar{g}_1^-\} \neq 0$.

We can summarize the results of this section in the form
of a theorem.

*Theorem 1.* If the steady state satisfies condition (18), then $H^{(2)}[\bar{g}_1^-] > 0$ for any $g_1$ with a non-trivial average $\bar{g}_1$.

4.3 Consequences for stability

It is quite obvious that the set $\mathcal{S}'$ of the generators under-
lying the spherically symmetric perturbations is 'closed with
respect to evolution' in the sense of Section 3.2, i.e., the
solution $g_1(t)$ to equation (25) generates a spherically sym-
metric perturbation if the initial condition $g_1(0)$ does. Then,
by the general statement (b) of Section 3.2, we can conclude
at once that any spherical steady state satisfying (18) is linearly
stable with respect to all the spherical perturbations.

4.4 Comparison with previous works

Of course, this result has been reported many times in the
literature for the case where $\Phi_e \equiv 0$, and it is then worth
explaining in a few words the differences between our proof
and previous ones. In fact, in the framework of the energy
methods, the aim consists, in any case, of showing the pos-
itiveness of some quadratic functional $W[q]$ over a set $\mathcal{Q}$ of
admissible functions $q$.
(a) In the first class of methods (Antonov 1961; Kulsrud & Mark 1970; Kandrup & Sygnet 1985), $\mathcal{Q}_1$ is taken to be the set of all the functions which vanish outside $\Omega_0$ and are antisymmetric in the change $v_r \to -v_r$, to which belongs the antisymmetric part $\bar{f}_1$ of any solution to the linearized VP. The relevant functional $W_1$ may be derived in a standard way (e.g., Laval, Mercier & Pellat 1965) from the equation of motion (second order in time!)

$$ \frac{1}{-f_{0E}} \frac{\partial^2 \bar{f}_1^-}{\partial t^2} = \mathcal{M} \bar{f}_1^- \qquad (58) $$

satisfied by $\bar{f}_1$, the self-adjoint character (with respect to the usual Hermitian product) of the operator $\mathcal{M}$ (Kandrup & Sygnet 1986) here being essential. It turns out that

$$ W_1[\bar{f}_1^-] = H^{(2)}[\bar{f}_1^-]. \qquad (59) $$

It appears that the 'energy principle' method and the symplectic one lead to formally identical functionals in the situation we are considering here. They can thus be treated by similar techniques, and the same results are obtained. However, the two methods are conceptually very different (Bartholomew 1971; Kandrup 1990, 1991). In the former, the derived 'energy functional' is closely related to the peculiar properties of the second-order equation one starts with (the existence of $W_1$ results from the self-adjointness of $\mathcal{M}$, which is just a consequence of the strong symmetries of the system), and it has no obvious physical meaning. The latter, on the contrary, is more transparent. It only deals with a first-order equation – the one satisfied by the generators $g_1$ of the symplectic perturbations (which form a restricted class of perturbations) – and $H^{(2)}[\bar{g}_1]$ is easily seen to actually be the energy of the perturbation, an interpretation which keeps its validity when one considers more general systems.

(b) In a second class of methods (Lynden-Bell & Sannit 1969; Ipser 1974; Ipser & Horwitz 1979), the relevant functional $W_2$ is the second variation of $E[f] + H_C[f]$. Here $E[f]$ is the energy of $f$ and

$$ H_C[f] := - \int C[f, L^2] \, d\xi := \iint_0^f E(f_0, L^2) \, df_0 \, d\xi, \qquad (60) $$

with $E(f_0, L^2)$ being the 'inverse' of the monotone function $f_0(E, L^2)$ (the results obtained by the authors quoted above are for isotropic systems, but there is no problem in generalizing them to the anisotropic case, as we do here). Formally, $W_2$ is easily shown to satisfy the relation

$$ W_2[-\{f_0, \bar{g}_1\}] = H^{(2)}[\bar{g}_1] \qquad (61) $$

for any function $g_1$. The admissible $\mathcal{Q}_2$ is taken to be the set of all the $\bar{f}_1$ satisfying the constraints

$$ \int \bar{f}_1 d\xi = 0 \quad \text{and} \quad \int \bar{f}_1 C'(f_0, L^2) d\xi = 0, \qquad (62) $$

which imply a conservation of mass and $H_C$ function to first order. The constraints which are retained from VP are thus quite weak, and, as a consequence, it has not been possible up to now to derive in this framework stability results as

complete as those furnished by the energy and the symplectic methods.

(c) The functional used by Sygnet et al. (1984) appears to be $W_2[\bar{f}_1]$, as in the previous method, but the set of admissible functions $\bar{f}_1$ is not clearly specified. The authors, however, effect without further justifications the change of dependant variable $\bar{f}_1 \to \mu$ defined by

$$ \bar{f}_1(r, E, L^2) = \left[ (rv_r)_{f_{0E}} \frac{\partial (rv_r \mu)}{\partial r} \right] (r, E, L^2) \qquad (63) $$

(their equation B3, written with our notation), where use is made of $(r, E, L^2)$ as independent variables [actually, this is a proper choice only in the region $(v_r > 0)$ or in the region $(v_r < 0)$]. However, such a change can be valid only for particular perturbations. First, it is easy to see that any $\bar{f}_1$ of the form equation (63) is symplectic, being generated by $\bar{g}_1(r, E, L^2) := (rv_r \mu)(r, E, L^2)$. Secondly, the new function $\mu$ introduced by this relation can be regular only if $\bar{g}_1=0$ for $rv_r=0$. In our work, on the contrary, these conditions are an explicit part of a coherent formalism: the symplectic nature of the perturbations is a basic element of the framework, while the vanishing of $\bar{g}_1$ for $rv_r=0$ is automatically ensured by restricting our attention to the antisymmetric part $\bar{g}_1^-$ of $\bar{g}_1$, the part of the functional related to $\bar{g}_1^+$ turning out to be trivially positive.

# 5 STABILITY WITH RESPECT TO ASPHERICAL PERTURBATIONS

## 5.1 A further decomposition

Using the relation

$$ \delta f_1 = - \{f_0, \delta g_1\} = -f_{0E} \{E, \delta g_1\} - f_{0L^2} \{L^2, g_1\}, \qquad (64) $$

we can write the quantity $H^{(2)}[\delta g_1]$ as a sum of two terms:

$$ H^{(2)}[\delta g_1] = H_1^{(2)}[\delta g_1] + H_2^{(2)}[\delta g_1], \qquad (65) $$

with

$$ H_1^{(2)}[\delta g_1] := \frac{1}{2} \int \frac{|\{f_0, \delta g_1\}|^2}{-f_{0E}} d\xi - \frac{1}{8\pi G} \int |\nabla_x \delta \Phi_1|^2 dx, \qquad (66) $$

$$ H_2^{(2)}[\delta g_1] := -\frac{1}{2} \int \frac{f_{0L^2}}{-f_{0E}} \{L^2, \delta g_1^*\} \{f_0, \delta g_1\} d\xi. \qquad (67) $$

## 5.2 Sign of $H_1^{(2)}[\delta g_1]$

We now set

$$ \delta f_1 = - \{f_0, \delta g_1\} := -f_{0E} \delta \Phi_1 + \delta \bar{f}_1, \qquad (68) $$

that is, we use the standard trick (introduced in the gravitational context by Lynden-Bell 1966) which consists of singling out the 'quasi-static part' in the variation $\delta f_1$ of the distribution function. Then, taking into account the relation

$$ \int \delta f_1 \delta \Phi_1 d\xi = - \frac{1}{4\pi G} \int |\nabla_x \delta \Phi_1|^2 dx, \qquad (69) $$
we can write

$$
\begin{align}
H_1^{(2)}[\delta g_1] &= \frac{1}{2} \int \frac{(\delta f_1)^2}{-f_{0E}} d\xi - \frac{1}{8\pi G} \int |\nabla_x \delta\Phi_1|^2 dx \nonumber \\
&= \frac{1}{2} \int \left[ \frac{(\delta f_1)^2}{-f_{0E}} - f_{0E}(\delta\Phi_1)^2 + (\delta f_1 + f_{0E}\delta\Phi_1)\delta\Phi_1 \right] \nonumber \\
&\quad \times d\xi - \frac{1}{8\pi G} \int |\nabla_x \delta\Phi_1|^2 dx \nonumber \\
&= \frac{1}{8\pi G} \int \left\{ |\nabla_x \delta\Phi_1|^2 - 4\pi G \left[ \int (-f_{0E}) dv \right] |\delta\Phi_1|^2 \right\} \nonumber \\
&\quad \times dx + \frac{1}{2} \int \frac{(\delta f_1)^2}{-f_{0E}} d\xi. \tag{70}
\end{align}
$$

Following Aly & Perez (1992), we now set

$$
\delta\Phi_1 := \Psi'_0 w. \tag{71}
$$

Then

$$
\begin{align}
H_1^{(2)}[\delta g_1] &= \frac{1}{2} \int \frac{(\delta \tilde{f}_1)^2}{-f_{0E}} d\xi + \frac{1}{8\pi G} \int \left\{ \Psi'_0^2 [ |\nabla_x w|^2 \right. \nonumber \\
&\quad \left. - 4\pi G \left( \int (-f_{0E}) dv \right) |w|^2 \right] + |w|^2 (\nabla_x \Psi'_0)^2 \nonumber \\
&\quad + \Psi'_0 \nabla_x \Psi'_0 \cdot \nabla_x |w|^2 \right\} dx \nonumber \\
&= \frac{1}{2} \int \frac{(\delta \tilde{f}_1)^2}{-f_{0E}} d\xi + \frac{1}{8\pi G} \int \left\{ (\Psi'_0)^2 \left[ |\nabla_x w|^2 \right. \right. \nonumber \\
&\quad \left. \left. - 4\pi G \left( \int (-f_{0E}) dv \right) |w|^2 \right] \right. \nonumber \\
&\quad \left. - |w|^2 \Psi'_0 \nabla_x^2 \Psi'_0 \right\} dx. \tag{72}
\end{align}
$$

The quantity $\nabla_x^2 \Psi'_0$ can be computed by differentiating equation (14) with respect to $r$, which gives

$$
\nabla_x^2 \Psi'_0 = 4\pi G (\rho'_0 + \rho'_c) + 2 \frac{\Psi'_0}{r^2} \quad (73)
$$

$$
= 4\pi G \left[ \int (\Psi'_0 f_{0E} + 2rv'_i f_{0L^2}) dv + \rho'_c \right] + 2 \frac{\Psi'_0}{r^2}.
$$

On the other hand, $w$ has zero average value on any spheri-
cal surface of centre $O$ and, then, by the so-called Wirtinger
inequality (Aly & Perez 1992), we have for any value of $r$

$$
\int \left( \left| \nabla_{xs} w \right|^2 - \frac{2}{r^2} \left| w \right|^2 \right) d\Omega \ge 0, \quad (74)
$$

with equality holding if and only if $w$ is of the form $w(r, \theta, \phi) = w_0(r) \hat{r} \cdot a$ for some constant vector $a$ (in equation 74, we have set $\nabla_{xs} := \nabla_x - \hat{r}\partial/\partial r$ and $d\Omega := \sin\theta d\theta d\phi$).

Using these results in equation (72), we eventually
obtain

$$
\begin{equation}
\begin{aligned}
H_1^{(2)}[\delta g_1] &= \frac{1}{2} \int \frac{(\delta \tilde{f}_1)^2}{-f_{0E}} d\xi + \frac{1}{8\pi G} \int \left\{ (\Psi'_0)^2 \left[ \left| \frac{\partial w}{\partial r} \right|^2 + |\nabla_{xs} w|^2 - 2 \frac{|w|^2}{r^2} \right] \right. \\
&\qquad \left. \times \left[ (-\rho'_c) + 2r \int v_i^2 (-f_{0L^2}) \, dv \right] \right\} dx \\
&\geq \frac{1}{2} \int \frac{(\delta \tilde{f}_1)^2}{-f_{0E}} d\xi + \frac{1}{8\pi G} \int \left\{ (\Psi'_0)^2 \left| \frac{\partial w}{\partial r} \right|^2 \right. \\
&\qquad \left. + 4\pi G \Psi'_0 \left[ (-\rho'_c) + 2r \int v_i^2 (-f_{0L^2}) \, dv \right] |w|^2 \right\} dx.
\end{aligned}
\tag{75}
\end{equation}
$$

Therefore, $H_1^{(2)}[\delta g_1] \ge 0$ for all $g_1$ if

$$
\rho'_{c} \le 0 \tag{76}
$$

and

$$
f_{0L^2} \le 0. \tag{77}
$$

Let us now assume that both conditions are satisfied, and that $H_1^{(2)}[\delta g_1] = 0$. Then, (i) $\delta \tilde{f}_1' = 0$; (ii) $\partial w / \partial r = 0$ and Wirtinger's inequality reduces to an equality, whence $w = \hat{r} \cdot a$; (iii) if either $\rho'_c \ne 0$ or $f_{0L^2} \ne 0$, we must also have $w = 0$ for some values of $r$ - then $w = 0$ everywhere, and $\delta f_1' = 0$; (iv) if $\rho'_c = 0$ (which implies $\rho_c = 0$ and $\Phi_c = 0$) and $f_{0L^2} = 0$, the previous results and the definition of $w$ imply $\delta f_1' = -f_{0E} \Psi'_0 \hat{r} \cdot a = -a \cdot \nabla_x f_0$, which means that the perturbation results simply from a rigid translation of the system of a vector $a$. But there are no such perturbations with $a \ne 0$ generated by the elements of $\mathcal{G}$ (see Section 3), and also in this case we have $w = 0$ and $\delta f_1' = 0$. We have thus proven eventually that, under conditions (76) and (77),

$$
H_1^{(2)}[\delta g_1] > 0 \quad \text{if} \quad \delta f_1' \neq 0. \tag{78}
$$

5.3 Sign of $H_2^{(2)}[\delta g_1]$

Clearly, the term $H_1^{(2)}[\delta g_1]$ can be of either sign, depending on the choice of the generator. We even know from the general theorem of Perez (1995), quoted in Section 3.2, that, by an appropriate choice of $g_1$, it can be made suffi- ciently negative to overcome the positiveness of $H_1^{(2)}[\bar{g}_1] + H_1^{(2)}[\delta g_1]$, which holds under conditions (76) and (77). Then, we could only attempt to classify the perturba- tions according to the sign of that quantity.

Even this more modest goal appears to be difficult to
reach, and we shall content ourselves here with the follow-
ing simple remark: there is a particular class of perturba-
tions, defined by the condition

$$
f_{0L^2}\{L^2, g_1\} = f_{0L^2}\{L^2, \delta g_1\} = 0, \quad (79)
$$

for which $H_1^{(2)}[\delta g_1]$ vanishes [note the use in equation (79) of the obvious relation $\{L^2, \bar{g}_1\} = 0$]. We feel that these ele-
ments of $\mathcal{G}$ – which thus commute with $L^2$ in the region of the phase-space where $f_{0L^2} \neq 0$ – play some role in the dynamics of the system, and then they deserve a particular name: we shall call them ‘preserving generators’, the associated perturbations being the ‘preserving perturbations’, and denote by $\mathcal{P}$ their set.

Clearly, $\mathcal{P}$ is a linear subspace of $\mathcal{G}$, and the general form of its elements (which do automatically satisfy condition 26) can be found as follows. Consider the part $\Omega_1$ of $\Omega_0$ where $f_{0L^2} \neq 0$. We first note that there are five obvious solutions to equation (79) in $\Omega_1$, namely $r, p_r, L_x = -\cot \theta \cos \phi_\phi - \sin \phi p_\theta, L_y = -\cot \theta \sin \phi_\phi + \cos \phi p_\theta$ and $L_z = p_\phi$. Moreover, they are linearly independent, which can be most easily checked by considering the Jacobian determinant $D(r, p_r, L_x, L_y, L_z)/D(r, \theta, \phi, p_r, p_\theta, p_\phi)$, which turns out to be of rank 5. As the general solution to equation (70) in $\Omega_1$ must depend on at most five integrals, we can thus conclude that any preserving generator is of the form

$$g_1(r, \theta, \phi, p_r, p_\theta, p_\phi) = g_1(r, p_r, p_\phi, -\cot \theta \cos \phi_\phi - \sin \phi p_\theta, -\cot \theta \sin \phi_\phi + \cos \phi p_\theta), \quad (80)$$

in $\Omega_1$, and arbitrary elsewhere. Of course, all the spherically symmetric $g_1$ belong to $\mathcal{P}$; but this set contains many other elements, as the spherically symmetric elements depend on three variables $[g_1 = g_1(r, p_r, L_x^2 + L_y^2 + L_z^2)]$, while the general preserving ones depend on five variables. Unfortunately, $\mathcal{P}$ is not ‘closed under evolution’ (see Section 3.2). Indeed, it results at once from the equation of evolution (25) that condition (79) is time invariant only if the potential $\Phi_1$ stays spherically symmetric, which can be the case only for spherically symmetric perturbations.

We can thus eventually summarize the results of Sections 5.2–5.3 in the following form.

**Theorem 2.** If the conditions (18), (76) and (77) are satisfied, then we have $H^{(2)}[g_1] > 0$ for all the non-trivial aspherical generators of $\mathcal{G}$ satisfying $f_{0L^2}\{L^2, \delta g_1\} = 0$.

## 5.4 Consequences for stability

From Theorem 2 and the general statements of Section 3.2, we can make the immediate conclusion that an isotropic system is stable with respect to any aspherical perturbations – and then is stable with respect to any perturbation, because of the result of Section 4.2 – if condition (76) is satisfied, i.e., if the external mass density $\rho_c$ is a non-increasing function of $r$ (which should be the case in the applications).

For anisotropic systems satisfying (76) and (77), our Theorem 2 shows that any non-trivial aspherical perturbation generated by a preserving $\delta g_1$ has positive energy – the same result holding true for any preserving $g_1$ because of Theorem 1. This completely new result is unfortunately more difficult to interpret. Strictly speaking, the general stability criteria established thus far do not apply here, as the preserving character of a perturbation is conserved in time only for spherically symmetric perturbations ($\{f_0, g_1\} = 0$). However, we feel that our result has some relevance for the stability problem. Numerical simulations giving some support to this idea are reported in Paper II.

## 5.5 Comparison with previous works

For isotropic systems, the stability result above has already been derived by the thermodynamic method and the energy principle (Antonov 1961; Lynden-Bell & Sannit 1969; Binney & Tremaine 1987, where it is referred to as the Antonov-Lebowitz theorem) in the case where there is no external potential (i.e., $\Phi_c \equiv 0$). The new proof we have presented here, however, is more simple and transparent, and it is self-contained (e.g., it does not use results on the stability of gaseous stars). Moreover, it does not appeal to some unproven assumptions on the completeness of the modes of some operator (assumptions which may certainly be proven, but with much effort). Our work thus answers an important question addressed by Kandrup & Sygnet (1985) concerning the possibility of finding a stability proof that does not use some peculiar completeness theorem, which certainly justified its presentation in detail.

Concerning the aspherical stability of an anisotropic system, there seems to be only one published result based on the consideration of an energy functional. This is by Gillon et al. (1976), and it states that any spherically symmetric steady state satisfying $\Phi_c \equiv 0$ and conditions (18) and (77) is stable with respect to any aspherical perturbations. Some doubts have been repeatedly cast in the literature on the validity of this result (e.g., Binney & Tremaine 1987), but no explicit rebuttal (or approval!) has been yet presented. Owing to its possible importance, we have thus conducted a detailed analysis of Gillon et al.'s complicated proof, in which a spherical system is approximated by a multiple water-bag model. Unfortunately, we have found two serious flaws and our conclusion does confirm the fear of the sceptics.

(a) There is an inconsistency in the calculations, as the conservation law for the basic quantity $A$ (their equation 49) is obtained by using both the non-linear and the linearized VP systems, without distinguishing carefully the order of the various terms. When the derivation is made correctly, it does not seem to lead to an interesting result.

(b) If we forget about this inconsistency and admit for a while the validity of their equation (49), then we must note that the ‘conserved’ quantity $A$ (their equation 52) takes a positive value $A_{eq} > 0$ at equilibrium (just make $a_{\pm k} = 0$ in their equation 52, where $a_{\pm k}$ are the functions describing an arbitrary perturbation). Therefore, to be a Lyapunov function, $A$ should satisfy the condition $A \geq A_{eq}$ for any perturbation rather than the stated one, $A \geq 0$. It is, however, easy to check that the former condition is certainly not fulfilled, as $A$ can be made larger or smaller than $A_{eq}$ by choosing adequately $a_{\pm k}$. Then, even under our optimistic assumption, no conclusion could be actually drawn from the sign of $A$.

# 6 EXTENSION TO OTHER SYSTEMS

## 6.1 Rotating spherical systems

Consider a spherically symmetric steady state characterized by the function $f_s(E, L^2)$ and the potential $\Psi_0$, and introduce an arbitrary function $f_a(E, L^2, L_z)$ such that

$$f_0(E, L^2, L_z) := f_s(E, L^2) + f_a(E, L^2, L_z) \geq 0, \quad (81)$$
and

$$f_a(E, L^2, -L_z) = -f_a(E, L^2, L_z). \quad (82)$$

As noted long ago by Lynden-Bell (1960), $f_0$ is also a steady-state distribution function associated with the potential $\Psi_0$ (the mass density corresponding to $f_0$ is the same as the one corresponding to $f_s$), and it describes a system which is rotating about the z-axis. We now discuss whether the stability results derived above can be extended to these more general distribution functions.

Proceeding as in Section 3.3, we obtain, after some algebra,

$$ \begin{aligned} H^{(2)}[\bar{g}_1] = & -\frac{1}{2} \int \{f_s, \bar{g}_1\} \{E, \bar{g}_1^*\} d\xi - \frac{1}{2} \int \{f_0, \delta g_1\} \\ & \times \{E, \delta g_1^*\} d\xi - \int \{f_a, \bar{g}_1\} \{E, \delta g_1^*\} d\xi \\ & - \frac{1}{8\pi G} \int |\nabla_x \bar{\Phi}_1|^2 dx - \frac{1}{8\pi G} \int |\nabla_x \delta\Phi_1|^2 dx, \end{aligned} \quad (83) $$

with

$$\bar{\Phi}_1 = K[ - \{f_s, \bar{g}_1\} - \{f_a, \delta g_1\} ], \quad (84)$$

$$\delta\Phi_1 = K[ - \{f_s, \delta g_1\} - \{f_a, \delta g_1\} - \{f_a, \overline{\delta g_1}\} ]. \quad (85)$$

Because of the presence of the non-spherically symmetric term $f_a$ in $f_0$, it is then no longer possible to separate $H^{(2)}[\bar{g}_1]$ into two terms depending respectively on $\bar{g}_1$ and $\delta g_1$, and the analysis becomes rather intricate.

However, there are two simple – but new – results which can be extracted from the expressions above when $\partial f_s / \partial E < 0$. Consider first a purely spherically symmetric perturbation generated by $\bar{g}_1$. Then

$$ \begin{aligned} H^{(2)}[\bar{g}_1] = H^{(2)}[\bar{g}_1] = & -\frac{1}{2} \int \{f_s, \bar{g}_1\} \{E, \bar{g}_1^*\} d\xi \\ & -\frac{1}{8\pi G} \int |\nabla_x \bar{\Phi}_1|^2 dx, \end{aligned} \quad (86) $$

$$\bar{\Phi}_1 = K[ - \{f_s, \bar{g}_1\} ], \quad (87)$$

that is, we recover the expressions valid for the spherical state $f_s$. It thus results at once from Section 4 that $H^{(2)}[\bar{g}_1] > 0$ for any non-trivial $\bar{g}_1$, and our rotating state turns out to be stable with respect to all the spherically symmetric perturbations (their character being preserved in time).

Consider now a system with $f_s = f_s(E)$ which is submitted to some axisymmetric perturbation ($\partial \delta g_1 / \partial \phi = 0$) satisfying $\delta g_1(-v_\phi) = \delta g_1(v_\phi)$ (a property also conserved in time): stability is then obtained in that case too. As is easily checked, the cross term in equation (83) [containing both $\bar{g}_1$ and $\delta g_1$] and the perturbed potential associated with $f_a$ vanish, and the results of Sections 4 and 5 imply that $H^{(2)}[\bar{g}_1] > 0$.

## 6.2 Systems with a non-decreasing distribution function

The stability of a spherical system when condition (18) is not satisfied has been considered numerically by Henon (1973). From an analytical point of view it appears to be quite a difficult problem, and here we shall content ourselves with making two points.

(a) Consider the case where $f_{0E} > 0$ in the region $\Omega_0$ where $f_0 > 0$. This of course implies that $f_0$ is discontinuous on the boundary $\partial\Omega_0 = \{\xi | E(\xi) = E_0 \le 0\}$, and there is a violation of our smoothness assumption. Consider a spherical perturbation generated by $\bar{g}_1$, and assume that $\{E, \bar{g}_1\}$ vanishes on $\partial\Omega_0$. It is easy to check that $H^{(2)}[\bar{g}_1]$ is still given by equation (39) [the boundary term related to the presence of the delta function $\delta(E - E_0)$ for the derivative $f_{0E}$ disappearing], and that $H^{(2)}[\bar{g}_1] < 0$ for any non-trivial $\bar{g}_1$. Then the system appears to be stable with respect to these particular perturbations (Section 3.2a) – a conclusion also obtained by Kandrup & Sygnet (1985) by means of their energy principle.

(b) Consider the case where the system is isotropic [$f_0 = f_0(E)$] and $f_{0E}$ may have a changing sign in $\Omega_0$, and let us assume that the system admits a neutral aspherical mode generated by $\delta g_1$. Then we have $\delta f_1 = 0$ and $H^{(2)}[\bar{g}_1] = 0$ (Perez 1995), which implies at once, by using the arguments at the end of Section 5.1, that $\delta f_1 = 0$. We have thus proved the following new result: whatever the sign of $f_{0E}$ is, an isotropic spherically symmetric system cannot admit an aspherical neutral mode.

## 7 CONCLUSION

In this paper we have considered a self-gravitating collisionless system, in which the particles are also submitted to an external potential $\Phi_c$. We have first presented some new considerations on the relations between the spherical symmetry of a steady state and the form of its distribution function $f_0$. In particular, we have argued that, under quite general assumptions, a system with an $f_0$ depending only on $E$ and $L^2$ is necessarily spherically symmetric.

We have thus reinvestigated from an analytical point of view the problem of the linear stability of this system in the consistent framework of the symplectic method. By effecting a systematic study of the sign of the energy of an arbitrary symplectic perturbation, we have reached in particular the following conclusions, valid when $f_{0E} < 0$.

(a) Any steady state is stable with respect to all the spherically symmetric perturbations.

(b) An isotropic system [$f_0 = f_0(E)$] is also stable with respect to all the aspherical perturbations, at least if the external mass density $\rho_c$ is a non-increasing function of the radial coordinate $r$. Of course, this result and the one in (a) are well known in the case where $\Phi_c \equiv 0$; but even if we restrict our attention to this particular case, our proofs are much simpler and complete than all the previous ones, bringing about some new important elements which were called for by many authors. In particular, our proof of the second quoted result shows more clearly the origin of the stability, and it is the first which does not appeal to some unproved assumption on the completeness of a set of modes.
(c) For an anisotropic system with $f_{\text{OL}^2} \le 0$ and $\rho'_c \le 0$, any non-trivial 'preserving perturbation' – i.e., one having a generator commuting with $L^2$, which of course does not mean that it is spherically symmetric – has a positive energy. The exact meaning of this result is still not completely clear, and we have tried to elucidate it by effecting numerical investigations. The results obtained thus far (in particular for the case where the distribution function is of the Ossipkov–Meritt type – see, e.g., Binney & Tremaine 1987 – and there is no external potential) are reported in the companion paper (Paper II). More general situations, currently under intensive study, will be analysed in forthcoming papers.

Finally,

(a) we have analysed in detail a challenged result of Gil- lon et al. (1976), according to which stability with respect to any perturbation holds if $f_{\text{OL}^2} \le 0$, and found the proof to contain some flaws, which we have clearly exposed;

(b) we have shown that a spherical rotating system is stable with respect to all the spherically symmetric perturba- tions and [when $f_s=f_s(E)$] to all the axisymmetric even ones;

(c) we have proved that an isotropic system with a deriva- tive $f_{\text{DE}}$, which may change its sign, cannot admit a non- trivial neutral mode (i.e., one differing from a mere global translation of the system).

REFERENCES

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

Antonov V. A., 1961, SvA, 4, 859

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

Bartholomew P., 1971, MNRAS, 151, 333

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

Dejonghe H., 1986, Phys. Rep., 133, 217

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

Gidas B., Ni W. M., Nirenberg L., 1981, Math. Analysis and Appli- cations, Part A: Advances in Math. Supplementary Studies, 7A, 369

Gillon D., Doremus J. P., Baumann G., 1976, A&A, 48, 467

Henon M., 1973, A&A, 24, 229

Ipser J. R., 1974, ApJ, 193, 463

Ipser J. R., Horwitz P., 1979, ApJ, 232, 863

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

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

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

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

Landau L. D., Lifshitz E., 1966, Cours de Physique Théorique – Tome 1: Mécanique. Mir, Moscou

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

Lynden-Bell D., 1960, MNRAS, 120, 204

Lynden-Bell D., 1966, in Contopoulos G., ed., The Theory of Orbits in the Solar System and in Stellar Systems. Academic Press, London and New York

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

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

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

Perez J., Alimi J. M., Aly J. J., Scholl H., 1996, MNRAS, 280, 700 (Paper II, this issue)

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