© ESO, 2013

1. Introduction

Core collapse in an N-body system is a problem relevant to the dynamics of globular clusters, because these are considered to be the only stellar systems that possess the necessary degree of collisionality to relax thermally. The discovery of some clusters with cuspy cores was interpreted as a sign that they indeed experienced the gravothermal catastrophe. In the context of globular clusters, the three main phenomena caused by collisionality are: core collapse, evaporation, and mass segregation. The focus of this paper is on core collapse.

Consider the gravitational N-body problem, where N ≫ 1. That is, consider a set of N classical point masses, each of mass m, mutually interacting through Newtonian gravity. The particles may or may not be confined within a spherical volume of radius R. We are interested in the following questions:

• Can the system reach some sort of equilibrium? Can thisequilibrium be called thermal?

• Can a model (for example, collisionless or fluid) of the N-body problem reach equilibrium? How does the equilibrium of a model relate to the equilibrium of the pure N-body problem and of a real stellar system?

• Are these equilibria stable?

This topic has been studied by many authors (for recent reviews see for example Heggie & Hut 2003; Binney & Tremaine 2008). Different models of the N-body problem (in particular gaseous models, collisionless models, and fluid models) admit equilibrium configurations that are spatially truncated self-gravitating isothermal spheres, but it is not clear to what extent they can be considered as representative of the equilibrium states of the pure N-body problem. Several studies have addressed the stability problem of isothermal spheres using thermodynamical methods, starting with Antonov (1962) and Lynden-Bell & Wood (1968) and continuing with a long list of papers (for example, Hachisu & Sugimoto 1978; Nakada 1978; Katz 1978; Inagaki 1980; Padmanabhan 1989; Chavanis 2002, 2003, see also Thirring 1970). In the thermodynamical approach the study of isothermal spheres is based on a form of entropy known as the Boltzmann entropy: ${\mathit{S}}_{\mathrm{b}}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\equiv }\mathrm{-}\mathit{k}\mathrm{\int }\mathit{f}\mathrm{\left(}{r}\mathit{,}{v}\mathrm{\right)}\ln \mathit{f}\mathrm{\left(}{r}\mathit{,}{v}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{v,}$(1)where f is the one-particle distribution function¹, r and v the position and velocity vector respectively, and k the Boltzmann constant.

Unfortunately, the thermodynamics of self-gravitating systems still depends on a number of unresolved issues, partly because of the long-range nature of the force and partly because of the divergent behavior of the force at short distance (e.g., see Padmanabhan 1990; Katz 2003; Chavanis 2006; Mukamel 2008; Campa et al. 2009; Bouchet et al. 2010). A critical analysis of the use of the Boltzmann entropy has been made by Miller (1973). We briefly comment on this point in Sect. 2.

Therefore, it would be interesting to study the gravothermal catastrophe by limiting the use of thermodynamical arguments as much as possible. In this paper we reconsider the gravothermal catastrophe and analyze the dynamical stability of self-gravitating fluids governed by the Euler equation by finding the normal modes and frequencies of spherical systems under various external conditions: constant total energy E and volume V (this case corresponds to the gravothermal catastrophe of Lynden-Bell & Wood 1968), constant temperature T and volume V (isothermal collapse), constant temperature T and boundary pressure P (isobaric collapse, Bonnor 1956; Ebert 1955). In contrast to thermodynamical arguments, the linear modal analysis automatically gives the time scale for the development of the instability. The constant {T,V} case has been addressed by Semelin et al. (2001), who studied the stability of the system numerically, and by Chavanis (2002), who found an analytical solution for the marginally stable perturbations using methods developed by Padmanabhan (1989). The constant {T,P} case has been addressed previously with synthetic arguments by Bonnor (1956) and Ebert (1955), by Yabushita (1968) who studied the stability of the system numerically, by Lombardi & Bertin (2001) who extended the analysis by Bonnor (1956) to the nonspherically symmetric case, and by Chavanis (2003) who gave an analytical solution for the case of marginally stable perturbations using the same method as he had used for the constant {T,V} case. The constant { E,V } case has been studied using a model based on the Smoluchowski-Poisson system (different from the Euler-Poisson system considered in this paper), by Chavanis et al. (2002).

In this paper, we extend these analyses of the constant {T,V} and {T,P} cases to the general calculation of eigenfrequencies and eigenfunctions for conditions outside those of marginal stability, and study the constant {E,V} case in detail. We use a Eulerian or Lagrangian representation of hydrodynamics as suggested by the boundary conditions to be imposed on the system and provide a unified treatment for all cases. Surprisingly, we find that the system becomes dynamically unstable in all the cases considered exactly at the same points as found by Lynden-Bell & Wood (1968) by means of a thermodynamical analysis, that is, for values of the density contrast (the ratio of the central density to the boundary density) ≈14, 32.1, and 709 respectively for the constant {T,P}, {T,V}, and {E,V} cases. These results suggest that the gravothermal catastrophe may reduce to the Jeans (1902) instability, rediscovered in the inhomogeneous context. Then we briefly outline the problem of the linear dynamical stability of a spherical fluid generalized to the two-component case, by perturbing the two-component spatially truncated isothermal sphere configurations previously considered by Taff et al. (1975), Lightman (1977), Yoshizawa et al. (1978), de Vega & Siebert (2002), who studied the stability of these systems with a thermodynamical approach, and Sopik et al. (2005), who also addressed the dynamical stability with a different model from the one used in this paper. We show that the onset of thermodynamical and dynamical stability occurs at the same point in the simplest case (constant {T,V}) and confirm that the component made of heavier particles is the primary driver of the instability. This result is likely to be related to a similar finding by Breen & Heggie (2012a,b) for two-component gravothermal oscillations. Finally we focus on the following puzzling phenomenon: the instability disappears in the collisionless model. The phenomenon is at variance with what happens in homogeneous systems (e.g., see Bertin 2000). In the last part of the paper we argue that the cause of the difference lies in the role of the angular momentum of the individual particles.

The paper is organized as follows. In Sect. 2 we comment on the meaning of the Boltzmann entropy. In Sect. 3 we write the basic hydrodynamic equations in the Eulerian and Lagrangian representation. In Sect. 4 we show the results of the linear analysis of the hydrodynamic equations by studying the properties of spherically-symmetric perturbations. We find the relevant normal modes and frequencies under different boundary conditions. In the analysis of the gravothermal catastrophe, we propose a modified expression for the energy, which can be considered as a simple way to incorporate the energy generation from binaries in the linear regime; we show that the catastrophe can indeed be halted if we consider such a modified expression. In the last section and in Appendix E, we briefly consider the two-component case. In Sect. 5.1 we discuss the relevant time scales. In Sect. 5.2 we discuss the difference between the collisionless model and the fluid model of the N-body problem. In Sect. 6 we draw our conclusions and identify some open questions and issues.

2. Thermodynamical approach and Boltzmann entropy

Antonov (1962) and Lynden-Bell & Wood (1968) based their stability analysis of self-gravitating spatially-truncated isothermal spheres on the use of the Boltzmann entropy (1). Starting from a kinetic description, they looked for stationary states of the Boltzmann entropy with respect to the distribution function f, at fixed² total mass M, total energy E, and volume V = 4πR³/3, where R is the radius of a spherical box. These stationary states are spatially truncated isothermal spheres. It is found that, for a given single-particle mass m, each spatially truncated isothermal sphere is identified by two dimensional scales (e.g., the central density ρ₀(0) and the temperature T) and one dimensionless parameter. A possible choice for the dimensionless parameter is Ξ = R/λ, where λ =  [kT/m4πGρ₀(0)] ^1/2 is reminiscent of the Jeans length (Jeans 1902); another equivalent choice is the density contrast ρ₀(0)/ρ₀(R), that is, the ratio of the central density to the density at the truncation radius R. Conversely, each choice of the three quantities {T,ρ(0),Ξ } identifies one particular isothermal sphere. This one-to-one correspondence may be lost for other choices of the quantities characterizing the system.

In this thermodynamical approach the necessary condition for stability is that the stationary state corresponds to a (local at least) maximum of the Boltzmann entropy functional. Antonov (1962) found that if Ξ < 34.4 (i.e., if ρ₀(0)/ρ₀(R) < 709) the equilibrium configurations are local maxima, whereas they are saddle points for Ξ > 34.4. Therefore, he concluded that self-gravitating isothermal spheres are unstable if Ξ > 34.4.

The work of Antonov (1962) was extended by Lynden-Bell & Wood (1968). These authors studied the thermodynamical stability of the system under various conditions by applying the relevant thermodynamical potentials (based on the Boltzmann entropy) and gave a physical interpretation of the instability in terms of negative specific heats, as a characteristic feature of self-gravitating systems, thus creating the paradigm of the gravothermal catastrophe. They found that the critical value of Ξ that corresponds to the onset of instability depends on the adopted conditions. For example, Ξ = 34.4 holds for a system at constant total energy E and volume V, Ξ = 8.99 for a system at constant temperature T and volume V, and Ξ = 6.45 for a system at constant temperature T and external pressure P. We show that the three points can be identified by means of a dynamical stability analysis of the fluid model.

We briefly comment here on the use of the Boltzmann expression for the entropy (1). If we refer to systems for which the traditional thermodynamic limit (N → ∞ and V → ∞ at N/V fixed) is well defined and leads to homogeneous equilibria, it is well known (Jaynes 1965) that the Boltzmann entropy corresponds to the entropy defined in phenomenological thermodynamics only when the interparticle forces do not affect the thermodynamic properties; that is, the Boltzmann entropy neglects the interparticle potential energy and the effect of the interparticle forces on the pressure. If the equation of state is different from that of an ideal gas, the Boltzmann entropy is in error by a non-negligible amount. The correct expression for the entropy (corresponding to phenomenological thermodynamics) of systems with a well-defined thermodynamic limit is the one given by Gibbs (for the definition of Gibbs entropy, see Jaynes 1965).

The generalization of the above result to systems with inhomogeneous equilibria, such as self-gravitating systems, is that the Boltzmann entropy is valid (i.e., it corresponds to entropy as defined in phenomenological thermodynamics) if and only if (1) the equation of state of the ideal gas applies locally to the equilibrium states; (2) hydrostatic equilibrium holds. Indeed, the most elementary way to understand the Boltzmann entropy in the case of self-gravitating systems is to think of a fluid in hydrostatic equilibrium, with the equation of state of an ideal gas p = ρkT/m, and consider it subject to reversible transformations between equilibrium states (i.e., between different truncated isothermal spheres in the spherically symmetric case). As shown by Lynden-Bell & Wood (1968) in their Appendix I, the classical thermodynamic entropy defined from the relation dS_c = dQ/T and calculated from such transformations coincides with S_b. For the general case of a system of particles interacting through an arbitrary two-body potential (not necessarily gravitational), it is possible to show that the stationary states of the Boltzmann entropy have the same density distribution as that of a fluid with the equation of state of an ideal gas in hydrostatic equilibrium.

The above considerations suggest that the validity of the Boltzmann entropy is strictly related to the validity of the equation of state of an ideal gas. In a kinetic description, the equation of state of an ideal gas is obtained by considering the local one-particle distribution function to be a Maxwellian (i.e., of the form f(r,v) = Aexp [ − mv²/2kT(r)], where A is a normalization constant) and by defining the pressure as p ≡ ∫mf(v²/3)   d³r   d³v. This definition of pressure, which is obtained by considering particles that would reverse their momentum when hitting an imaginary wall, ignores the effects of interparticle forces. If, for example, particles repelling one another were confined in a box, they would exert some pressure on the walls of the box even if at rest. This contribution is completely neglected by the equation of state of the ideal gas. The neglected pressure is similar to that of rigid spheres when packed too closely, thus the Boltzmann entropy cannot give a correct result for a gas of almost rigid spheres when their density is too high.

In the case of particles interacting through gravity, the neglected contribution is attractive and should therefore decrease the pressure compared to that of an ideal gas. Therefore, it is unlikely that the true thermal equilibrium state of an N-body system when t → ∞ (t is time) has an effective equation of state³ that is not affected by the attractive nature of the gravitational force: thus the Boltzmann entropy may not be applicable to find the true thermal equilibrium of a self-gravitating N-body system, since it assumes the equation of state of an ideal gas. Moreover, because gravitational forces are long-ranged, each particle feels the influence of all other distant particles. Thus, it may be that, strictly speaking, in the case of self-gravitating systems an equation of state cannot be defined in terms of local quantities.

A different way to justify the use of the Boltzmann entropy is to show that it can be obtained from the Gibbs microcanonical entropy by means of the so-called mean field approximation (Katz 2003; Padmanabhan 1990). However, the range of applicability of this approximation is still not clear, also because a small-scale cut-off is needed to avoid divergences in the microcanonical entropy (Padmanabhan 1990). The Boltzmann entropy is also the H quantity of the Boltzmann H-theorem; from this point of view, a critical analysis of the Boltzmann entropy in the context of the gravothermal catastrophe has been performed by Miller (1973).

The discussion in this section supports the hypothesis that self-gravitating isothermal spheres are not true thermal equilibrium states of the pure N-body problem (i.e., the states that would be found in a spherical box containing N gravitating particles eventually, after an infinite amount of time), but are only metastable states, and that the significance of these metastable states is still not completely understood (for example, see Padmanabhan 1990; Chavanis 2006, and references therein).

3. Basic equations of the dynamical approach

In this section we derive the linearized hydrodynamic equations that govern the evolution of a fluid system for small deviations from the truncated isothermal sphere equilibrium configurations. We assume spherical symmetry; that is, we only consider radial perturbations. The configurations are known to be stable against nonradial perturbations (Semelin et al. 2001; Chavanis 2002; Binney & Tremaine 2008).

Consider a self-gravitating fluid governed by the Navier-Stokes and continuity equations, together with the equation of state of an ideal gas (for example, see Landau & Lifshitz 1987): $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{c}\end{array}\frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{u}\mathrm{·}{\nabla }\mathrm{)}{u}\mathrm{=}\mathrm{-}\frac{{\nabla }\mathit{p}}{\mathit{\rho }}\mathrm{-}{\nabla }\mathrm{\Phi }\mathrm{+}\frac{\mathit{\eta }}{\mathit{\rho }}{\mathrm{\nabla }}^{\mathrm{2}}{u}\mathrm{+}\frac{\mathit{\zeta }\mathrm{+}\frac{\mathit{\eta }}{\mathrm{3}}}{\mathit{\rho }}{\nabla }\mathrm{(}{\nabla }\mathrm{·}{{u}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{c}\end{array}\mathit{p}\mathrm{=}\mathit{\rho }\frac{\mathit{kT}}{\mathit{m}}\mathit{,}\end{array}$(2)where ρ is the density, u the fluid velocity, p the local pressure, Φ the gravitational potential, η and ζ are viscosity coefficients, m is the one-particle mass, T the temperature, t time, and k the Boltzmann constant. We keep track of the viscosity terms until the end of Sect. 3.1; then (from Sect. 3.2 on), for simplicity, we set η = ζ = 0. In this paper, we do not perform an analysis of the effects of viscosity. In particular we do not discuss whether it has a stabilizing or destabilizing effect: hopefully, the role of viscosity will be studied in detail in a subsequent paper.

Under the assumption of spherical symmetry, the gravitational potential obeys the following equations: $\begin{array}{ccc}{\nabla }\mathrm{\Phi }& \mathrm{=}& \frac{\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathit{r̂}\mathit{,}\\ \mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& {\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}\mathit{\rho }\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{4}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{s,}\end{array}$(3)which are equivalent to the Poisson equation.

The hydrostatic equilibria of such a fluid system are spatially truncated isothermal spheres (Chandrasekhar 1967) and are briefly described in Appendix A. As mentioned in Sect. 2, each self-gravitating truncated isothermal sphere is identified by two dimensional scales (for example the central density ρ₀(0) and the temperature T) and one dimensionless parameter Ξ = R/λ, which is the value of the dimensionless radius ξ ≡ r/λ at the truncation radius. The scale λ was defined at the beginning of Sect. 2. Conversely, each choice of the three quantities { ρ₀(0),T,Ξ } identifies a particular isothermal sphere.

We denote unperturbed quantities by subscript 0 and the perturbations by subscript 1. Thus the unperturbed density profile (the density profile of a truncated isothermal sphere) is, as a function of the dimensionless radius ξ (see Appendix A), $\{\begin{array}{c}\\ & \mathrm{if}\mathit{\xi }\mathrm{\le }\mathrm{\Xi }\\ & \mathrm{if}\mathit{\xi }\mathit{>}\mathrm{\Xi }\mathit{,}\end{array}$(4)where ρ₀(0) is the unperturbed central density and ψ(ξ) is the regular solution to the Emden equation (the symbol ′ denotes derivative with respect to the argument ξ): $\frac{\mathrm{d}}{\mathrm{d}\mathit{\xi }}\mathrm{\left(}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}{\mathit{\xi }}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\mathit{,\psi }\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0.}$(5)When we perturb the hydrodynamic equations it is convenient to work in the Eulerian or Lagrangian representation of hydrodynamics, depending on which boundary conditions we impose. In the Eulerian representation the independent variables are the time t and the position vector r. In the Lagrangian representation the independent variables are the time t and the original position r₀ that the fluid element under consideration had at the initial time t₀.

3.1. Eulerian representation

In this section we derive the linearized perturbation equations in the Eulerian representation, so we write each quantity as the sum of an unperturbed part (characterizing the truncated isothermal sphere) and a perturbation: $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\\ \mathit{u}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& {\mathit{u}}_{\mathrm{1}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\\ \mathit{T}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& {\mathit{T}}_{\mathrm{0}}\mathrm{+}{\mathit{T}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}\end{array}$(6)where u is the radial component of the fluid velocity. Recall that we only consider radial perturbations. We substitute Eq. (6)into Eq. (2)and expand to first order in the perturbed quantities. We allow the temperature to vary in time while remaining uniform in space: this constraint will be used to impose the condition of constant total energy (see Sect. 4.1.2). To find the normal modes of the system, we look for solutions of the following form: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \begin{array}{c}\mathrm{˜}\\ {\mathit{\rho }}_{\mathrm{1}}\end{array}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega t}}\\ {\mathit{u}}_{\mathrm{1}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \begin{array}{c}\mathrm{˜}\\ {\mathit{u}}_{\mathrm{1}}\end{array}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega t}}\\ {\mathit{T}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{1}}\end{array}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega t}}\mathit{.}\end{array}$(7)In the following, we drop the symbol ~ to keep the notation simpler. After some manipulations (see Appendix B.1) and using the unperturbed density profile (4)we obtain the following linearized equation: $\begin{array}{ccc}\mathit{Lf}& \mathrm{=}& \mathrm{ℒ}\mathit{f}\mathrm{+}\frac{\mathrm{i}\mathit{\omega }}{\mathrm{4}\mathit{\pi G\lambda }}\frac{\mathit{k}{\mathit{T}}_{\mathrm{1}}}{\mathit{m}}{\mathit{\psi }}^{\mathrm{\prime }}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\\ & \mathrm{+}& \mathrm{i}\mathit{\omega }\frac{\mathit{m}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{k}{\mathit{T}}_{\mathrm{0}}}\{\mathit{\eta }\frac{\mathrm{1}}{{\mathit{\xi }}^{\mathrm{2}}}\mathrm{\left[}{\mathit{\xi }}^{\mathrm{2}}\mathrm{\right(}\mathit{f}{\mathrm{e}}^{\mathit{\psi }}{\mathrm{)}}^{\mathrm{\prime }}{\mathrm{]}}^{\mathrm{\prime }}\mathrm{+}\mathrm{(}\mathit{\zeta }\mathrm{+}\frac{\mathit{\eta }}{\mathrm{3}}\mathrm{)}{\left[\frac{\mathrm{1}}{{\mathit{\xi }}^{\mathrm{2}}}{\mathrm{(}{\mathit{\xi }}^{\mathrm{2}}\mathit{f}{\mathrm{e}}^{\mathit{\psi }}\mathrm{)}}^{\mathrm{\prime }}\right]}^{\mathrm{\prime }}\}\mathit{,}\end{array}$(8)where f(ξ) = ρ₀(ξ)u₁(ξ) is the unknown function, $\mathit{L}\mathrm{=}\frac{{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$(9)represents the dimensionless (squared) eigenfrequency and ℒ is the following differential operator: $\mathrm{ℒ}\mathrm{\equiv }\mathrm{-}\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\mathit{\xi }}^{\mathrm{2}}}\mathrm{-}\left(\frac{\mathrm{2}}{\mathit{\xi }}\mathrm{+}{\mathit{\psi }}^{\mathrm{\prime }}\right)\frac{\mathrm{d}}{\mathrm{d}\mathit{\xi }}\mathrm{+}\left(\frac{\mathrm{2}}{{\mathit{\xi }}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{2}{\mathit{\psi }}^{\mathrm{\prime }}}{\mathit{\xi }}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\right)\mathit{.}$(10)Equation (8)is the equation to be solved to find the normal modes and frequencies of the system, under the appropriate boundary conditions and constraints. The boundary conditions for Eq. (8)are defined in the following way. The absence of sinks and sources of mass, together with the assumption of spherical symmetry, implies f(0) = 0. Since in the Eulerian representation we only consider systems in a spherical box of constant volume, the radial velocity at the edge must be zero, which implies f(Ξ) = 0. Thus the boundary conditions are, $\mathit{f}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathit{f}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}\mathrm{=}\mathrm{0.}$(11)To solve the equations, we still need to specify the function T₁(t). This specification distinguishes between the constant energy and the constant temperature case. Once this specification is made, Eq. (8)is an eigenvalue problem: for fixed ω, the equation admits a solution only for discrete values of L. These solutions are the radial normal modes of the system.

