Dynamical stability of infinite homogeneous self-gravitating systems: application of the Nyquist method

We consider a self-gravitating gas described by the Euler-Poisson system

$${\partial\rho\over\partial t}+\nabla\cdot(\rho{\bf u})=0,$$

(1)

$${\partial{\bf u}\over\partial t}+({\bf u}\cdot\nabla){\bf u}=-{1\over\rho}\nabla p-\nabla\Phi,$$

(2)

$$\Delta\Phi=4\pi G\rho.$$

(3)

These are the usual equations used to study the dynamical stability of a gaseous medium like a molecular cloud or a star. We assume that the gas is a perfect gas in local thermodynamic equilibrium (L.T.E). Therefore, at each point, the distribution function of the particles (atoms or molecules) is of the form

$$f({\bf r},{\bf v},t)=\biggl{[}{1\over 2\pi T({\bf r},t)}\biggr{]}^{3/2}\rho({\bf r},t)e^{-{[{\bf v}-{\bf u}({\bf r},t)]^{2}\over 2T({\bf r},t)}}.$$

(4)

Note that the kinetic temperature $T({\bf r},t)=\frac{1}{3}\langle w^{2}\rangle$ where ${\bf w}={\bf v}-{\bf u}({\bf r},t)$ is the relative velocity⁸⁸8Throughout this paper, by an abuse of language, we shall define the kinetic temperature by $T({\bf r},t)\equiv\frac{1}{d}\langle({\bf v}-{\bf u}({\bf r},t))^{2}\rangle$ where $d$ is the dimension of space. Therefore, the kinetic temperature is a measure of the local velocity dispersion of the particles in one dimension. It is related to the real kinetic temperature by $T({\bf r},t)=k_{B}T_{real}({\bf r},t)/m$. We do that in order to extend this definition to the case of collisionless stellar systems described by the Vlasov equation where the individual mass of the stars does not appear.. Introducing the pressure $p={1\over 3}\int fw^{2}d{\bf v}$, we find that the local equation of state is

$$p({\bf r},t)=\rho({\bf r},t)T({\bf r},t).$$

(5)

The system of equations (1)-(3) is not closed since it must be supplemented by an equation for the energy or, equivalently, for the temperature $T({\bf r},t)$. Here, we shall restrict ourselves to the case of a barotropic gas for which the pressure depends only on the density

$$p({\bf r},t)=p[\rho({\bf r},t)].$$

(6)

The local velocity of sound is

$$c_{s}^{2}({\bf r},t)=p^{\prime}(\rho({\bf r},t)).$$

(7)

A steady state of the Euler-Poisson system satisfies ${\bf u}={\bf 0}$ and the condition of hydrostatic equilibrium

$$\nabla p+\rho\nabla\Phi={\bf 0}.$$

(8)

Combining the condition of hydrostatic equilibrium (8) and the equation of state $p=p(\rho)$, we find that the density is a function $\rho=\rho(\Phi)$ of the gravitational potential given by

$\displaystyle\int^{\rho}{p^{\prime}(\rho^{\prime})\over\rho^{\prime}}d\rho^{\prime}=-\Phi.$

(9)

We note that $p^{\prime}(\rho)=-\rho/\rho^{\prime}(\Phi)$. Since $p^{\prime}(\rho)=c_{s}^{2}\geq 0$, we find that $\rho^{\prime}(\Phi)\leq 0$. The density is a monotonically decreasing function of the gravitational potential. The total energy of the star is

	$\displaystyle{\cal W}[\rho,{\bf u}]=\int\rho\int_{0}^{\rho}{p(\rho^{\prime})\over\rho^{{}^{\prime}2}}d\rho^{\prime}d{\bf r}+{1\over 2}\int\rho\Phi d{\bf r}+\int\rho{{\bf u}^{2}\over 2}d{\bf r},$
			(10)

including the internal energy, the potential energy and the kinetic energy of the mean motion. The mass $M[\rho]$ and the energy ${\cal W}[\rho,{\bf u}]$ are conserved by the barotropic Euler-Poisson system: $\dot{M}=\dot{\cal W}=0$. Therefore, a minimum of the energy functional at fixed mass determines a steady state of the barotropic Euler-Poisson system that is nonlinearly dynamically stable. Here, we have stability in the sense of Lyapunov. This means that the size of the perturbation is bounded by the size of the initial perturbation for all times. We are led therefore to considering the minimization problem

$$\min_{\rho,{\bf u}}\ \left\{{\cal W}[\rho]\quad|\quad M[\rho]=M\right\}.$$

(11)

The critical points of energy at fixed mass are determined by the Euler-Lagrange equation $\delta{\cal W}-\mu\delta M=0$ where $\mu$ is a Lagrange multiplier. This yields the condition of hydrostatic equilibrium (8). Then, a critical point of energy at fixed mass is a (local) energy minimum iff

$$\delta^{2}{\cal W}=\int\frac{p^{\prime}(\rho)}{2\rho}(\delta\rho)^{2}\,d{\bf r}+\frac{1}{2}\int\delta\rho\delta\Phi\,d{\bf r}\geq 0$$

