Issue 
A&A
Volume 552, April 2013



Article Number  A37  
Number of page(s)  19  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201220665  
Published online  21 March 2013 
Gravothermal catastrophe: The dynamical stability of a fluid model^{⋆}
^{1} Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy
email: mattia.sormani@physics.ox.ac.uk
^{2} Università degli Studi di Milano, Dipartimento di Fisica, via Celoria 16, 20133 Milano, Italy
email: giuseppe.bertin@unimi.it
Received: 30 October 2012
Accepted: 21 January 2013
A reinvestigation of the gravothermal catastrophe is presented. By means of a linear perturbation analysis, we study the dynamical stability of a spherical selfgravitating isothermal fluid of finite volume and find that the conditions for the onset of the gravothermal catastrophe, under different external conditions, coincide with those obtained from thermodynamical arguments. This suggests that the gravothermal catastrophe may reduce to Jeans instability, rediscovered in an inhomogeneous framework. We find normal modes and frequencies for the fluid system and show that instability develops on the dynamical time scale. We then discuss several related issues. In particular, (1) for perturbations at constant total energy and constant volume, we introduce a simple heuristic term in the energy budget to mimic the role of binaries. (2) We outline the analysis of the twocomponent case and show how linear perturbation analysis can also be carried out in this more complex context in a relatively straightforward way. (3) We compare the behavior of the fluid model with that of the collisionless sphere. In the collisionless case the instability seems to disappear, which is at variance with the linear Jeans stability analysis in the homogeneous case. We argue that a key ingredient for understanding the difference lies in the role of the detailed angular momentum in a collisionless system. A spherical stellar system is expected to undergo the gravothermal catastrophe only in the presence of some collisionality, which suggests that the instability is dissipative and not dynamical. Finally, we briefly comment on the meaning of the Boltzmann entropy and its applicability to the study of the dynamics of selfgravitating inhomogeneous gaseous systems.
Key words: hydrodynamics / instabilities / galaxies: clusters: general / gravitation
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2013
1. Introduction
Core collapse in an Nbody 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 Nbody 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 Nbody problem reach equilibrium? How does the equilibrium of a model relate to the equilibrium of the pure Nbody 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 Nbody problem (in particular gaseous models, collisionless models, and fluid models) admit equilibrium configurations that are spatially truncated selfgravitating isothermal spheres, but it is not clear to what extent they can be considered as representative of the equilibrium states of the pure Nbody problem. Several studies have addressed the stability problem of isothermal spheres using thermodynamical methods, starting with Antonov (1962) and LyndenBell & 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: (1)where f is the oneparticle distribution function^{1}, r and v the position and velocity vector respectively, and k the Boltzmann constant.
Unfortunately, the thermodynamics of selfgravitating systems still depends on a number of unresolved issues, partly because of the longrange 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 selfgravitating 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 LyndenBell & 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 SmoluchowskiPoisson system (different from the EulerPoisson 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 LyndenBell & 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 twocomponent case, by perturbing the twocomponent 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 twocomponent 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 sphericallysymmetric 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 twocomponent 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 Nbody problem. In Sect. 6 we draw our conclusions and identify some open questions and issues.
2. Thermodynamical approach and Boltzmann entropy
Antonov (1962) and LyndenBell & Wood (1968) based their stability analysis of selfgravitating spatiallytruncated 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^{2} total mass M, total energy E, and volume V = 4πR^{3}/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 singleparticle mass m, each spatially truncated isothermal sphere is identified by two dimensional scales (e.g., the central density ρ_{0}(0) and the temperature T) and one dimensionless parameter. A possible choice for the dimensionless parameter is Ξ = R/λ, where λ = [kT/m4πGρ_{0}(0)] ^{1/2} is reminiscent of the Jeans length (Jeans 1902); another equivalent choice is the density contrast ρ_{0}(0)/ρ_{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 onetoone 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}(0)/ρ_{0}(R) < 709) the equilibrium configurations are local maxima, whereas they are saddle points for Ξ > 34.4. Therefore, he concluded that selfgravitating isothermal spheres are unstable if Ξ > 34.4.
The work of Antonov (1962) was extended by LyndenBell & 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 selfgravitating 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 nonnegligible amount. The correct expression for the entropy (corresponding to phenomenological thermodynamics) of systems with a welldefined 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 selfgravitating 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 selfgravitating 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 LyndenBell & 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 twobody 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 oneparticle distribution function to be a Maxwellian (i.e., of the form f(r,v) = Aexp [ − mv^{2}/2kT(r)], where A is a normalization constant) and by defining the pressure as p ≡ ∫mf(v^{2}/3) d^{3}r d^{3}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 Nbody system when t → ∞ (t is time) has an effective equation of state^{3} 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 selfgravitating Nbody system, since it assumes the equation of state of an ideal gas. Moreover, because gravitational forces are longranged, each particle feels the influence of all other distant particles. Thus, it may be that, strictly speaking, in the case of selfgravitating 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 socalled mean field approximation (Katz 2003; Padmanabhan 1990). However, the range of applicability of this approximation is still not clear, also because a smallscale cutoff is needed to avoid divergences in the microcanonical entropy (Padmanabhan 1990). The Boltzmann entropy is also the H quantity of the Boltzmann Htheorem; 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 selfgravitating isothermal spheres are not true thermal equilibrium states of the pure Nbody 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 selfgravitating fluid governed by the NavierStokes and continuity equations, together with the equation of state of an ideal gas (for example, see Landau & Lifshitz 1987): (2)where ρ is the density, u the fluid velocity, p the local pressure, Φ the gravitational potential, η and ζ are viscosity coefficients, m is the oneparticle 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: (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 selfgravitating truncated isothermal sphere is identified by two dimensional scales (for example the central density ρ_{0}(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}(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), (4)where ρ_{0}(0) is the unperturbed central density and ψ(ξ) is the regular solution to the Emden equation (the symbol ′ denotes derivative with respect to the argument ξ): (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_{0} that the fluid element under consideration had at the initial time t_{0}.
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: (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: (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: (8)where f(ξ) = ρ_{0}(ξ)u_{1}(ξ) is the unknown function, (9)represents the dimensionless (squared) eigenfrequency and ℒ is the following differential operator: (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, (11)To solve the equations, we still need to specify the function T_{1}(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): (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_{0} of the fluid element under consideration at the initial time t_{0}. Thus in Eq. (2)we have to perform the following change of independent variables: (13)In the Lagrangian representation, each quantity is a function of the new independent variables r_{0} 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):
As for the Eulerian case we separate each quantity in an unperturbed part and a perturbation: (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: (16)In the following we shall drop the symbol ~ for simplicity of notation. After introducing the dimensionless radius ξ_{0} ≡ r_{0}/λ, using the density profile (4)of the unperturbed state and after some manipulations we obtain the following equation for ρ_{1}(ξ_{0}) (see Appendix B.2.2 for an outline of the calculations): (17)where L = ω^{2}/4πGρ_{0}(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_{1} = 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 selfgravitating 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_{1} = 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 (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 (19)where is an arbitrary point in the domain of f, Eq. (18)can be recast in the form of a Schrödinger equation: (20)where U(ξ) is the effective potential given by (21)The boundary conditions become (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 LyndenBell & 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 (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 (see also Semelin et al. 2001; Chavanis 2002).
4.1.2. Constant {E, V} case (gravothermal catastrophe)
Here we consider a selfgravitating fluid sphere at constant total energy E and volume V. In this case the instability has been named gravothermal catastrophe by LyndenBell & Wood (1968). The total energy E of the fluid is defined as (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 SmoluchowskiPoisson 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 , we have (25)By substituting T = T_{0} + T_{1}, ρ = ρ_{0} + ρ_{1} in Eq. (25), only keeping firstorder quantities and imposing that the energy remains constant, we obtain the following expression for T_{1}: (26)Here Φ_{0} is the gravitational potential of the unperturbed density distribution ρ_{0}. After some manipulations (see Appendix B.3) we obtain (27)By substituting Eq. (27)in Eq. (8), recalling the definition λ = [kT_{0}/4πGρ_{0}(0)m]^{1/2}, and neglecting viscosity, we obtain (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 (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}(0)/ρ_{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: (30)where (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 (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: (33)Using the expression of C (31)and dividing by b, we obtain (34)A second condition on a/b follows from the boundary condition f(Ξ) = 0: (35)Substituting a/b from (35)in (34)we obtain (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 “corehalo” structure, that is, an oscillation in ρ_{1}: the density perturbation is positive in the inner part (nucleus), negative in the middle (emptying area), and then positive again (halo). The corehalo structure has been physically interpreted in the framework of the gravothermal catastrophe given by LyndenBell & 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 corehalo 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: (37)which contains an additional term proportional to the dimensionless parameter υ and to the central density ρ(0). The quantity 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_{1} (27)obtained previously, to be substituted in Eq. (8), is replaced with the expression for T_{1} that is obtained from Eq. (37)by substituting T = T_{0} + T_{1}, ρ = ρ_{0} + ρ_{1} 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 Nbody 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 LyndenBell & Eggleton (1980) a phenomenological energygeneration 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 Nbody simulations (Makino 1996). We comment further on this topic in Sect. 5.2.
4.1.3. Constant {T, V}: twocomponent case
We briefly studied the problem of dynamical stability by means of a linear modal analysis of a twocomponent ideal fluid. Each component is assumed to interact with the other only through the common gravitational potential. The hydrostatic equilibrium configurations are spatially truncated twocomponent 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 singleparticle masses m_{A} and m_{B}, with m_{B} > m_{A}. Twocomponent isothermal spheres are characterized by two additional dimensionless parameters (with respect to the onecomponent case): the ratio of the singleparticle 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}/ρ_{0} and ρ_{1B}/ρ_{0} 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}/ρ_{0}, ρ_{1B}/ρ_{0} 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 SmoluchowskiPoisson system of equations. Breen & Heggie (2012a,b) found that the heavier component dominates gravothermal oscillations in twocomponent 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_{1} = 0 and Eq. (17)becomes: (38)Boundary conditions are as follows. The condition u_{1}(0) = 0 translates into the condition (using Euler Eq. (14), under the assumption that ∂u/∂r_{0} 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: (39)From mass continuity and by imposing (which is true in the Eulerian representation and thus we expect to be true in the present case), we also have the condition ρ_{1}(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 LyndenBell & 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 nonsphericallysymmetric 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 (40)where ℒ is defined as (41)From Emden Eq. (5), the operator ℒ has the following properties: (42)Hence, an analytical solution of Eq. (40)is the following (for a derivation in the Eulerian representation, see also Chavanis 2003): (43)This solution satisfies the boundary conditions ρ_{1}(0) = 2 and , as required by Eq. (39). The zeros of G_{TP}(ξ_{0}) allow us to identify the marginally stable normal modes that satisfy the correct boundary conditions (39). The first zero of G_{TP}(ξ_{0}) 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 selfgravitating 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: (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 (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: (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 selfgravitating isothermal sphere
Isothermal spheres are equilibrium configurations for several idealized models of the pure Nbody problem. In particular, selfgravitating isothermal spheres can be studied as stationary states of the collisionless Boltzmann equation: (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^{4}. 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^{5} (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 selfgravitating 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 selfgravitating 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^{6}. 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 Nbody 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 Nbody 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 (48)Therefore, for a stellar system of N stars and total mass M we can define a radius related to detailed angular momentum as (49)where A is the sum of the squares of the angular momenta of the individual stars, that is, the following quantity: (50)Thus we may argue that the typical time scale for core collapse should correlate with (51)The main mechanism through which A varies with time is expected to be that of twobody collisions. Thus t_{cc} should be similar to the twobody relaxation time.
By monitoring the quantity A(t) in Nbody 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 twobody relaxation time.), its relevance to the studies of core collapse in the gaseous model (LyndenBell & Eggleton 1980; Sugimoto & Bettwieser 1983), and its connection with the phenomenon of the gyrogravothermal catastrophe (Hachisu 1979).
6. Discussion and conclusions
In this paper we have studied systematically the dynamical stability of a selfgravitating 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; LyndenBell & 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 twocomponent 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 selfgravitating 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 Nbody 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 particleparticle interaction, for nonNewtonian cases (Padmanabhan 1989). Potentials that exhibit a softening at small radii, such as , where r_{0} 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 selfgravitating fluid isothermal spheres.
The equations for the hydrostatic equilibrium of a spherically symmetric fluid with the equation of state of an ideal gas are 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πρ_{0}(r)r^{2}. By making the change of variable ρ_{0}(r) = ρ_{0}(0)e^{− ψ(r)}, where ρ_{0}(0) is the central density, and introducing the dimensionless radius ξ = r/λ, where λ = [kT/4πGρ_{0}(0)m]^{1/2}, we obtain the differential equation for ψ(ξ), recorded in the main text as Eq. (5). (Because the constant ρ_{0}(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ξ^{2} and ψ′(ξ) ~ 2/ξ as ξ → ∞.
From Eq. the mass enclosed within the radius ξ is (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}(0) and one dimensionless parameter Ξ. The density profile of a spatially truncated isothermal sphere is given by ρ_{0}(r) = ρ_{0}(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 Then we look for solutions of the form (7). From Eq. (B.1)we obtain ρ_{1}: (B.3)and thus eliminate it from Eq. (B.2)to find the radial component of the NavierStokes equation: (B.4)where u_{1} is the radial component of the velocity. By defining f(r) ≡ ρ_{0}(r)u_{1}(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:
Now we perform a change of variables. In the Lagrangian representation, each quantity is described as a function of the new independent variables r_{0} and t, where r_{0} is the position of the fluid element at t = t_{0}. The standard rules to transform derivatives of a generic function f are: (B.6)The partial derivatives of r_{0} are obtained from the following relation, which expresses the condition that two fluid shells do not cross each other: (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 (B.8)By taking the partial derivative of Eq. (B.7)with respect to t we obtain (B.9)From the continuity equation, Eq. (B.9)then becomes (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_{0},t) = M(r_{0},t = t_{0}), we obtain (B.11)From the usual rules of derivation, we have: (B.12)By applying Eqs. (B.10)and (B.8), Eq. (B.12)can be written as (B.13)Now we wish to obtain one equation involving only the unknown ρ_{1}, by combining the three Eqs. (B.11)and (B.13). We assume the modal dependence (16). We eliminate r_{1} from (B.13)and the second of (B.11). We then eliminate u_{1} from the resulting equation and the first of (B.11), to obtain the following equation {we also used hydrostatic equilibrium to replace }: (B.14)By referring to the dimensionless radius ξ_{0} ≡ r_{0}/λ, 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 Φ_{0} of the unperturbed density distribution in the following way: (B.15)From Eq. (B.15), by recalling that we only consider perturbations that do not change the total mass (), we obtain (B.16)From Eq. (B.3), by setting f = ρ_{0}u_{1} 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 singleparticle Lagrangian as (C.1)where V(r,t) is a timedependent 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 (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^{2}/r^{3} − GmM(r,t)/r^{2}. The following inequality holds: (C.3)By multiplying by ṙ and integrating both sides of Eq. (C.2), we obtain (C.4)Since ṙ^{2}(t)/2 is positive, the righthand side of Eq. (C.5)must be positive. By taking r(t) ≤ r(t_{0}) (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 (C.5)For given values of r(t_{0}) and ṙ(t_{0}), 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; ρ_{1}/ρ_{0} and u_{1} 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 LyndenBell & 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 corehalo 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 corehalo 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 twocomponent case
In this appendix we summarize the equations of the linear analysis for the twocomponent 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 twocomponent selfgravitating 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 (E.1)\vspace*{1.2mm}where ρ_{A0} and ρ_{B0} are respectively the density profiles of the lighter and heavier component, ξ ≡ r/λ_{2} is the dimensionless radial coordinate, where λ_{2} ≡ [kT(1/m_{A} + 1/m_{B})/4πGρ_{0}(0)]^{1/2} and we denote by ρ_{0}(ξ) ≡ ρ_{A0}(ξ) + ρ_{B0}(ξ) the total unperturbed density; Ξ is the value of ξ at the truncation radius, β ≡ m_{B}/m_{A} is the ratio of the singleparticle masses, ψ is the solution of the following generalization of the Emden Eq. (5): 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 twocomponent 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: (E.4)Here L = ω^{2}/4πGρ_{0}(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: (E.5)Similarly to the onecomponent 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 twocomponent 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 β = m_{B}/m_{A} = 3 at different values of the total mass ratio M_{B}/M_{A}. Even for low values of M_{B}/M_{A}, the density perturbation of the heavier component dominates, thus suggesting that the heavier component is the more important driver of the instability. 
Acknowledgments
We would like to thank Marco Lombardi, Francesco Pegoraro, and Steven N. Shore for many interesting comments and discussions.
References
 Antonov, V. A. 1962, in Vest. Leningrad Univ. 7 135 (in Russian). English translation in Dynamics of star clusters, eds. J. Goodman, & P. Hut (Dordrecht: Reidel), IAU Symp. 113, 525 [Google Scholar]
 Batt, J., Morrison, P. J., & Rein, G. 1995, Archive for Rational Mechanics and Analysis, 130, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Bertin, G. 2000, Dynamics of Galaxies (Cambridge University Press) [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
 Bonnor, W. B. 1956, MNRAS, 116, 351 [CrossRef] [Google Scholar]
 Bouchet, F., Gupta, S., & Mukamel, D. 2010, Phys. A Stat. Mech. Appl., 389, 4389 [Google Scholar]
 Breen, P. G., & Heggie, D. C. 2012a, MNRAS, 425, 2493 [NASA ADS] [CrossRef] [Google Scholar]
 Breen, P. G., & Heggie, D. C. 2012b, MNRAS, 420, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Campa, A., Dauxois, T., & Ruffo, S. 2009, Phys. Rep., 480, 57 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Casetti, L., & Nardini, C. 2012, Phys. Rev. E, 85, 061105 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1967, An introduction to the study of stellar structure (New York: Dover) [Google Scholar]
 Chavanis, P. H. 2002, A&A, 381, 340 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P. H. 2003, A&A, 401, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P. H. 2006, Int. J. Mod. Phys. B, 20, 3113 [Google Scholar]
 Chavanis, P. H., & Ispolatov, I. 2002, Phys. Rev. E, 66, 036109 [Google Scholar]
 Chavanis, P.H., Rosier, C., & Sire, C. 2002, Phys. Rev. E, 66, 036105 [Google Scholar]
 de Vega, H. J., & Siebert, J. A. 2002, Phys. Rev. E, 66, 016112 [NASA ADS] [CrossRef] [Google Scholar]
 Ebert, R. 1955, Z. Astrophys., 37, 217 [Google Scholar]
 Hachisu, I. 1979, PASJ, 31, 523 [Google Scholar]
 Hachisu, I., & Sugimoto, D. 1978, Progr. Theor. Phys., 60, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D., & Hut, P. 2003, The Gravitational MillionBody Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge University Press) [Google Scholar]
 Inagaki, S. 1980, PASJ, 32, 213 [NASA ADS] [Google Scholar]
 Ispolatov, I., & Cohen, E. G. D. 2001, Phys. Rev. E, 64, 056103 [Google Scholar]
 Jaynes, E. T. 1965, Am. J. Phys., 33, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Jeans, J. H. 1902, Roy. Soc. London Philosoph. Trans. Ser. A, 199, 1 [Google Scholar]
 Kandrup, H. E. 1990, ApJ, 351, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Kandrup, H. E., & Sygnet, J. F. 1985, ApJ, 298, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Katz, J. 1978, MNRAS, 183, 765 [NASA ADS] [Google Scholar]
 Katz, J. 2003, Found. Phys., 33, 223 [CrossRef] [MathSciNet] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1987, Fluid Mech. (Course of Theoretical Physics), 2nd edn. (ButterworthHeinemann), 6 [Google Scholar]
 Lightman, A. P. 1977, ApJ, 215, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Lombardi, M., & Bertin, G. 2001, A&A, 375, 1091 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LyndenBell, D., & Eggleton, P. P. 1980, MNRAS, 191, 483 [NASA ADS] [Google Scholar]
 LyndenBell, D., & Wood, R. 1968, MNRAS, 138, 495 [NASA ADS] [Google Scholar]
 Makino, J. 1996, ApJ, 471, 796 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, R. H. 1973, ApJ, 180, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Mukamel, D. 2008, in Dynamics and Thermodynamics of Systems with Long Range Interactions: Theory and Experiments, eds. A. Campa, A. Giansanti, G. Morigi, & F. S. Labini, ASP Conf. Ser., 970, 22 [Google Scholar]
 Nakada, Y. 1978, PASJ, 30, 57 [NASA ADS] [Google Scholar]
 Padmanabhan, T. 1989, ApJS, 71, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Padmanabhan, T. 1990, Phys. Rep., 188, 285 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Semelin, B., Sánchez, N., & de Vega, H. J. 2001, Phys. Rev. D, 63, 084005 [NASA ADS] [CrossRef] [Google Scholar]
 Sopik, J., Sire, C., & Chavanis, P.H. 2005, Phys. Rev. E, 72, 026105 [NASA ADS] [CrossRef] [Google Scholar]
 Sugimoto, D., & Bettwieser, E. 1983, MNRAS, 204, 19P [NASA ADS] [Google Scholar]
 Taff, L. G., Hansen, C. J., Ross, R. R., & van Horn, H. M. 1975, ApJ, 197, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Thirring, W. 1970, Z. Phys., 235, 339 [Google Scholar]
 Yabushita, S. 1968, MNRAS, 140, 109 [NASA ADS] [Google Scholar]
 Yoshizawa, M., Inagaki, S., Nishida, M. T., et al. 1978, PASJ, 30, 279 [NASA ADS] [Google Scholar]
All Figures
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. 

In the text 
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). 

In the text 
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 text 
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 the text 
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 LyndenBell & Wood (1968). 

In the text 
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 corehalo structure described in Sect. 4.1.2 disappears between L = 0.021 and L = 0.022. 

In the text 
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 corehalo structure described in Sect. 4.1.2 disappears when the velocity has no internal zeros, between L = 0.021 and L = 0.022. 

In the text 
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. 

In the text 
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 twocomponent 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 β = m_{B}/m_{A} = 3 at different values of the total mass ratio M_{B}/M_{A}. Even for low values of M_{B}/M_{A}, the density perturbation of the heavier component dominates, thus suggesting that the heavier component is the more important driver of the instability. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.