The operator ℒ has the following properties, which can be proved directly from the Emden Eq. (5): $\begin{array}{ccc}\mathrm{ℒ}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}{\mathit{\psi }}^{\mathrm{\prime }}\\ \mathrm{ℒ}\mathrm{\left(}\mathit{\xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}{\mathit{\psi }}^{\mathrm{\prime }}\mathit{.}\end{array}$(12)These properties allow us to obtain analytical solutions in some cases. They are equivalent to those found by Padmanabhan (1989) in his review of the thermodynamical analysis of Antonov (1962) and later also used by Chavanis (2002, 2003) to obtain analytical solutions of the present hydrodynamic problem for the constant { T,V } (Sect. 4.1.1) and constant {T,P} (Sect. 4.2.1) cases.

Viscosity disappears when ω = 0, which is the situation of marginal stability. Thus viscosity does not modify the points of the onset of instability (Semelin et al. 2001; Chavanis 2002).

3.2. Lagrangian representation

In this section we present the linearized perturbation equations in the Lagrangian representation in the inviscid case. In this representation, the independent variable is the position r₀ of the fluid element under consideration at the initial time t₀. Thus in Eq. (2)we have to perform the following change of independent variables: $\{\begin{array}{c}\\ \\ \end{array}$(13)In the Lagrangian representation, each quantity is a function of the new independent variables r₀ and t.

The calculations are summarized in Appendix B.2.1. For linear perturbations, the resulting continuity and Euler equations in the Lagrangian representation are (neglecting the term quadratic in the velocity):

$\begin{array}{ccc}& & \left(\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{-}\frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}\mathit{u}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\right)\mathit{\rho }\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathit{\rho }{\mathit{u}}^{\mathrm{)}}\mathrm{=}\mathrm{0}\\ & & \left(\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{-}\frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}\mathit{u}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\right)\mathit{u}\mathrm{=}\mathrm{-}\frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{\partial \rho }}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{kT}}{\mathit{m}}\mathrm{-}\frac{\mathit{GM}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{·}\end{array}$(14)

As for the Eulerian case we separate each quantity in an unperturbed part and a perturbation: $\begin{array}{c}\\ \mathit{\rho }\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\mathrm{+}{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\\ \mathit{u}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{u}}_{\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\\ \mathit{r}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{r}}_{\mathrm{0}}\mathrm{+}{\mathit{r}}_{\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\\ \mathit{T}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{0}}\mathrm{+}{\mathit{T}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}\end{array}$(15)We substitute Eq. (15)in Eqs. (14)and expand to first order in quantities with subscript 1. Then to find the normal modes we take: $\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{\rho }}_{\mathrm{1}}\end{array}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega t}}\\ {\mathit{u}}_{\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{u}}_{\mathrm{1}}\end{array}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega t}}\\ {\mathit{r}}_{\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{r}}_{\mathrm{1}}\end{array}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega t}}\mathit{.}\end{array}$(16)In the following we shall drop the symbol ~ for simplicity of notation. After introducing the dimensionless radius ξ₀ ≡ r₀/λ, using the density profile (4)of the unperturbed state and after some manipulations we obtain the following equation for ρ₁(ξ₀) (see Appendix B.2.2 for an outline of the calculations): $\mathrm{-}\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}\left[\frac{\left(\frac{\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}\mathrm{\left(}\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\right)}}{\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}}}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{1}}}{{\mathit{T}}_{\mathrm{0}}}\right)}{\frac{\mathrm{4}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{3}}}\mathrm{+}\frac{\mathit{L}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}}}\right]\mathrm{=}\mathrm{0}\mathit{,}$(17)where L = ω²/4πGρ₀(0) as for the Eulerian case. Boundary conditions are discussed in Sect. 4.2.1. Solving Eq. (17)allows us to find the normal modes in the Lagrangian representation. We have allowed the temperature to vary in time in Eq. (17), but in the following we only analyze the isothermal case with T₁ = 0.