(12)

for all perturbations $\delta\rho$ that conserve mass: $\delta M=0$.

We consider the linear dynamical stability of an infinite homogeneous gaseous medium described by the hydrodynamic equations (1), (2), (3) and (6). This is the classical Jeans (1902,1929) problem. Linearizing Eqs. (1)-(3) around a solution $\rho=\Phi={\rm const.}$ and ${\bf u}={\bf 0}$, making the Jeans swindle and decomposing the perturbations in normal modes of the form $e^{i({\bf k}\cdot{\bf r}-\omega t)}$, we obtain the classical dispersion relation (Binney & Tremaine 1987):

$\displaystyle\omega^{2}=c_{s}^{2}k^{2}-4\pi G\rho,$

(13)

where $c_{s}^{2}=dp/d\rho$ is the velocity of sound. Since $\omega^{2}$ is real, the complex pulsation $\omega$ is either real or purely imaginary. The condition $\omega=0$ determines a critical wavenumber

$\displaystyle k_{c}^{2}=\frac{4\pi G\rho}{c_{s}^{2}},$

(14)

called the Jeans wavenumber for a barotropic gas. The dispersion relation can be rewritten

$\displaystyle\frac{\omega^{2}}{4\pi G\rho}=\frac{k^{2}}{k_{c}^{2}}-1.$

(15)

For $k>k_{c}$, we have $\omega=\omega_{r}=\pm\sqrt{\omega^{2}}$ so that the perturbation oscillates like $e^{-i\omega_{r}t}$ with a pulsation

$\displaystyle\omega_{r}=\pm\sqrt{c_{s}^{2}k^{2}-4\pi G\rho},$

(16)

without attenuation. This corresponds to gravity-modified sound waves. In that case the system is stable. For $k<k_{c}$, we have $\omega=\omega_{i}=\pm i\sqrt{-\omega^{2}}$ so that the perturbation grows like $e^{\omega_{i}t}$ with a growth rate

$\displaystyle\omega_{i}=\sqrt{4\pi G\rho-c_{s}^{2}k^{2}},$

(17)

without oscillation (the second mode is attenuated exponentially rapidly so that the growing mode dominates). In that case, the system is unstable. This is the so-called Jeans instability. The growth rate is maximum for $k=0$ and its value is $\omega_{i}(k=0)=\sqrt{4\pi G\rho}$.

In conclusion, a perturbation with wavenumber $k$ is stable if

$\displaystyle k>k_{c}=\frac{\sqrt{4\pi G\rho}}{c_{s}},$

(18)

and unstable otherwise. These results are valid for an arbitrary barotropic equation of state $p=p(\rho)$. We now specialize on particular cases, namely isothermal and polytropic equations of state.

We first consider an isothermal equation of state

$$p({\bf r},t)=\rho({\bf r},t)T,$$

(19)

where the temperature $T({\bf r},t)=T$ is a uniform. At hydrostatic equilibrium, according to Eq. (9), the relation between the density and the gravitational potential is given by the Boltzmann law

$$\rho=A^{\prime}e^{-{\Phi\over T}},$$

(20)

where $A^{\prime}$ is a constant. The energy functional (10) reads

	$\displaystyle{\cal W}[\rho,{\bf u}]=k_{B}T\int{\rho}\ln{\rho}d{\bf r}+{1\over 2}\int\rho\Phi d{\bf r}+\int\rho{{\bf u}^{2}\over 2}d{\bf r},$
			(21)

and it can be viewed as a Boltzmann free energy $F=E-TS$ where $E$ is the macroscopic energy and $S$ the Boltzmann entropy. The velocity of sound (7) is uniform:

$$c_{s}^{2}=T.$$

(22)

Finally, for an infinite homogeneous isothermal gas, the Jeans instability criterion takes the classical form

$\displaystyle k^{2}<k^{2}_{c}=\frac{4\pi G\rho}{T}.$

(23)

The equation of state of a polytropic gas⁹⁹9The polytropic equation of state corresponds to an adiabatic transformation for which the specific entropy $s({\bf r},t)=s$ is uniform. In that case, $\gamma=c_{p}/c_{v}$ is the ratio of specific heats at constant pressure and constant volume. For a monoatomic gas $\gamma=5/3$. This law describes convective regions of stellar interior. An approximate polytropic equation of state with index $\gamma\simeq 3.25$ also holds in the radiative region of a star. Finally, classical white dwarf stars are described by a polytropic equation of state with index $\gamma=5/3$ ($n=3/2$). This equation of state is valid for a completely degenerate gas of electrons described by the Fermi-Dirac statistics at $T=0$ (Chandrasekhar 1942). is

$$p({\bf r},t)=K\rho({\bf r},t)^{\gamma},$$

(24)

where $K$ is a constant. The polytropic index $n$ is defined by $\gamma=1+1/n$. For $n\rightarrow+\infty$, we recover the isothermal case with $\gamma=1$. At hydrostatic equilibrium, according to Eq. (9), the relation between the density and the gravitational potential can be written

$\displaystyle\rho=\biggl{[}\lambda-{\gamma-1\over K\gamma}\Phi\biggr{]}^{1\over\gamma-1},$