Fig. 1 Effective potential U (see Eqs. (20)and (21)) for the constant { T,V } case. To find the frequencies of normal modes, the Schrödinger Eq. (20)has to be solved with boundary conditions (22), which means that the effective potential is U(ξ) for ξ ≤ Ξ and taken to be infinite for ξ > Ξ. States with negative energies, which exist for Ξ ≥ 8.99, correspond to unstable modes.

4. Modal stability of a self-gravitating inviscid fluid sphere under different boundary conditions

Here we analyze the equations obtained in the previous section by imposing different boundary conditions.

4.1. Eulerian representation

In this section we analyze Eq. (8)by imposing two kinds of boundary conditions: constant temperature T and volume V, or constant total energy E and volume V.

4.1.1. Constant {T, V} case (isothermal collapse)

Here we consider a fluid at constant temperature T contained in a sphere of fixed radius R. The condition of constant temperature is satisfied by imposing T₁ = 0 in Eq. (8). The condition of constant volume has been discussed in Sect. 3.1 and leads to the boundary condition f(Ξ) = 0. Thus Eq. (8)for the constant { T,V } case becomes $\begin{array}{ccc}\mathrm{ℒ}\mathit{f}\mathrm{=}\mathit{Lf}& & \end{array}$(18)with the boundary conditions (11).

Equation (18)is an eigenvalue equation that, at given Ξ, admits solutions only for discrete values of L. If the lowest value of L at fixed Ξ is positive, then all modes are stable (because ω is always real) and the system is stable. If the lowest value of L at a given value of Ξ is negative, then unstable modes are present and the system is unstable.

By means of the standard transformation $\mathit{f}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\exp \left[\mathrm{-}{\mathrm{\int }}_{\mathit{\xi ̅}}^{\mathit{\xi }}\left(\frac{\mathrm{1}}{{\mathit{\xi }}^{\mathrm{\prime }}}\mathrm{+}\frac{{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\xi }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathrm{2}}\right)\mathrm{d}{\mathit{\xi }}^{\mathrm{\prime }}\right]\mathit{,}$(19)where is an arbitrary point in the domain of f, Eq. (18)can be recast in the form of a Schrödinger equation: $\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{f}}^{\mathrm{\prime \prime }}\end{array}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathit{U}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{=}\mathit{L}\begin{array}{c}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathit{,}$(20)where U(ξ) is the effective potential given by $\mathit{U}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{\equiv }\frac{\mathrm{2}}{{\mathit{\xi }}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{(}\mathit{\xi }{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{2}}{\mathit{\xi }}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{(}\mathit{\xi }\mathrm{)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathit{.}$(21)The boundary conditions become $\begin{array}{c}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathit{f}\end{array}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}\mathrm{=}\mathrm{0.}$(22)The effective potential is shown in Fig. 1. From the boundary conditions (22)we see that choosing a specific Ξ requires us to consider a potential that is infinite for ξ ≥ Ξ. The existence of negative eigenvalues implies that the system is unstable. From the form of the potential it is clear that for low values of Ξ negative eigenvalues (i.e., unstable modes) do not exist. Negative eigenvalues only appear for sufficiently high values of Ξ. For the present problem this occurs at Ξ ≥ 8.99. When Ξ → ∞, an infinite number of negative eigenvalues appear, precisely at the same points where new unstable modes appear in the thermodynamical approach (see Katz 1978).

Fig. 2 Minimum value of the eigenvalue L at given Ξ (the dimensionless radius characterizing the system) as a function of Ξ, for the constant { T,V } case. If the minimum L is negative then the system is unstable. The system becomes unstable at Ξ = 8.99, the same value obtained from the thermodynamical approach. From the figure we can read the typical time scale of the instability. As Ξ → ∞, the asymptotic value L∞ =  −0.042 is the same as for the other cases (see Figs. 3 and 4).

We now consider the solutions to Eq. (18)in greater detail. Figure 2 shows the minimum value of L at given dimensionless radius Ξ as a function of Ξ. We see that L is negative, that is, the system is unstable, for Ξ > 8.99, which is the same point found by Lynden-Bell & Wood (1968) in the thermodynamical approach. Higher modes, that is, higher values of L at given Ξ, would be represented by lines above the plotted curve. These lines would intersect the Ξ axis at some points, which are the zeros of the analytical solution G_TV described below (see Eq. (23)). In Appendix D.1 a few density and velocity profiles of numerically calculated eigenfunctions are shown. Most of them are for modes of minimum L at given Ξ. Density profiles of higher modes exhibit oscillations not present in the lowest mode.

For the case of marginal stability (L = 0), with the help of properties (12), the relevant eigenfunction can be expressed analytically (Chavanis 2002). The function ${\mathit{G}}_{\mathrm{TV}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{\equiv }{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{-}\mathit{\xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}$(23)is indeed a solution to Eq. (18)with L = 0, which satisfies G_TV(0) = 0. The values of Ξ for which G_TV(Ξ) = 0 are those for which the boundary conditions (11)are satisfied, so they are the values of Ξ at which each new unstable mode appears. The first zero of G_TV occurs where the first unstable mode appears, at Ξ = 8.99. From the asymptotic behavior of ψ it is possibile to obtain an asymptotic approximation of the zeros: they approximately follow a geometric progression of ratio ${\mathrm{e}}^{\mathrm{2}\mathit{\pi }\mathit{/}\sqrt{\mathrm{7}}}$ (see also Semelin et al. 2001; Chavanis 2002).

4.1.2. Constant {E, V} case (gravothermal catastrophe)

Here we consider a self-gravitating fluid sphere at constant total energy E and volume V. In this case the instability has been named gravothermal catastrophe by Lynden-Bell & Wood (1968). The total energy E of the fluid is defined as $\mathit{E}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{NkT}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\int }\mathit{\rho }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{\Phi }\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}\mathit{.}$(24)It has two terms, which represent the thermal and gravitational contributions. The condition of constant energy is imposed in the following way. When the fluid is perturbed, its gravitational energy changes as a consequence of the redistribution of matter. We suppose that the temperature varies in time, while remaining uniform in space, so as to keep the total energy (thermal plus gravitational) constant. For a different model, based on the Smoluchowski-Poisson system of equations, the same nonstandard assumption has been made by Chavanis et al. (2002). The thermal energy expression at time t is thus given by 3NkT(t)/2. In doing so, we are assuming infinite thermal conductivity (see also Sect. 5.1 for the relation of this fact to the relevant time scales).

To reduce Eq. (8)to the constant { E,V } case we need to find the expression for the temperature as a function of the density distribution at fixed total energy, in the linear regime of small perturbations. Starting from Eq. (24)for the total energy and recalling that $\mathrm{\Phi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{G}{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{3}}{\mathit{r}}^{\mathrm{\prime }}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathrm{|}{r}\mathrm{-}{{r}}^{\mathrm{\prime }}\mathrm{|}$, we have $\mathit{E}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{NkT}\mathrm{-}\frac{\mathit{G}}{\mathrm{2}}\mathrm{\int }\frac{\mathit{\rho }\mathrm{\left(}{r}\mathrm{\right)}\mathit{\rho }\mathrm{\left(}{{r}}^{\mathrm{\prime }}\mathrm{\right)}}{\mathrm{|}{r}\mathrm{-}{{r}}^{\mathrm{\prime }}\mathrm{|}}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{\mathit{r}}^{\mathrm{\prime }}\mathit{.}$(25)By substituting T = T₀ + T₁, ρ = ρ₀ + ρ₁ in Eq. (25), only keeping first-order quantities and imposing that the energy remains constant, we obtain the following expression for T₁: ${\mathit{T}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{\mathrm{\int }{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}{r}\mathrm{\right)}{\mathrm{\Phi }}_{\mathrm{0}}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}}{\frac{\mathrm{3}}{\mathrm{2}}\mathit{Nk}}\mathrm{·}$(26)Here Φ₀ is the gravitational potential of the unperturbed density distribution ρ₀. After some manipulations (see Appendix B.3) we obtain ${\mathit{T}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{i}\mathit{\omega }}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Xi }}\mathrm{\left\{}\mathrm{d}\mathrm{\right[}{\mathit{\xi }}^{\mathrm{2}}\mathit{f}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{]}\mathit{/}\mathrm{d}\mathit{\xi }\mathrm{\}}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\xi }}{\frac{\mathrm{3}}{\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{\lambda }\hspace{0.17em}{\mathrm{\Xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}{\mathit{T}}_{\mathrm{0}}\mathit{.}$(27)By substituting Eq. (27)in Eq. (8), recalling the definition λ = [kT₀/4πGρ₀(0)m]^1/2, and neglecting viscosity, we obtain $\mathit{Lf}\mathrm{=}\mathrm{ℒ}\mathit{f}\mathrm{-}{\mathit{\psi }}^{\mathrm{\prime }}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\frac{\mathrm{1}}{\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\Xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Xi }}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\frac{\mathrm{d}}{\mathrm{d}\mathit{\xi }}\mathrm{\left[}{\mathit{\xi }}^{\mathrm{2}}\mathit{f}\mathrm{\right(}\mathit{\xi }{\mathrm{\left)}}^{\mathrm{\right]}}\hspace{0.17em}\mathrm{d}\mathit{\xi }\mathit{.}$(28)As discussed in Sect. 3.1 the boundary conditions are given by Eq. (11). By integrating the last term of Eq. (28)by parts under the boundary conditions (11), we obtain $\mathit{Lf}\mathrm{=}\mathrm{ℒ}\mathit{f}\mathrm{+}{\mathit{\psi }}^{\mathrm{\prime }}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\frac{\mathrm{1}}{\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\Xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Xi }}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathit{f}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\xi }\mathit{.}$(29)Equation (29)(or (28)) contains an integral global constraint. This (29)is to be solved to find the normal modes and corresponds to Eq. (18). The numerical procedure followed to solve Eq. (29)is relatively straightforward and is thus not reported here.