(25)

where $\lambda$ is a constant. The energy functional (10) can be put in the form

	$\displaystyle{\cal W}[\rho,{\bf u}]={K\over\gamma-1}\int(\rho^{\gamma}-\rho)d{\bf r}+{1\over 2}\int\rho\Phi d{\bf r}+\int\rho{{\bf u}^{2}\over 2}d{\bf r}.$
			(26)

In the limit $\gamma\rightarrow 1$, we recover the results (20) and (21) obtained for isothermal distributions. For a polytropic gas, the local temperature $T({\bf r},t)$ given by Eq. (5) reads

$$T({\bf r},t)=K\rho({\bf r},t)^{\gamma-1}.$$

(27)

We note that the temperature is position dependent (while the specific entropy $s$ is uniform) and related to the density by a power law (this is the local version of the usual isentropic law $TV^{\gamma-1}={\rm const.}$ in thermodynamics). The polytropic index $n$ is related to the gradients of temperature and density according to

$${\nabla\ln T}={1\over n}\nabla\ln\rho.$$

(28)

Combining Eqs. (25) and (27), we note that, for a polytropic gas at equilibrium, the relation between the kinetic temperature and the gravitational potential is linear so that

$$\nabla T=-\frac{\gamma-1}{\gamma}\nabla\Phi.$$

(29)

The coefficient of proportionality is related to the adiabatic index $\gamma$. The velocity of sound (7) can be expressed as

$$c_{s}({\bf r},t)^{2}=\gamma T({\bf r},t),$$

(30)

and it is usually position dependent. However, when we consider an infinite and homogeneous distribution (making the Jeans swindle), the velocity of sound and the temperature are uniform. In that case, we can speak of “the temperature $T$ of the polytropic gas” and write the Jeans instability criterion as

$\displaystyle k^{2}<k^{2}_{c}=\frac{4\pi G\rho}{\gamma T}.$

(31)

We thus find that the critical Jeans wavenumber $k_{c}^{(poly)}$ of a polytropic gas is related to the critical Jeans wavenumber $k_{c}^{(iso)}$ of an isothermal gas with the same temperature $T$ by

$\displaystyle k_{c}^{(poly)}=\frac{1}{\sqrt{\gamma}}k_{c}^{(iso)},$

(32)

where the proportionality factor is the polytropic index $\gamma$. The positivity of $c_{s}^{2}=dp/d\rho$ implies that $\gamma\geq 0$, i.e. $n\geq 0$ or $n\leq-1$ (for $-1<n<0$, the gas is always unstable, even in the absence of gravity, since $\omega^{2}<0$). For $\gamma=0$ ($n=-1$), the gas is unstable to all wavelengths since $k_{c}^{(poly)}=+\infty$. For $0<\gamma<1$ ($n<-1$), $k_{c}^{(poly)}>k_{c}^{(iso)}$. For $\gamma=1$ ($n=\infty$), $k_{c}^{(poly)}=k_{c}^{(iso)}$ (isothermal case). For $\gamma>1$ ($n>0$), $k_{c}^{(poly)}<k_{c}^{(iso)}$. For $\gamma=+\infty$ ($n=0$), the gas is stable to all wavelengths since $k_{c}^{(poly)}=0$ (solid medium). For an isentropic gas for which $\gamma=c_{p}/c_{v}$, the Mayer relation $c_{p}-c_{v}=k_{B}$ implies that $\gamma>1$ ($n>0$). Therefore, for an isentropic gas, the gravitational instability is always delayed (i.e., it occurs for larger wavelengths) with respect to an isothermal gas with the same temperature.

We consider a self-gravitating system described by the Vlasov-Poisson system

$${\partial f\over\partial t}+{\bf v}\cdot{\partial f\over\partial{\bf r}}+{\bf F}\cdot{\partial f\over\partial{\bf v}}=0,$$

(33)

$$\Delta\Phi=4\pi G\int fd{\bf v},$$

(34)

where ${\bf F}=-\nabla\Phi$ is the self-consistent gravitational field produced by the particles. The Vlasov description assumes that the evolution of the system is encounterless. This is a very good approximation for most astrophysical systems like galaxies and dark matter because the relaxation time due to close encounters is in general larger than the age of the universe by several orders of magnitude (Binney & Tremaine 1987). The Vlasov equation admits an infinite number of stationary solutions given by the Jeans theorem (Jeans 1915). For example, distribution functions of the form $f=f(\epsilon)$ which depend only on the individual energy $\epsilon={v^{2}}/{2}+\Phi({\bf r})$ of the stars are particular steady states of the Vlasov equation. They describe spherical stellar systems.

The Vlasov equation conserves the mass $M[f]=\int f\,d{\bf r}d{\bf v}$, the energy $E[f]=\frac{1}{2}\int fv^{2}\,d{\bf r}d{\bf v}+\frac{1}{2}\int\rho\Phi\,d{\bf r}$ and the Casimirs $I_{h}=\int h(f)\,d{\bf r}d{\bf v}$ for any continuous function $h$. Let us introduce the functionals