Figure 3 shows the minimum value of L at given Ξ. We see that the system is unstable for Ξ > 34.36, which corresponds to a density contrast ρ₀(0)/ρ₀(Ξ) > 709 (Antonov 1962). As Ξ increases, new unstable modes appear, similar to the behavior observed in the constant { T,V } case. As Ξ → ∞, the asymptotic value of the minimum L is found numerically to be the same as in the constant { T,V } case, L_∞ ≃  − 0.042. Such asymptotic value is the also same in the constant { T,P } case (see Sect. 4.2.1).

Fig. 3 Minimum value of the eigenvalue L at given Ξ, the dimensionless radius characterizing the system, as a function of Ξ, for the constant { E,V } case. When the minimum L is negative the system is unstable. The system becomes unstable at Ξ = 34.36, which corresponds to a density contrast ≃709, the same value obtained from the thermodynamical approach. From the figure we can read the typical time scale of the instability. As Ξ → ∞, the asymptotic value of L, L∞ =  −0.042, is the same as for the other cases (see Figs. 2 and 4).

In the case of marginal stability (L = 0), Eq. (29)is equivalent to what is found and solved analytically by Padmanabhan (1989) in the thermodynamical approach. Here, for completeness, we record how the solution is derived by adapting the method of Padmanabhan (1989) to our choice of variables and unknowns.

To analytically solve Eq. (29)for L = 0 we rewrite it in the following form: $\mathrm{ℒ}\mathit{f}\mathrm{=}\mathrm{-}\mathit{C}{\mathit{\psi }}^{\mathrm{\prime }}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\mathit{,}$(30)where $\mathit{C}\mathrm{=}\frac{\mathrm{1}}{\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\Xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Xi }}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathit{f}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\xi }$(31)is a constant (with respect to ξ) that is characterized by a global dependence on f. Using the properties (12), we look for a solution of the form ${\mathit{G}}_{\mathrm{EV}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{=}\mathit{a}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{+}\mathit{b\xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathit{,}$(32)where a and b are real numbers. Since Eq. (30)is linear in f, only the ratio a/b is relevant. Susbtituting (32)into (30)we obtain the first condition on a and b: $\mathit{a}\mathrm{+}\mathit{b}\mathrm{=}\mathit{C}\mathit{.}$(33)Using the expression of C (31)and dividing by b, we obtain $\frac{\mathit{a}}{\mathit{b}}\mathrm{+}\mathrm{1}\mathrm{=}\frac{\mathrm{1}}{\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\Xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Xi }}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\left(\frac{\mathit{a}}{\mathit{b}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{+}\mathit{\xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\right)\hspace{0.17em}\mathrm{d}\mathit{\xi }\mathit{.}$(34)A second condition on a/b follows from the boundary condition f(Ξ) = 0: $\frac{\mathit{a}}{\mathit{b}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}\mathrm{+}\mathrm{\Xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}\mathrm{=}\mathrm{0.}$(35)Substituting a/b from (35)in (34)we obtain $\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\Xi }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\right(}\mathrm{\Xi }\mathrm{)}\mathrm{-}\mathrm{\Xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\Xi }}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\left(\mathit{\xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathrm{-}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\frac{\mathrm{\Xi }{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}}{{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}}\right)\hspace{0.17em}\mathrm{d}\mathit{\xi }\mathit{.}$(36)The values of Ξ satisfying Eq. (36)are those for which Eq. (29)admits a solution for L = 0. In particular, the lowest value of Ξ that satisfies Eq. (36)determines the threshold of instability: by numerically solving the algebraic Eq. (36), this minimum value is found to be Ξ = 34.36 (as shown by Antonov 1962). As for the constant { T,V } case, the number of unstable modes for each Ξ coincides with the results of the thermodynamical analysis (see Katz 1978). Once a value of Ξ is obtained, the value of a/b and then an analytical solution G_EV is determined.

Density profiles for modes of minimum L at given Ξ are shown in Appendix D.2. Since the equations involved are equivalent, the density profile for the marginally stable perturbation (L = 0) is the same as found by Padmanabhan (1989). He pointed out that it has a “core-halo” structure, that is, an oscillation in ρ₁: the density perturbation is positive in the inner part (nucleus), negative in the middle (emptying area), and then positive again (halo). The core-halo structure has been physically interpreted in the framework of the gravothermal catastrophe given by Lynden-Bell & Wood (1968), using the concept of negative specific heat. However, this interpretation is not applicable to the present context.

We found that for modes of minimum L at given Ξ the core-halo structure is present if L < 0.021, disappears between L = 0.021 and L = 0.022, and is absent for L > 0.022, as illustrated in the figures presented in Appendix D.2. Higher modes always exhibit one ore more oscillations, as in the constant { T,V } case.

We have thus shown that for the present case the behavior of the eigenvalues L as a function of Ξ is similar to the constant { T,V } case, with the difference that the instability threshold is higher in the constant { E,V } case. The interpretation of this in the context of the fluid model is as follows. When the fluid is compressed, the gravitational potential energy decreases and the temperature increases in order to maintain the total energy constant. Therefore, the tendency toward collapse is weakened, and instability can take place only at higher values of Ξ (than when the temperature remains constant).

In fact, if instead of expression (24)for the total energy we consider a modified expression in such a way that the temperature increase is greater, the collapse can be halted completely. Consider the following heuristic expression for the total energy: $\begin{array}{ccc}\mathit{E}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{NkT}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\int }\mathit{\rho }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{\Phi }\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}\mathrm{-}\mathit{\upsilon }{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{3}}\mathit{\rho }\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{,}& & \end{array}$(37)which contains an additional term proportional to the dimensionless parameter υ and to the central density ρ(0). The quantity ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}\mathit{k}{\mathit{T}}_{\mathrm{0}}\mathit{/}\mathit{m}$ is the square of the unperturbed thermal speed. We repeated the analysis of this section with such a modified expression for the energy and found the new value of Ξ for the threshold of the instability. In this analysis the expression for T₁ (27)obtained previously, to be substituted in Eq. (8), is replaced with the expression for T₁ that is obtained from Eq. (37)by substituting T = T₀ + T₁, ρ = ρ₀ + ρ₁ and by imposing that the variation in the total energy E is zero. We found that when υ > 0 (i.e., when the fluid is compressed the new term contributes to make the temperature increase) the (linear) instability is postponed to higher values of Ξ for small υ and is completely halted for υ > 5 × 10^-4. We also verified that if υ < 0 (i.e., when the fluid is compressed the new term contributes in the opposite direction), instability occurs at lower values of Ξ.

In the case of the pure N-body problem, it is believed that collapse can be halted by energy “generation” through binaries, giving rise to a phenomenon called gravothermal oscillations (for a review, see Heggie & Hut 2003). Gravothermal oscillations were discovered by Sugimoto & Bettwieser (1983) using a gaseous model. These authors introduced in the model of Lynden-Bell & Eggleton (1980) a phenomenological energy-generation term (to represent the role of binaries), proportional to a power of the density and found that collapse can be halted and reversed. Similarly, our υ term shows in a simple manner that the instability can be halted in the linear regime with a suitable term in the total energy budget that mimics effects presumed to be associated with the presence of binaries.

Gravothermal oscillations have been confirmed by N-body simulations (Makino 1996). We comment further on this topic in Sect. 5.2.

4.1.3. Constant {T, V}: two-component case

We briefly studied the problem of dynamical stability by means of a linear modal analysis of a two-component ideal fluid. Each component is assumed to interact with the other only through the common gravitational potential. The hydrostatic equilibrium configurations are spatially truncated two-component isothermal spheres, as considered by Taff et al. (1975), Lightman (1977), and Yoshizawa et al. (1978), see also de Vega & Siebert (2002) and Sopik et al. (2005).

We call the single-particle masses m_A and m_B, with m_B > m_A. Two-component isothermal spheres are characterized by two additional dimensionless parameters (with respect to the one-component case): the ratio of the single-particle masses m_B/m_A and the ratio M_B/M_A of the total masses associated with the two components. The reader is referred to Appendix E for a description of the equations used.

We considered the constant { T,V } case, which is the simplest from the mathematical point of view, and found that the equivalence of the thermodynamical and dynamical approaches still holds. It is possible to show analytically that the onset of dynamical instability takes place at exactly the same points as those found with the thermodynamical approach. The analysis in the thermodynamical approach was performed by generalizing in a straightforward manner the analysis of Chavanis (2002).

An interesting result of this analysis is that the instability appears to be driven by the heavier component, in the following sense. Consider the ratios ρ_1A/ρ₀ and ρ_1B/ρ₀ of the density perturbation of each component to the total local unperturbed density. The ratio referring to the heavier component can be higher even if the total mass of the heavier component is very small. For example, for m_B/m_A = 3, we found that the ratios ρ_1A/ρ₀, ρ_1B/ρ₀ were comparable for M_B/M_A ≃ 0.066 (see Fig. E.1). For higher values of M_B/M_A the heavier component dominates. Independent indications that the heavier component dominates the collapse have been found by Sopik et al. (2005) for a different model based on the Smoluchowski-Poisson system of equations. Breen & Heggie (2012a,b) found that the heavier component dominates gravothermal oscillations in two-component clusters. This is likely to be related to the present analysis.

4.2. Lagrangian representation

4.2.1. Constant {T, P} case (isobaric collapse)

In this section we analyze perturbations at constant temperature T and constant boundary pressure P. Therefore, T₁ = 0 and Eq. (17)becomes: $\mathrm{-}\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}\left[\frac{\frac{\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}\mathrm{\left(}\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\right)}}{\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}}}}{\frac{\mathrm{4}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{3}}}\mathrm{+}\frac{\mathit{L}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}}}\right]\mathrm{=}\mathrm{0.}$(38)Boundary conditions are as follows. The condition u₁(0) = 0 translates into the condition ${\mathit{\rho }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$ (using Euler Eq. (14), under the assumption that ∂u/∂r₀ does not diverge and ∂u/∂t =  −iωu). The condition of constant pressure requires that a fixed Lagrangian fluid shell (which follows the fluid during the motion and thus does not have a fixed position in space) feels constant pressure. Constant pressure requires constant density because of the ideal gas equation of state. Therefore, the relevant boundary conditions are: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \mathrm{0}\\ {\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}& \mathrm{=}& \mathrm{0.}\end{array}$(39)From mass continuity and by imposing ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\ne \mathrm{0}$ (which is true in the Eulerian representation and thus we expect to be true in the present case), we also have the condition ρ₁(0) ≠ 0.

The numerical procedure to solve Eq. (38)is relatively straightforward and is thus not reported here. Figure 4 represents the minimum value of L at given Ξ. The system becomes unstable for Ξ > 6.45, which is the same condition as found by Bonnor (1956), Ebert (1955), and Lynden-Bell & Wood (1968). As Ξ increases, other unstable modes appear, similarly to the cases described previously. As Ξ → ∞, the asymptotic value of the minimum L is found numerically to be the same as for the constant { E,V } and { T,V } cases, L_∞ =  −0.042, but, in contrast to the constant { E,V } and { T,V } cases, it is reached from below. There is a minimum at L ≈  − 0.043 for Ξ ≈ 15.

Fig. 4 Minimum value of the eigenvalue L at given Ξ, the dimensionless radius characterizing the system, as a function of Ξ, for the constant { E,V } case. When the minimum L is negative the system is unstable. The system becomes unstable at Ξ = 6.45, the same value as found in the thermodynamical approach. From the figure we can read the typical time scale of the instability. As Ξ → ∞, the asymptotic value L∞ =  −0.042 is the same as for the other cases (see Figs. 2 and 3).

In Appendix D.3 density profiles of normal modes are shown. Most of them are for modes of minimum L at given Ξ. Higher modes present oscillations. It would be interesting to show whether the analysis of this section might be generalized to the nonspherically-symmetric case by following Lombardi & Bertin (2001).

As for previous cases, we can obtain analytical solutions for the marginally stable perturbations. For L = 0, Eq. (38)reads as $\mathrm{ℒ}\left(\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\right)\mathrm{=}\mathrm{0}\mathit{,}$(40)where ℒ is defined as $\begin{array}{ccc}\mathrm{ℒ}\mathrm{\equiv }\mathrm{-}\mathrm{1}\mathrm{+}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}\left[\frac{\frac{\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}}{\frac{\mathrm{d}}{\mathrm{d}{\mathit{\xi }}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}}}{\frac{\mathrm{4}}{{\mathit{\xi }}_{\mathrm{0}}^{\mathrm{3}}}}\right]\mathrm{=}\mathrm{0.}& & \end{array}$(41)From Emden Eq. (5), the operator ℒ has the following properties: $\begin{array}{ccc}\mathrm{ℒ}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{\prime }\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}}{\mathrm{2}}\\ \mathrm{ℒ}\mathrm{\left(}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}}{\mathrm{4}}\mathrm{·}\end{array}$(42)Hence, an analytical solution of Eq. (40)is the following (for a derivation in the Eulerian representation, see also Chavanis 2003): ${\mathit{G}}_{\mathrm{TP}}\mathrm{\left(}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\mathrm{-}{\mathit{\psi }}^{\mathrm{\prime }\mathrm{2}}\mathit{.}$(43)This solution satisfies the boundary conditions ρ₁(0) = 2 and ${\mathit{\rho }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$, as required by Eq. (39). The zeros of G_TP(ξ₀) allow us to identify the marginally stable normal modes that satisfy the correct boundary conditions (39). The first zero of G_TP(ξ₀) is at Ξ = 6.45. Other zeros correspond to the values of Ξ at which new unstable modes appear and are found to be the same as in the standard thermodynamical approach.

5. Different behavior of a collisionless self-gravitating sphere

5.1. Time scales

In this section we briefly discuss the typical time scales that characterize gaseous and fluid models and compare them to those of globular clusters. To a large extent, our discussion follows that of Inagaki (1980).

Let τ_LTE be the local relaxation time, τ_GTE the global relaxation time, and τ_d the dynamical time scale. For a gaseous system, we identify τ_LTE with the time (close to the inverse mean collision frequency) needed to reach a local Maxwellian distribution, τ_GTE with the time in which thermal conduction balances the temperatures of different parts of the system, and τ_d with the sound travel time. Gaseous models generally assume the following ordering: ${\mathit{\tau }}_{\mathrm{LTE}}\mathrm{\ll }{\mathit{\tau }}_{\mathrm{d}}\mathrm{\ll }{\mathit{\tau }}_{\mathrm{GTE}}\mathit{.}$(44)For a globular cluster, because the mean free path is very large and stars can cross the system many times before actually colliding, the ordering of time scales is different. Stars do not collide significantly with neighboring stars, following the mechanisms that usually characterize thermal conduction; instead, they tend to release their energy through the entire cluster (also because of the long range nature of the force). Thus for globular clusters we have ${\mathit{\tau }}_{\mathrm{d}}\mathrm{\ll }{\mathit{\tau }}_{\mathrm{LTE}}\mathrm{\simeq }{\mathit{\tau }}_{\mathrm{GTE}}\mathit{.}$(45)In the fluid model analyzed in this paper we assumed infinite thermal conductivity, because the temperature was always taken to be and to remain uniform. Moreover, we implicitly assumed that the distribution function is locally Maxwellian. Therefore, our fluid model follows the ordering: $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\tau }}_{\mathrm{LTE}}\mathrm{\ll }\mathit{\tau }d\\ \begin{array}{c}\end{array}{\mathit{\tau }}_{\mathrm{GTE}}\mathrm{\ll }{\mathit{\tau }}_{\mathrm{d}}\mathit{.}\end{array}$(46)Assumptions (44)and (46)are not applicable to the situation of a globular cluster (45). This difference and the related limitations must be kept in mind if we wish to apply these models to understand the evolution of globular clusters.

5.2. Dynamics of a collisionless self-gravitating isothermal sphere