$$S=-\int C(f)\,d{\bf r}d{\bf v},$$

(35)

where $C$ is a convex function ($C^{\prime\prime}>0$). These particular Casimirs will be called “pseudo entropies”¹⁰¹⁰10The functionals $H=-\int C(\bar{f})\,d{\bf r}d{\bf v}$ defined in terms of the coarse-grained distribution $\bar{f}$ are called generalized $H$-functions (Tremaine et al. 1986)..

The two constraints problem: since the Vlasov equation conserves $M$, $E$ and $S$, the maximization problem

$$\max_{f}\ \left\{S[f]\quad|\quad E[f]=E,\quad M[f]=M\right\},$$

(36)

determines a steady state of the Vlasov-Poisson system that is nonlinearly dynamically stable. The critical points of pseudo entropy at fixed mass and energy are determined by the variational principle $\delta S-\beta\delta E-\alpha\delta M=0$, where $\beta$ (pseudo inverse temperature) and $\alpha$ (pseudo chemical potential) are Lagrange multipliers. This yields $C^{\prime}(f)=-\beta\epsilon-\alpha$. Since $C$ is convex, this relation can be reversed to give $f=F(\beta\epsilon+\alpha)$ where $F=(C^{\prime})^{-1}(-x)$. We have $f^{\prime}(\epsilon)=-\beta/C^{\prime\prime}(f)$. Therefore, if $\beta>0$, the critical points of pseudo entropy at fixed mass and energy determine distributions of the form $f=f(\epsilon)$ with $f^{\prime}(\epsilon)<0$: the distribution function is a monotonically decreasing function of the energy¹¹¹¹11More generally, the solutions of (36) are of the form $f=f(\epsilon)$ where $f$ is monotonic, decreasing at positive temperatures and increasing at negative temperatures. For realistic stellar systems, the DF should decrease close to the escape energy $\epsilon=0$. Therefore, for the class of distributions considered, $f$ must decrease everywhere implying $\beta>0$.. Then, a critical point of pseudo entropy at fixed mass and energy is a (local) maximum iff

$$\delta^{2}{S}=-\frac{1}{2}\int C^{\prime\prime}(f)(\delta f)^{2}\,d{\bf r}d{\bf v}-\frac{1}{2}\beta\int\delta\rho\delta\Phi\,d{\bf r}<0,$$

(37)

for all perturbations $\delta f$ that conserve mass and energy at first order: $\delta M=\delta E=0$.

The one constraint problem: the minimization problem

$$\min_{f}\ \left\{F[f]=E[f]-TS[f]\quad|\quad M[f]=M\right\},$$

(38)

where $T=1/\beta>0$ is prescribed, determines a steady state of the Vlasov-Poisson system that is nonlinearly dynamically stable. The functional $F[f]$ is a pseudo free energy. The critical points of pseudo free energy at fixed mass are determined by the variational principle $\delta F+\alpha T\delta M=0$, where $\alpha$ (pseudo chemical potential) is a Lagrange multiplier. This yields $C^{\prime}(f)=-\beta\epsilon-\alpha$. Then, a critical point of pseudo free energy at fixed mass is a (local) minimum iff

$$\delta^{2}{F}=\frac{1}{2}T\int C^{\prime\prime}(f)(\delta f)^{2}\,d{\bf r}d{\bf v}+\frac{1}{2}\int\delta\rho\delta\Phi\,d{\bf r}>0,$$

(39)

for all perturbations $\delta f$ that conserve mass: $\delta M=0$.

The no constraint problem: the maximization problem

$$\max_{f}\ \left\{G[f]=S[f]-\beta E[f]-\alpha M[f]\right\},$$

(40)

where $\beta$ and $\alpha$ are prescribed, determines a steady state of the Vlasov-Poisson system that is nonlinearly dynamically stable. The functional $G[f]$ is a pseudo grand potential. The critical points of pseudo grand potential $G$, satisfying $\delta G=0$, are given by $C^{\prime}(f)=-\beta\epsilon-\alpha$. Then, a critical point of pseudo grand potential is a (local) maximum iff

$$\delta^{2}{G}=-\frac{1}{2}\int C^{\prime\prime}(f)(\delta f)^{2}\,d{\bf r}d{\bf v}-\frac{1}{2}\beta\int\delta\rho\delta\Phi\,d{\bf r}<0,$$

(41)

for all perturbations $\delta f$.

The optimization problems (36), (38) and (40) have the same critical points (canceling the first order variations). Furthermore, a solution of (40) is always a solution of the more constrained dual problem (38). Indeed, if inequality (41) is true for all perturbations, it is true a fortiori for all perturbations that conserve mass. Similarly, a solution of (38) is always a solution of the more constrained dual problem (36). Indeed, if inequality (39) is true for all perturbations that conserve mass, it is true a fortiori for all perturbations that conserve mass and energy. However, the reciprocal is wrong. A solution of (38) is not necessarily a solution of (40), and a solution of (36) is not necessarily a solution of (38). This is similar to the notion of ensemble inequivalence in thermodynamics (Ellis et al. 2000, Bouchet & Barré 2005, Chavanis 2006a). Indeed, the two constraints problem (36) is similar to a condition of microcanonical stability, the one constraint problem (38) is similar to a condition of canonical stability, and the no constraint problem (40) is similar to a condition of grand canonical stability. The implication $(\ref{erica1})\Rightarrow(\ref{vh4})\Rightarrow(\ref{vh3})$ is similar to the fact that grand canonical stability implies canonical stability which itself implies microcanonical stability (but not the converse) in thermodynamics. Therefore, (40) provides just a sufficient condition of nonlinear dynamical stability that is less refined than (38), and (38) provides just a sufficient condition of nonlinear dynamical stability that is less refined than (36).