Isothermal spheres are equilibrium configurations for several idealized models of the pure N-body problem. In particular, self-gravitating isothermal spheres can be studied as stationary states of the collisionless Boltzmann equation: $\frac{\mathit{\partial f}\mathrm{\left(}{r}\mathit{,}{v}\mathit{,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}\frac{\mathit{\partial f}\mathrm{\left(}{r}\mathit{,}{v}\mathit{,t}\mathrm{\right)}}{\mathit{\partial }{r}}\mathrm{-}\frac{\mathit{\partial }\mathrm{\Phi }\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}}{\mathit{\partial }{r}}\mathrm{·}\frac{\mathit{\partial f}\mathrm{\left(}{r}\mathit{,}{v}\mathit{,t}\mathrm{\right)}}{\mathit{\partial }{v}}\mathrm{=}\mathrm{0.}$(47)Therefore, it is natural to ask whether such collisionless isothermal spheres are or can be unstable. It can be shown (Binney & Tremaine 2008) that the unbounded collisionless isothermal sphere is linearly stable, in contrast to the results of the fluid counterpart⁴. It is generally believed that the collisionless isothermal sphere is also stable in the nonlinear regime, although to our knowledge a rigorous proof of this statement is still lacking. Moreover, if only spherically symmetric perturbations are considered, it is possible to show, through an argument based on the conservation of the detailed angular momentum⁵ (Appendix C), that a collisionless sphere (bounded or unbounded) cannot collapse, because each particle has a minimum radius it can attain. In contrast, the fluid system analyzed in Sects. 3 and 4 is only unstable with respect to spherically symmetric perturbations, with the subsequent nonlinear evolution presumably leading to a collapse. The above considerations support the hypothesis that self-gravitating collisionless isothermal spheres (bounded or unbounded) are stable with respect to all kinds of perturbations and cannot collapse (Kandrup & Sygnet 1985; Kandrup 1990; Batt et al. 1995).

The different behavior, between the collisionless and the fluid cases, emerged in this paper is actually quite surprising, because for the homogeneous case the collisionless and fluid models behave in the same way with respect to Jeans instability (e.g., see Bertin 2000), in the sense that both the fluid and the collisionless systems are linearly unstable under the same criterion for instability.

If the self-gravitating collisionless isothermal sphere were unstable, its instability would develop on the dynamical time scale; in turn, it is commonly believed that the gravothermal catastrophe may only occur if the system is at least weakly collisional, and thus it is thought to develop on the collision time scale.

In the spherically symmetric case a difference between the two models is the following. While in the fluid model each fluid element is sustained against gravity by pressure of the inner parts, in the collisionless model stars are sustained by their individual angular momentum relative to the center, that is, by their velocity dispersion⁶. In moving from a kinetic description to a fluid description, all the information about microscopic velocities and detailed angular momentum is lost: in particular, each fluid element has zero angular momentum (see also Appendix C). In this respect, the real situation of a globular cluster resembles the collisionless case more; strictly speaking, stars are not sustained by pressure, but rather by their velocity dispersion, because the mean free paths are long, and stars cross the cluster many times before feeling the effects of collisions.

If we come back to the original pure N-body problem, it is therefore natural to ask: what is, in this case, the role of detailed angular momentum? We may argue that a quantity related to the detailed angular momentum should characterize the process of core collapse. If the collapse can happen in an N-body system, it may be accompanied by a slow decrease in the sum of the magnitudes of the angular momenta of the individual particles slowly sinking toward the center.

For a particle with angular momentum J that moves in a gravitational potential generated by a mass M a radius related to angular momentum (which is the radius of the circular orbit if M is a point mass) can be defined as $\mathit{d}\mathrm{\equiv }\frac{{\mathit{J}}^{\mathrm{2}}}{\mathit{G}{\mathit{m}}^{\mathrm{2}}\mathit{M}}\mathrm{·}$(48)Therefore, for a stellar system of N stars and total mass M we can define a radius related to detailed angular momentum as ${\mathit{R}}_{\mathrm{am}}\mathrm{\equiv }\frac{\mathit{A}}{\mathit{NG}{\mathit{m}}^{\mathrm{2}}\mathit{M}}\mathit{,}$(49)where A is the sum of the squares of the angular momenta of the individual stars, that is, the following quantity: $\mathit{A}\mathrm{\equiv }\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{{J}}_{\mathit{i}}^{\mathrm{2}}\mathit{.}$(50)Thus we may argue that the typical time scale for core collapse should correlate with ${\mathit{t}}_{\mathrm{cc}}\mathrm{\equiv }\frac{\mathit{A}}{\mathrm{d}\mathit{A}\mathit{/}\mathrm{d}\mathit{t}}\mathrm{·}$(51)The main mechanism through which A varies with time is expected to be that of two-body collisions. Thus t_cc should be similar to the two-body relaxation time.

By monitoring the quantity A(t) in N-body simulations, it would be interesting to test the role of detailed angular momentum in the mechanism of gravothermal oscillations (Makino 1996) (A(t) may reach an equilibrium value, around which the system gravothermally oscillates. This would be consistent with the fact that the typical time scale of gravothermal oscillations is the two-body relaxation time.), its relevance to the studies of core collapse in the gaseous model (Lynden-Bell & Eggleton 1980; Sugimoto & Bettwieser 1983), and its connection with the phenomenon of the gyro-gravothermal catastrophe (Hachisu 1979).

6. Discussion and conclusions

In this paper we have studied systematically the dynamical stability of a self-gravitating isothermal fluid sphere by means of a linear modal analysis with respect to spherically symmetric perturbations. In this sense, we studied the Jeans instability in the inhomogeneous context of a sphere of finite size. Within a unified framework, by imposing the boundary conditions of constant { T,V } (isothermal collapse), constant { E,V } (gravothermal catastrophe), and constant { T,P } (isobaric collapse), we have proved that the onset of dynamical instability occurs exactly at the points identified in the thermodynamical approach (see Antonov 1962; Lynden-Bell & Wood 1968; Bonnor 1956; Ebert 1955) and, by adapting derivations from other authors (Padmanabhan 1990; Chavanis 2002, 2003), we provided an analytic expression for the eigenfunctions of the marginally stable modes. Indeed, as noted in the Introduction, some results along these lines have been obtained previously by other authors. The main new results obtained in this paper are the following:

• Using the fluid model based on the Euler equation, we extendedprevious studies of the constant {T,V} and {T,P} cases to the constant { E,V } case (gravothermal catastrophe), proving that the onset of Jeans instability occurs exactly at the same point identified in the thermodynamical approach also in this case. For this constant { E,V } case, we introduced a heuristic term to incorporate effects akin to the stabilizing role of binaries.

• For all the three cases described above, we calculated numerically eigenfrequencies and eigenfunctions of the relevant modes also outside the conditions of marginal stability. The time scale for the instability that we found is the dynamical time scale. The excitation of higher modes has been illustrated in a simple way, by referring to an effective potential that governs the structure of the linear modal analysis.

• We found that for all the cases treated in our investigation, as the dimensionless radius of the isothermal sphere Ξ becomes larger and larger, the value of the dimensionless growth rate of the most unstable mode tends to a universal asymptotic constant value, independent of the adopted boundary conditions.

• We briefly showed that the correspondence between the stability in the dynamical and in the thermodynamical approach also holds for the two-component case and found indications that the heavier component is the more important driver of the instability.

• As a general discussion, we commented on the meaning and applicability of the Boltzmann entropy for self-gravitating systems and argued that the main difference between the dynamical behavior of a fluid and a collisionless sphere, in relation to their application as models of the pure N-body problem or a real weakly collisional stellar system, such as globular clusters, is to be ascribed to the role of the detailed angular momentum behavior in the collisionless and weakly collisional cases.

The role of the viscosity and its consequences outside the condition of marginal stability have not been examined; hopefully, this issue will be addressed in a future paper.

Another interesting question is how the instability depends on the particle-particle interaction, for non-Newtonian cases (Padmanabhan 1989). Potentials that exhibit a softening at small radii, such as $\mathrm{1}\mathit{/}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, where r₀ is a constant (Chavanis & Ispolatov 2002; see also Casetti & Nardini 2012), or that decline with a different power law at large radii, such as 1/r^α, with α ≠ 1 (Ispolatov & Cohen 2001), have indeed been considered. Based on the present article, we may argue that the results obtained from the thermodynamical approach would be reinterpreted and clarified as the analog of the Jeans instability, by studying a fluid model for the potential considered, with the equation of state of a perfect gas.

In general, this paper strengthens the view that the applicability of different idealized models to describe the process of core collapse in systems made of a finite number of stars is more subtle than commonly reported and that, in general, the study of Jeans instability of inhomogeneous stellar systems still leaves a number of questions open.

Online material

Appendix A: Fluid truncated isothermal spheres

In this appendix we summarize the properties of spatially truncated self-gravitating fluid isothermal spheres.

The equations for the hydrostatic equilibrium of a spherically symmetric fluid with the equation of state of an ideal gas are $\begin{array}{ccc}\frac{\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{d}{\mathit{p}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{d}\mathit{r}}\mathit{,}\\ \mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{r}\\ \mathrm{0}\end{array}\mathrm{4}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{d}\mathit{s,}\\ {\mathit{p}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{kT}\mathit{/}\mathit{m}\mathit{.}\end{array}$We express M(r) by means of the condition of hydrostatic equilibrium (A.1), differentiate it to obtain dM, and equate this to dM = 4πρ₀(r)r². By making the change of variable ρ₀(r) = ρ₀(0)e^− ψ(r), where ρ₀(0) is the central density, and introducing the dimensionless radius ξ = r/λ, where λ = [kT/4πGρ₀(0)m]^1/2, we obtain the differential equation for ψ(ξ), recorded in the main text as Eq. (5). (Because the constant ρ₀(0) is interpreted as the central density, we are considering the boundary condition ψ(0) = 0; the density is taken to be regular at the origin, so that the second boundary condition is ψ′(0) = 0.) The solution of Eq. (5)(called Emden equation) is a monotonic increasing function characterized by logarithmic behavior ψ(ξ) ~ lnξ² and ψ′(ξ) ~ 2/ξ as ξ → ∞.

From Eq. the mass enclosed within the radius ξ is $\mathit{M}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{=}\frac{\mathit{kT\lambda }}{\mathit{Gm}}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathit{.}$(A.4)From Eq. (A.4)and the asymptotic behavior of ψ, it is clear that a solution with finite total mass is only obtained by truncating the system at a dimensionless radius ξ = Ξ. A truncated isothermal sphere is then identified by two scales T and ρ₀(0) and one dimensionless parameter Ξ. The density profile of a spatially truncated isothermal sphere is given by ρ₀(r) = ρ₀(0)e^− ψ(r) where ψ is the solution to Eq. (5).

Appendix B: Linearization of the hydrodynamic equations

Appendix B.1: Eulerian representation

Here we record the calculations leading to the linearized Eq. (8). The unperturbed density profile is given by Eq. (4). We substitute Eqs. (6)in Eq. (2)and expand to first order in quantities with subscript 1 to obtain $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{1}}}{\mathit{\partial t}}\mathrm{+}{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}}{{u}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }{{u}}_{\mathrm{1}}}{\mathit{\partial t}}\mathrm{=}\left(\frac{{\nabla }{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{-}\frac{{\nabla }{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}}\right)\frac{\mathit{k}{\mathit{T}}_{\mathrm{0}}}{\mathit{m}}\mathrm{-}\frac{{\nabla }{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathit{k}{\mathit{T}}_{\mathrm{1}}}{\mathit{m}}\mathrm{-}\mathrm{4}\mathit{\pi G}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{s,t}\mathrm{\right)}{\mathit{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{s}}{{\mathit{r}}^{\mathrm{2}}}\\ & & \mathrm{+}\frac{\mathit{\eta }}{{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{\nabla }}^{\mathrm{2}}{{u}}_{\mathrm{1}}\mathrm{+}\frac{\mathit{\zeta }\mathrm{+}\frac{\mathit{\eta }}{\mathrm{3}}}{{\mathit{\rho }}_{\mathrm{0}}}{\nabla }\mathrm{(}{\nabla }\mathrm{·}{{u}}_{\mathrm{1}}\mathrm{)}\mathit{.}\end{array}$Then we look for solutions of the form (7). From Eq. (B.1)we obtain ρ₁: ${\mathit{\rho }}_{\mathrm{1}}\mathrm{=}\frac{{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}}{{u}}_{\mathrm{1}}\mathrm{\right)}}{\mathit{i\omega }}\mathit{,}$(B.3)and thus eliminate it from Eq. (B.2)to find the radial component of the Navier-Stokes equation: $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{u}}_{\mathrm{1}}& \mathrm{=}& \left\{\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathit{\partial r}}{{\mathit{\rho }}_{\mathrm{0}}}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{u}}_{\mathrm{1}}\mathrm{)}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{u}}_{\mathrm{1}}\mathrm{)}\right]\right\}\frac{\mathit{k}{\mathit{T}}_{\mathrm{0}}}{\mathit{m}}\\ & & \mathrm{-}\mathrm{i}\mathit{\omega }\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\partial r}}\frac{\mathit{k}{\mathit{T}}_{\mathrm{1}}}{\mathit{m}}\mathrm{-}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{u}}_{\mathrm{1}}\\ & & \mathrm{+}\mathrm{i}\mathit{\omega }\left\{\mathit{\eta }\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}\frac{\mathit{\partial }{\mathit{u}}_{\mathrm{1}}}{\mathit{\partial r}}\right)\mathrm{+}\mathrm{(}\mathit{\zeta }\mathrm{+}\frac{\mathit{\eta }}{\mathrm{3}}\mathrm{)}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}{\mathit{u}}_{\mathrm{1}}\mathrm{\right)}\right]\right\}\mathit{,}\\ & & \end{array}$(B.4)where u₁ is the radial component of the velocity. By defining f(r) ≡ ρ₀(r)u₁(r) and introducing the dimensionless radius ξ = r/λ we obtain Eq. (8).

Appendix B.2: Lagrangian representation

Appendix B.2.1: Change of variables

Here we show the change of variables leading from Eqs. (2)(Eulerian representation) to Eqs. (14)(Lagrangian representation). Let us assume spherical symmetry and neglect viscosity. By dropping the nonlinear term (u·∇)u, the Euler equation and the continuity Eq. (2)become, in the Eulerian representation:

$\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathit{\rho u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial u}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial \rho }\mathit{/}\mathit{\partial r}}{\mathit{\rho }}\frac{\mathit{kT}}{\mathit{m}}\mathrm{-}\frac{\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{·}\end{array}$(B.5)

Now we perform a change of variables. In the Lagrangian representation, each quantity is described as a function of the new independent variables r₀ and t, where r₀ is the position of the fluid element at t = t₀. The standard rules to transform derivatives of a generic function f are: $\begin{array}{ccc}\frac{\mathit{\partial f}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial r}}& \mathrm{=}& \frac{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial r}}\frac{\mathit{\partial f}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\\ \frac{\mathit{\partial f}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial t}}& \mathrm{=}& \frac{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial t}}\frac{\mathit{\partial f}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{\partial f}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{·}\end{array}$(B.6)The partial derivatives of r₀ are obtained from the following relation, which expresses the condition that two fluid shells do not cross each other: ${\mathrm{\int }}_{\mathrm{0}}^{{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{4}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}\mathit{\rho }\mathrm{\left(}\mathit{s,t}\mathrm{\right)}\mathrm{4}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathit{.}$(B.7)By taking the partial derivative with respect to r of Eq. (B.7)(each side of the equation is considered as a function of r, t) we obtain ${\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\frac{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial r}}\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathit{.}$(B.8)By taking the partial derivative of Eq. (B.7)with respect to t we obtain ${\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\frac{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathit{s}}^{\mathrm{2}}\frac{\mathit{\partial \rho }\mathrm{\left(}\mathit{s,t}\mathrm{\right)}}{\mathit{\partial t}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathit{.}$(B.9)From the continuity equation, Eq. (B.9)then becomes ${\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\frac{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{\mathit{r}}^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathit{u}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathit{.}$(B.10)From Eqs. (B.10), (B.8), and (B.6)we obtain the equations of hydrodynamics in the Lagrangian representation (14).

Appendix B.2.2: Linearization

Here we approximate Eqs. (14)to first order for small perturbations around the hydrostatic equilibrium states. We substitute Eqs. (15)in Eqs. (14)and expand to first order in quantities with subscript 1. By noting that M(r₀,t) = M(r₀,t = t₀), we obtain $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{1}}}{\mathit{\partial t}}\mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{u}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{c}\end{array}\frac{\mathit{\partial }{\mathit{u}}_{\mathrm{1}}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathit{k}{\mathit{T}}_{\mathrm{0}}}{\mathit{m}}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\left(\mathrm{2}\frac{{\mathit{r}}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{1}}\mathit{/}\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}{\mathit{\partial }{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathit{\partial }{\mathit{r}}_{\mathrm{0}}}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{1}}}{{\mathit{T}}_{\mathrm{0}}}\right)\mathrm{+}\mathrm{2}\frac{\mathit{GM}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}{\mathit{r}}_{\mathrm{1}}\mathit{.}\end{array}$(B.11)From the usual rules of derivation, we have: $\frac{\mathit{\partial r}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathit{,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathit{/}\mathit{\partial t}}{\mathit{\partial }{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathit{/}\mathit{\partial r}}\mathrm{·}$(B.12)By applying Eqs. (B.10)and (B.8), Eq. (B.12)can be written as $\frac{\mathit{\partial }{\mathit{r}}_{\mathrm{1}}}{\mathit{\partial t}}\mathrm{=}{\mathit{u}}_{\mathrm{1}}\mathit{.}$(B.13)Now we wish to obtain one equation involving only the unknown ρ₁, by combining the three Eqs. (B.11)and (B.13). We assume the modal dependence (16). We eliminate r₁ from (B.13)and the second of (B.11). We then eliminate u₁ from the resulting equation and the first of (B.11), to obtain the following equation {we also used hydrostatic equilibrium to replace $\mathrm{\left[}\mathrm{\right(}\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathrm{d}{\mathit{r}}_{\mathrm{0}}\mathrm{)}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}\mathrm{]}\mathrm{(}\mathit{k}{\mathit{T}}_{\mathrm{0}}\mathit{/}\mathit{m}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{[}\mathit{GM}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\mathit{/}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{]}$}: $\mathrm{-}{\mathit{\rho }}_{\mathrm{1}}\mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}{\mathit{r}}_{\mathrm{0}}}\left[\frac{\mathit{GM}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\left(\frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{1}}\mathit{/}\mathrm{d}{\mathit{r}}_{\mathrm{0}}}{\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathrm{d}{\mathit{r}}_{\mathrm{0}}}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{1}}}{{\mathit{T}}_{\mathrm{0}}}\right)}{\mathrm{4}\frac{\mathit{GM}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}\mathrm{+}{\mathit{\omega }}^{\mathrm{2}}}\right]\mathrm{=}\mathrm{0.}$(B.14)By referring to the dimensionless radius ξ₀ ≡ r₀/λ, we obtain Eq. (17).

Appendix B.3: Temperature expression for the constant {E, V} case

Here we show the steps leading from Eq. (26)to Eq. (27).