The most refined problem: the minimization problem

$$\min_{f}\ \left\{E[f]\quad|{\rm\ all\ the\ Casimirs\ }I_{h}\right\},$$

(42)

determines a steady state of the Vlasov-Poisson system that is nonlinearly dynamically stable. A distribution function is a critical point of energy for symplectic perturbations (i.e. perturbations that conserve all the Casimirs) iff $f({\bf r},{\bf v})$ is a steady state of the Vlasov equation (Bartholomew 1971, Kandrup 1991). Furthermore, if we consider spherical stellar systems for which $f=f(\epsilon)$, it can be shown (Bartholomew 1971, Kandrup 1991) that a critical point of energy for symplectic perturbations is a (local) minimum iff

$$\delta^{2}{E}=-\frac{1}{2}\int\frac{(\delta f)^{2}}{f^{\prime}(\epsilon)}\,d{\bf r}d{\bf v}+\frac{1}{2}\int\delta\rho\delta\Phi\,d{\bf r}>0,$$

(43)

for all perturbations $\delta f$ that conserve energy and all the Casimirs at first order: $\delta E=\delta I_{h}=0$. Such symplectic (physically accessible) perturbations are of the form $\delta f=D\delta g$ where $D={\bf v}\cdot\nabla_{\bf r}-\nabla\Phi\cdot\nabla_{\bf v}$ is the advective operator in phase space and $\delta g({\bf r},{\bf v})$ is any perturbation. Inequality (43) corresponds to the Antonov (1960) criterion of dynamical stability that was obtained by investigating the linear dynamical stability of a steady state of the Vlasov equation. Since this is the most constrained criterion, this is the most refined one. We note that inequality (43) is equivalent to inequalities (37), (39) and (41) if we use the identity $f^{\prime}(\epsilon)=-\beta/C^{\prime\prime}(f)$ derived above. However, the classes of perturbations to consider are different: in (41), we must consider all perturbations, in (39) we must consider perturbations that conserve mass, in (37) we must consider perturbations that conserve mass and energy, and in (43) we must consider perturbations that conserve energy and all the Casimirs. Therefore, we have the implications $(\ref{erica1})\Rightarrow(\ref{vh4})\Rightarrow(\ref{vh3})\Rightarrow(\ref{erica3})$. The connection between (42) and the preceding optimization problems can be understood as follows. The minimization problem (42), conserving all the Casimirs, is clearly more refined than the minimization problem

$$\min_{f}\ \left\{E[f]\quad|\quad M[f]=M,\quad S[f]=S\right\},$$

(44)

where only the mass and one Casimir of the form (35) is conserved. Now, it is easy to show that the minimization problem (44) is equivalent to the maximization problem (36). Hence, we have the chain of relations

$$(\ref{erica1})\Rightarrow(\ref{vh4})\Rightarrow(\ref{vh3})\Leftrightarrow(\ref{erica5})\Rightarrow(\ref{erica3}).$$

(45)

To summarize, the minimization problem (42) is the most refined stability criterion because it tells that, in order to settle the dynamical stability of a stellar system, we just need considering symplectic (i.e. dynamically accessible) perturbations. Of course, if inequality (43) is satisfied by a larger class of perturbations, as implied by problems (40), (38), (36) and (44), the system will be stable a fortiori. Therefore, we have the implications (45). Problems (40), (38), (36) and (44) provide sufficient (but not necessary) conditions of dynamical stability. A steady state can be stable according to (38) while it does not satisfy (40), or it can be stable according to (36) and (44) while it does not satisfy (38), or it can be stable according to (42) while it does not satisfy (36) and (44). Therefore, the criterion (42) is more refined than (44) or (36), which is itself more refined than (38), which is itself more refined than (40). Let us give an astrophysical illustration of all that (Chavanis 2006a). Using (38), we can show that stellar polytropes with $3/2\leq n\leq 3$ are dynamically stable because they are minima of $F$ at fixed mass $M$. By contrast, stellar polytropes with $3<n<5$ are not minima of $F$ at fixed mass. However, using (36), we can show that stellar polytropes with $3/2<n<5$ are dynamically stable because they are maxima of $S$ at fixed mass and energy. However, not all stellar systems of the form $f=f(\epsilon)$ are maxima of $S$ at fixed mass and energy. But, using (42), we can show that all spherical galaxies $f=f(\epsilon)$ with $f^{\prime}(\epsilon)<0$ are dynamically stable. This last statement has been shown for linear stability by Kandrup (1991). A rigorous mathematical proof has been given recently by Lemou et al. (2009) for nonlinear stability.