From the Poisson and Emden equations, the dimensionless potential ψ is related to the gravitational potential Φ₀ of the unperturbed density distribution in the following way: $\mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\left[{\mathrm{\Phi }}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}{\mathrm{\Phi }}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\right]\mathit{/}\mathrm{(}\mathit{k}{\mathit{T}}_{\mathrm{0}}\mathit{/}\mathit{m}\mathrm{)}\mathit{.}$(B.15)From Eq. (B.15), by recalling that we only consider perturbations that do not change the total mass (${}^{\mathrm{\int }}\mathit{\rho }_{\mathrm{1}}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}\mathrm{=}\mathrm{0}$), we obtain ${\mathit{T}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{\mathrm{\int }{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{r}}{\frac{\mathrm{3}}{\mathrm{2}}\mathit{Nm}}{\mathit{T}}_{\mathrm{0}}\mathit{.}$(B.16)From Eq. (B.3), by setting f = ρ₀u₁ and ξ = r/λ, we obtain Eq. (27).

Appendix C: Detailed angular momentum conservation and collapse

Here, for completeness, we show in detail that the angular momentum barrier prevents a spherically symmetric collisionless system from collapsing. We only consider perturbations that do not break the assumed spherical symmetry of the system.

Consider a collisionless system of particles of individual mass m and total mass M. From the assumption of spherical symmetry, each particle is confined to a plane and we can write the single-particle Lagrangian as $\mathrm{ℒ}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{m}\mathrm{(}\mathit{ṙ}\mathrm{2}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{\theta ̇}{\mathrm{2}}^{\mathrm{)}}\mathrm{-}\mathit{V}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathit{,}$(C.1)where V(r,t) is a time-dependent potential, t time, r the distance from the center, and θ the angular coordinate. The Lagrangian (C.1)conserves the angular momentum of the particle, that is, , where C is a constant. Then the equation of the motion is $\begin{array}{c}\mathrm{¨}\\ \mathit{r}\end{array}\mathrm{=}\frac{{\mathit{C}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{1}}{\mathit{m}}\frac{\mathit{\partial V}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial r}}\mathrm{\equiv }\frac{{\mathit{C}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{-}\frac{\mathit{GM}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}$(C.2)where M(r,t) is the total mass contained in the sphere of radius r. Equation (C.2)is the equation of the motion of a particle moving in one dimension and subject to the force F_r = mC²/r³ − GmM(r,t)/r². The following inequality holds: $\mathrm{0}\mathrm{\le }\mathit{M}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{\le }\mathit{M}\mathit{.}$(C.3)By multiplying by ṙ and integrating both sides of Eq. (C.2), we obtain $\frac{\mathit{ṙ}\mathrm{2}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{2}}\mathrm{=}\frac{\mathit{ṙ}\mathrm{2}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}}\left[\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\right]\mathrm{-}{\mathrm{\int }}_{\mathit{r}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}^{\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\frac{\mathit{GM}\mathrm{\left(}\mathit{s,t}\mathrm{\right(}\mathit{s}\mathrm{\left)}\mathrm{\right)}}{{\mathit{s}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathit{.}$(C.4)Since ṙ²(t)/2 is positive, the righthand side of Eq. (C.5)must be positive. By taking r(t) ≤ r(t₀) (we are not interested in the case in which r(t) is greater than the initial radius) and using Eqs. (C.5)and (C.3), we find $\begin{array}{ccc}\frac{\mathit{ṙ}\mathrm{2}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{2}}& \mathrm{=}& \frac{\mathit{ṙ}\mathrm{2}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}}\left[\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\right]\mathrm{+}{\mathrm{\int }}_{\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}}^{\mathit{r}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}\frac{\mathit{GM}\mathrm{\left(}\mathit{s,t}\mathrm{\right(}\mathit{s}\mathrm{\left)}\mathrm{\right)}}{{\mathit{s}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{s}\\ & \mathrm{\le }& \\ & \mathrm{=}& \frac{\mathit{ṙ}\mathrm{2}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}}\left[\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\right]\mathrm{+}\mathit{GM}\left[\frac{\mathrm{1}}{\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}\right]\mathrm{·}\end{array}$(C.5)For given values of r(t₀) and ṙ(t₀), the quantity appearing in the last line of Eq. (C.5)tends to − ∞ as r(t) → 0. Hence, for given initial conditions the particle cannot reach arbitrarily low values of r(t).

Appendix D: Density and velocity profiles of the linear modes

In this appendix we show density and velocity profiles of the normal modes for the linear stability analysis presented in Sect. 4; ρ₁/ρ₀ and u₁ are meant to be on arbitrary scales.

Appendix D.1: Constant {T, V} profiles

Fig. D.1 Relative density perturbation profiles ρ1(ξ)/ρ0(ξ) and velocity profiles of the normal modes for the constant { T,V } case, obtained by solving Eq. (18). The profiles should be truncated at a value ξ = Ξ where the velocity profile vanishes, to satisfy boundary conditions (11). The vertical dotted line indicates where the system should be truncated to obtain the mode of lowest L for fixed Ξ: for L =  −0.02 and L = 0 only this mode is entirely displayed, while for L = 0.02 two modes are displayed, depending on which zero of the velocity profile is chosen. In the case L = 0, the total density ρ(t) = ρ0 + ρ1(t) at the point ξ = 4.07 remains unchanged, that is, unperturbed; this is one of the relevant points listed by Lynden-Bell & Wood (1968).

Appendix D.2: Constant {E, V} profiles

Fig. D.2 Relative density perturbation profiles ρ1(ξ)/ρ0(ξ) of the normal modes for the constant { E,V } case, obtained by solving Eq. (28). The modes should be truncated at a value ξ = Ξ where the corresponding velocity profile shown in Fig. D.3 vanishes, as marked by the vertical dotted lines, in order to satisfy the boundary conditions (11). Zeros are displayed. Only modes of lowest L at given Ξ are shown. Note that the core-halo structure described in Sect. 4.1.2 disappears between L = 0.021 and L = 0.022.

Fig. D.3 Velocity profiles of the normal modes for the constant { E,V } case, obtained by solving Eq. (28). The modes should be truncated at a value ξ = Ξ corresponding to the vertical dotted lines, in order to satisfy the boundary conditions (11). Other zeros are displayed. Only modes of lowest L for fixed Ξ are shown. Note that the core-halo structure described in Sect. 4.1.2 disappears when the velocity has no internal zeros, between L = 0.021 and L = 0.022.

Appendix D.3: constant {T, P} profiles

Fig. D.4 Relative density perturbation profiles ρ1(ξ0)/ρ0(ξ0) of the normal modes for the constant { T,P } case, obtained by solving Eq. (38), calculated and displayed here in the Lagrangian representation. The modes should be truncated at a value ξ0 = Ξ where the density profile vanishes, in order to satisfy the boundary conditions (39). The first zero, which represents the mode of minimum L at given Ξ, is indicated by the vertical dotted line. A mode of higher L for fixed Ξ is shown in the L = 0.03 case; for other cases higher modes can be identified in a similar way.

Appendix E: Equations for the two-component case

In this appendix we summarize the equations of the linear analysis for the two-component case and show some examples of the density profiles associated with the modes that characterize the onset of the instability. We denote by subscripts A and B the lighter and the heavier component, respectively.

The unperturbed states are the two-component self-gravitating truncated isothermal spheres considered by Taff et al. (1975), Lightman (1977), Yoshizawa et al. (1978), de Vega & Siebert (2002), Sopik et al. (2005). The density profiles can be written as $\{\begin{array}{c}\\ & \mathrm{if}\mathit{\xi }\mathrm{\le }\mathrm{\Xi }\\ & \\ & \mathrm{if}\mathit{\xi }\mathit{>}\mathrm{\Xi }\mathit{,}\end{array}$(E.1)\vspace*{1.2mm}where ρ_A0 and ρ_B0 are respectively the density profiles of the lighter and heavier component, ξ ≡ r/λ₂ is the dimensionless radial coordinate, where λ₂ ≡ [kT(1/m_A + 1/m_B)/4πGρ₀(0)]^1/2 and we denote by ρ₀(ξ) ≡ ρ_A0(ξ) + ρ_B0(ξ) the total unperturbed density; Ξ is the value of ξ at the truncation radius, β ≡ m_B/m_A is the ratio of the single-particle masses, ψ is the solution of the following generalization of the Emden Eq. (5): $\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{d}\mathit{\xi }}\mathrm{(}{\mathit{\xi }}^{\mathrm{2}}{\mathit{\psi }}^{\mathrm{\prime }}\mathrm{)}& \mathrm{=}& {\mathit{\xi }}^{\mathrm{2}}\left[\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{\alpha }}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{)}\mathit{\psi }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{\alpha }}}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{\psi }}\right]\mathit{,}\\ \mathit{\psi }\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& {\mathit{\psi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$where α = ρ_A0(0)/ρ_B0(0) is the ratio of the unperturbed central densities. The symbol ′ denotes derivative with respect to the argument ξ.

The linearized hydrodynamical equations, governing the evolution of the two-component fluid system for small deviations from the unperturbed states described above, are obtained by generalizing in a straightforward manner the steps leading

from Eqs. (2)to Eq. (18). The result, which generalizes Eq. (18), is the following system of equations that governs the evolution of radial perturbations: $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{L}{\mathit{f}}_{\mathrm{A}}\mathrm{=}\left[\mathrm{-}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{\beta }}\right){\mathit{\psi }}^{\mathrm{\prime }}\left(\frac{\mathrm{2}}{\mathit{\xi }}{\mathit{f}}_{\mathrm{A}}\mathrm{+}{\mathit{f}}_{\mathrm{A}}^{\mathrm{\prime }}\right)\mathrm{-}{\mathit{f}}_{\mathrm{A}}^{\mathrm{\prime \prime }}\mathrm{-}\frac{\mathrm{2}}{\mathit{\xi }}{\mathit{f}}_{\mathrm{A}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{2}}{{\mathit{\xi }}^{\mathrm{2}}}{\mathit{f}}_{\mathrm{A}}\right]\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{\beta }}}\\ \begin{array}{c}\end{array}\mathrm{-}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{\alpha }}}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{)}\mathit{\psi }}\mathrm{(}{\mathit{f}}_{\mathrm{A}}\mathrm{+}{\mathit{f}}_{\mathrm{B}}\mathrm{)}\\ \begin{array}{c}\end{array}\mathit{L}{\mathit{f}}_{\mathrm{B}}\mathrm{=}\left[\mathrm{-}\left(\mathrm{1}\mathrm{+}\mathit{\beta }\right){\mathit{\psi }}^{\mathrm{\prime }}\left(\frac{\mathrm{2}}{\mathit{\xi }}{\mathit{f}}_{\mathrm{B}}\mathrm{+}{\mathit{f}}_{\mathrm{B}}^{\mathrm{\prime }}\right)\mathrm{-}{\mathit{f}}_{\mathrm{B}}^{\mathrm{\prime \prime }}\mathrm{-}\frac{\mathrm{2}}{\mathit{\xi }}{\mathit{f}}_{\mathrm{B}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{2}}{{\mathit{\xi }}^{\mathrm{2}}}{\mathit{f}}_{\mathrm{B}}\right]\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{\beta }}\\ \begin{array}{c}\end{array}\mathrm{-}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{\alpha }}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{\psi }}\mathrm{(}{\mathit{f}}_{\mathrm{A}}\mathrm{+}{\mathit{f}}_{\mathrm{B}}\mathrm{)}\mathit{.}\end{array}$(E.4)Here L = ω²/4πGρ₀(0) represents the dimensionless (squared) eigenfrequency, f_A(ξ) ≡ ρ_A0(ξ)u_A1(ξ) and f_B(ξ) ≡ ρ_B0(ξ)u_B1(ξ), where u_A1 and u_B1 are the radial velocity perturbations of the two components. The boundary conditions are: $\begin{array}{c}\\ {\mathit{f}}_{\mathrm{A}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{B}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}\\ {\mathit{f}}_{\mathrm{A}}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{B}}\mathrm{\left(}\mathrm{\Xi }\mathrm{\right)}\mathrm{=}\mathrm{0.}\end{array}$(E.5)Similarly to the one-component case, the two conditions at the center follow from requiring regularity and spherical symmetry, while the two conditions at the truncation radius satisfy the requirement that the radial velocities must vanish at the edge.

The system (E.4)for L = 0 is equivalent to the system that can be obtained by generalizing in a straightforward manner the thermodynamical analysis of Chavanis (2002). The latter analysis can be used to find the points for the onset of instability. This proves that the onset of instability occurs at the same values of Ξ in the dynamical and in the thermodynamical approach.

In Fig. E.1 we show the density profiles for the marginally stable modes (L = 0) in three different situations, that is, with β = 3 and three different values of M_B/M_A. The density perturbation of the heavier component is greater than the density perturbation of the lighter component even for low values of M_B/M_A, indicating that the heavier component is the more important driver of the instability.

Fig. E.1 Relative density perturbation profiles ρ1A(ξ)/ρ0(ξ) (lighter component, dotted line) and ρ1B(ξ)/ρ0(ξ) (heavier component, solid line) of the normal modes for the two-component constant { T,V } case, obtained by solving Eq. (E.4). The plots show marginally stable modes (L = 0) of minimum L at given Ξ. They represent the density profiles that characterize the onset of the instability for fixed value of β = mB/mA = 3 at different values of the total mass ratio MB/MA. Even for low values of MB/MA, the density perturbation of the heavier component dominates, thus suggesting that the heavier component is the more important driver of the instability.