Remark: similar results are obtained in 2D fluid mechanics based on the Euler-Poisson system (see Chavanis 2009). Criterion (42) is equivalent to the so-called Kelvin-Arnol’d energy principle, criterion (40) is equivalent to the standard Casimir-energy method (see Holm et al. 1985) introduced by Arnol’d (1966) and criterion (36) is equivalent to the refined stability criterion given by Ellis et al. (2002).

To any stellar system with $f=f(\epsilon)$ and $f^{\prime}(\epsilon)<0$, we can associate a corresponding barotropic star with the same equilibrium density distribution (Lynden-Bell & Sanitt 1969). Indeed, defining the density and the pressure by $\rho=\int fd{\bf v}$ and $p={1\over 3}\int fv^{2}d{\bf v}$, we get $\rho=\rho(\Phi)$ and $p=p(\Phi)$. Eliminating the potential $\Phi$ between these two expressions, we find that $p=p(\rho)$. Furthermore, taking the gradient of the pressure and using the chain of identities

	$\displaystyle\nabla p={1\over 3}\int{\partial f\over\partial{\bf r}}v^{2}d{\bf v}={1\over 3}\nabla\Phi\int f^{\prime}(\epsilon)v^{2}d{\bf v}$
	$\displaystyle={1\over 3}\nabla\Phi\int{\partial f\over\partial{\bf v}}\cdot{\bf v}\,d{\bf v}=-\nabla\Phi\int f\,d{\bf v}=-\rho\nabla\Phi,$		(46)

we obtain the condition of hydrostatic equilibrium (8). Finally, we define the kinetic temperature of an isotropic stellar system by the relation

$\displaystyle\frac{3}{2}T({\bf r})=\frac{1}{2}\langle v^{2}\rangle.$

(47)

Therefore, the quantity

$\displaystyle T({\bf r})=\frac{1}{3}\langle v^{2}\rangle={{1\over 3}\int fv^{2}d{\bf v}\over\int fd{\bf v}}={p({\bf r})\over\rho({\bf r})}.$

(48)

measures the velocity dispersion (in one direction) of an isotropic stellar system. In general, the kinetic temperature is position dependent.

Finally, we can show that the variational principles $(\ref{vh4})$ and $(\ref{ep10})$ are equivalent, i.e. a stellar system is a minimum of pseudo free energy at fixed mass iff the corresponding barotropic gas is a minimum of energy at fixed mass (see Appendix A). This leads to a nonlinear version of the Antonov (1960b) first law: “a stellar system with $f=f(\epsilon)$ and $f^{\prime}(\epsilon)<0$ is nonlinearly dynamically stable with respect to the Vlasov-Poisson system if the corresponding barotropic gas is nonlinearly dynamically stable with respect to the Euler-Poisson system”. However, the reciprocal is wrong because, as we have already indicated, (38) provides just a sufficient condition of nonlinear dynamical stability with respect to the Vlasov equation. A galaxy can be dynamically stable according to criterion (36) [or even more generally (42)] while it fails to satisfy criterion (38). Therefore, the nonlinear Antonov first law is similar to a notion of ensembles inequivalence between microcanonical and canonical ensembles in thermodynamics (Chavanis 2006a).

We shall study the linear dynamical stability of a spatially homogeneous stationary solution of the Vlasov equation described by a distribution function $f=f({\bf v})$. Like for an infinite homogeneous gas, we make the Jeans swindle. Linearizing the Vlasov equation around this steady state, taking the Laplace-Fourier transform of this equation and writing the perturbations in the form $e^{i({\bf k}\cdot{\bf r}-\omega t)}$, we obtain the classical dispersion relation (Binney & Tremaine 1987):

$\displaystyle\epsilon({\bf k},\omega)\equiv 1+{4\pi G\over k^{2}}\int{{\bf k}\cdot{\partial f\over\partial{\bf v}}\over{\bf k}\cdot{\bf v}-\omega}d{\bf v}=0,$

(49)

where $\epsilon({\bf k},\omega)$ is similar to the “dielectric function” of plasma physics (with the sign $+$ instead of $-$). For a given distribution $f({\bf v})$, this equation determines the complex pulsation(s) $\omega=\omega_{r}+i\omega_{i}$ of a perturbation with wavevector ${\bf k}$. Since the time dependence of the perturbation is $\delta f\sim e^{-i\omega_{r}t}e^{\omega_{i}t}$, the system is linearly stable if $\omega_{i}<0$ and linearly unstable if $\omega_{i}>0$. The condition of marginal stability corresponds to $\omega_{i}=0$.

If we take the wavevector ${\bf k}$ along the $z$-axis¹²¹²12If the distribution function is isotropic, there is no restriction in making this choice., the dispersion relation becomes

$\displaystyle\epsilon({k},\omega)\equiv 1+{4\pi G\over k^{2}}\int_{C}{{k}{\partial f\over\partial{v_{z}}}\over{k}{v_{z}}-\omega}d{v}_{x}d{v}_{y}d{v}_{z}=0,$

(50)

where the integral must be performed along the Landau contour $C$ (Binney & Tremaine 1987). We define

$\displaystyle g(v_{z})=\int fdv_{x}dv_{y}.$

(51)

In the following, we shall note $v$ instead of $v_{z}$ and $f$ instead of $g$. With these conventions, the dispersion relation (50) can be rewritten

$\displaystyle\epsilon({k},\omega)\equiv 1+{4\pi G\over k^{2}}\int_{C}{f^{\prime}(v)\over{v}-\frac{\omega}{k}}d{v}=0.$

(52)

This is the fundamental equation of the problem. For future reference, let us recall that

$\displaystyle\int_{C}{{f^{\prime}(v)}\over v-\frac{\omega}{k}}dv=\int_{-\infty}^{+\infty}{{f^{\prime}(v)}\over v-\frac{\omega}{k}}dv,\quad(\omega_{i}>0),$

(53)

	$\displaystyle\int_{C}{{f^{\prime}(v)}\over v-\frac{\omega}{k}}dv=P\int_{-\infty}^{+\infty}{{f^{\prime}(v)}\over v-\frac{\omega}{k}}dv+i\pi f^{\prime}\left(\frac{\omega}{k}\right),\quad(\omega_{i}=0),$
			(54)

	$\displaystyle\int_{C}{{f^{\prime}(v)}\over v-\frac{\omega}{k}}dv=\int_{-\infty}^{+\infty}{{f^{\prime}(v)}\over v-\frac{\omega}{k}}dv+2\pi if^{\prime}\left(\frac{\omega}{k}\right),\quad(\omega_{i}<0),$
			(55)

where $P$ denotes the principal value.

In general, the dispersion relation (52) cannot be solved explicitly to obtain $\omega(k)$ except in some very simple cases ¹³¹³13There exists less simple cases were explicit solutions can be obtained. We should mention in this respect the extensive set of closed-form solutions for generalized Lorentzian distributions due to Summers & Thorne (1991) and Thorne & Summers (1991). They correspond to Tsallis (polytropic) distributions of the form (158) with negative index $n<-1$.. For example, for cold systems described by the distribution function $f=\rho\delta({v}-{v}_{0})$, we obtain after an integration by parts

$\displaystyle\omega=v_{0}k\pm i\sqrt{4\pi G\rho}.$

(56)

In particular, when $v_{0}=0$, we get

$\displaystyle\omega=\pm i\sqrt{4\pi G\rho}.$

(57)

The system is unstable to all wavenumbers and the perturbation grows with a growth rate $\omega_{i}=\sqrt{4\pi G\rho}$. We also note that the dispersion relation (57) coincides with the dispersion relation (17) of a self-gravitating gas with $c_{s}=0$.

For $\omega_{i}=0$, the real and imaginary parts of the dielectric function $\epsilon(k,\omega_{r})=\epsilon_{r}(k,\omega_{r})+i\epsilon_{i}(k,\omega_{r})$ are given by

$\displaystyle\epsilon_{r}({k},\omega_{r})=1+\frac{4\pi G}{k^{2}}P\int_{-\infty}^{+\infty}{f^{\prime}(v)\over{v}-{\omega_{r}}/{k}}d{v},$

(58)

$\displaystyle\epsilon_{i}({k},\omega_{r})=\frac{4\pi^{2}G}{k^{2}}f^{\prime}\left({\omega_{r}}/{k}\right).$

(59)

The condition of marginal stability corresponds to $\epsilon(k,\omega)=0$ and $\omega_{i}=0$, i.e. $\epsilon_{r}({k},\omega_{r})=\epsilon_{i}({k},\omega_{r})=0$. We obtain therefore the equations

$\displaystyle 1+\frac{4\pi G}{k^{2}}P\int_{-\infty}^{+\infty}{f^{\prime}(v)\over{v}-{\omega_{r}}/{k}}d{v}=0,$

(60)

$\displaystyle f^{\prime}\left({\omega_{r}}/{k}\right)=0.$

(61)

The second condition (61) imposes that the phase velocity $\omega_{r}/k=v_{ext}$ is equal to a velocity where $f(v)$ is extremum ($f^{\prime}(v_{ext})=0$). The first condition (60) then determines the wavenumber(s) $k_{c}$ corresponding to marginal stability. It can be written

$\displaystyle 1+\frac{4\pi G}{k_{c}^{2}}\int_{-\infty}^{+\infty}{f^{\prime}(v)\over{v}-v_{ext}}d{v}=0,$

(62)

where the principal value is not needed anymore. The wavenumber(s) corresponding to marginal stability are therefore given by

$\displaystyle k_{c}=\left(-4\pi G\int_{-\infty}^{+\infty}\frac{f^{\prime}(v)}{v-v_{ext}}\,dv\right)^{1/2}.$

(63)

Finally, the pulsation(s) corresponding to marginal stability are $\omega_{r}=v_{ext}k_{c}$ and we have $\delta f\sim e^{-i\omega_{r}t}$. Note that the distribution $f(v)$ can be relatively arbitrary. There can be pure oscillations ($\omega=\omega_{r}\neq 0$) only if $f(v)$ has some extrema at $v\neq 0$. If $f(v)$ has a single maximum at $v=0$, then $\omega_{r}=0$ (implying $\omega=0$) and the condition of marginal stability becomes

$\displaystyle k_{c}=\left(-4\pi G\int_{-\infty}^{+\infty}\frac{f^{\prime}(v)}{v}\,dv\right)^{1/2}.$

(64)

To determine whether the distribution $f=f(v)$ is stable or unstable for a perturbation with wavenumber $k$, one possibility is to solve the dispersion relation (52) and determine the sign of the imaginary part of the complex pulsation. This can be done analytically in some simple cases (see Secs. 3.8, 4 and 5). We can also apply the Nyquist method introduced in plasma physics. This is a graphical method that does not require to solve the dispersion relation. The details of the method are explained by Nicholson (1992) and we just recall how it works in practice. In the $\epsilon$-plane, taking $\omega_{i}=0$, we plot the Nyquist curve¹⁴¹⁴14This curve is also called an hodograph. $(\epsilon_{r}(k,\omega_{r}),\epsilon_{i}(k,\omega_{r}))$ parameterized by $\omega_{r}$ going from $-\infty$ to $+\infty$ (for a given wavenumber $k$). This curve is closed and always rotates in the counterclockwise sense. If the Nyquist curve does not encircle the origin, the system is stable (for the corresponding wavenumber $k$). If the Nyquist curve encircles the origin one or more times, the system is unstable. The number $N$ of tours around the origin gives the number of zeros of $\epsilon(k,\omega)$ in the upper half $\omega$-plane, i.e. the number of unstable modes with $\omega_{i}>0$. The Nyquist method by itself does not give the growth rate of the instability.

Let us consider the asymptotic behavior of $(\epsilon_{r}(k,\omega_{r}),\epsilon_{i}(k,\omega_{r}))$ defined by Eqs. (58) and (59) for $\omega_{r}\rightarrow\pm\infty$. Since $f(v)$ is positive and tends to zero for $v\rightarrow\pm\infty$, we conclude that $\epsilon_{i}(k,\omega_{r})\rightarrow 0$ for $\omega_{r}\rightarrow\pm\infty$ and that $\epsilon_{i}(k,\omega_{r})>0$ for $\omega_{r}\rightarrow-\infty$ while $\epsilon_{i}(k,\omega_{r})<0$ for $\omega_{r}\rightarrow+\infty$. On the other hand, integrating by parts in Eq. (58), we obtain

$\displaystyle\epsilon_{r}(k,\omega_{r})=1+\frac{4\pi G}{k^{2}}{P}\int_{-\infty}^{+\infty}{{f(v)}\over(v-\omega_{r}/k)^{2}}dv,$

(65)

provided that $f(v)$ decreases sufficiently rapidly. Therefore, for $\omega_{r}\rightarrow\pm\infty$, we obtain at leading order

$\displaystyle\epsilon_{r}(k,\omega_{r})\simeq 1+\frac{4\pi G\rho}{\omega_{r}^{2}},\qquad(\omega_{r}\rightarrow\pm\infty).$

(66)

In particular, $\epsilon_{r}(k,\omega_{r})\rightarrow 1^{+}$ for $\omega_{r}\rightarrow\pm\infty$. From these results, we conclude that the behavior of the Nyquist curve close to the limit point $(1,0)$ is like the one represented in Fig. 3. In addition, according to Eq. (59), the Nyquist curve crosses the $x$-axis at each value of $\omega_{r}/k$ corresponding to an extremum of $f(v)$. For $\omega_{r}/k=v_{ext}$, where $v_{ext}$ is a velocity at which the distribution is extremum $(f^{\prime}(v_{ext})=0)$, the imaginary part of the dielectric function $\epsilon_{i}(k,kv_{ext})=0$ and the real part of the dielectric function

$$\epsilon_{r}(k,kv_{ext})=1+\frac{4\pi G}{k^{2}}\int_{-\infty}^{\infty}\frac{f^{\prime}(v)}{v-v_{ext}}\,dv.$$

(67)

Subtracting the value $f^{\prime}(v_{ext})=0$ in the numerator of the integrand, and integrating by parts, we obtain

$$\epsilon_{r}(k,kv_{ext})=1-\frac{4\pi G}{k^{2}}\int_{-\infty}^{\infty}\frac{f(v_{ext})-f(v)}{(v-v_{ext})^{2}}\,dv.$$

(68)

If $v_{Max}$ denotes the velocity corresponding to the global maximum of the distribution, we clearly have

	$\displaystyle\epsilon_{r}(k,kv_{Max})=1-\frac{4\pi G}{k^{2}}\int_{-\infty}^{\infty}\frac{f(v_{Max})-f(v)}{(v-v_{Max})^{2}}\,dv<1.$
			(69)

Furthermore, $\epsilon_{r}(k,kv_{Max})<0$ for sufficiently small $k$. Therefore, by tuning $k$ appropriately, we can always make the Nyquist curve encircle the origin. We conclude that a spatially homogeneous self-gravitating system is always unstable to some wavelengths.