4 Dynamical instability of isothermal gas spheres in general relativity

   
4.1 The condition of marginal stability

The stability of relativistic stars against radial perturbations is usually studied with the equation of pulsations derived by Chandrasekhar (1964). Then, the stability of stars can be reduced to a Sturm-Liouville problem, which is usually solved numerically (see, e.g., Bardeen et al. 1966). The method is quite general but it is not the most efficient in the present situation. Indeed, by formulating the problem differently, using the form of the equation of pulsations given by Yabushita (1973) and making use of the Milne variables, it is possible to study the stability of bounded isothermal spheres analytically, like in the Newtonian case.

Let us consider a small perturbation $\delta\epsilon$ around a configuration of isothermal equilibrium and let us introduce the mass perturbation $f(\xi)$ within the sphere of radius $\xi$ such that

$${\delta\epsilon\over\epsilon_{0}}={1\over 4\pi\xi^{2}}{{\rm d}f\over {\rm d}\xi}\cdot$$ (163)

Yabushita (1973, 1974) has shown that the equation for radial pulsations satisfied by the function f can be written

 

$${{\rm d}^{2}f\over {\rm d}\xi^{2}}+\biggl (-{2\over\xi}+{{\rm d}\... ...\rbrack f={\sigma^{2}c^{4}\over 4\pi G\epsilon_{0}(1+q)}{\rm e}^{\lambda-\nu}f,$$ (164)

where $\sigma$ is the period of oscillations defined by $\delta\epsilon\sim {\rm e}^{\sigma ct}$. On the other hand, the velocity profile of the perturbation is given by

$$\delta u^{1}=-{\sigma c\over q+1}{1\over 4\pi \epsilon}{f\over r^{2}}{\rm e}^{-\nu/2}.$$ (165)

If the system is confined within a box, we must require that $\delta u^{1}=0$ at r=R. The function $f(\xi)$ must therefore satisfy the boundary conditions

$$f(0)=f(\alpha)=0.$$ (166)

These boundary conditions are equivalent to the conservation of the mass-energy M. In the Newtonian limit $q\rightarrow 0$, the equation of pulsations (164) reduces to the one considered by Yabushita (1968) and Chavanis (2001). Clearly, the condition of instability corresponds to $\sigma^{2}>0$. We shall be particularly interested by the point of marginal stability $(\sigma=0)$ in the series of equilibrium at which the system becomes unstable. We then have to solve the equation

 

$${{\rm d}^{2}F\over {\rm d}\xi^{2}}+\biggl (-{2\over\xi}+{{\rm d}\... ...biggl \lbrack 1+{2q\xi\over 1+q}{{\rm d}\psi\over {\rm d}\xi}\biggr\rbrack F=0,$$ (167)

with $F(0)=F(\alpha)=0$. If we denote by ${\cal L}$ the differential operator that appears in Eq. (167), it is possible to show that (Yabushita 1974)

 

$${\cal L}(M(\xi))={\cal L}(\xi^{3}{\rm e}^{-\psi})={{\rm d}\psi\ov... ...}+{2q\xi\over 1+q}{\rm e}^{\lambda-\psi}M(\xi){{\rm d}\psi\over{\rm d}\xi}\cdot$$ (168)

Therefore, the eigenfunction associated with $\sigma=0$ is

$$F(\xi)=c_{1}(\xi^{3}{\rm e}^{-\psi}-M(\xi)).$$ (169)

In the classical limit, we recover the result $\xi^{3}{\rm e}^{-\psi}-\xi^{2}\psi'$ derived from the analysis of the second order variations of the free energy (in the canonical ensemble) or from the Navier-Stokes equations (Chavanis 2001). The condition F(0)=0 is automatically satisfied. The condition $F(\alpha)=0$determines the values of $\alpha$ at which a new mode of stability is lost. According to Eq. (169), we get

$$M(\alpha)=\alpha^{3}{\rm e}^{-\psi(\alpha)}.$$ (170)

Using Eq. (144) and introducing the Milne variables (139), it is straightforward to check that this relation is equivalent to

$$pv_{0}={1\over u_{0}}-q-1,$$ (171)

which is precisely Eq. (148). Therefore, a new eigenvalue $\sigma^{2}$ becomes positive at the points where $\chi$ is extremum. In particular, the first instability sets in precisely at the point of maximum mass. Isothermal configurations with a central density $\alpha>\alpha_{1}$, a density contrast ${\cal R}>{\cal R}_{\rm c}$ or a central redshift $z_{0}>z_{0}^{\rm c}$ are unstable with respect to radial perturbations. Similar results were obtained by Misner & Zapolsky (1964) who numerically solved the Chandrasekhar equation for radial pulsations applied to their neutron star model. Their method was, however, limited to the fundamental mode and they could not demonstrate that secondary mass peaks correspond to the change of sign for the next eigenvalues (see their discussion). Our method gives an unambiguous definite answer to that question, at least in the framework of our box model. It has to be noted that these results differ from those obtained by Yabushita (1974) who found that instability occurs before the first mass peak (except in the Newtonian case q=0). Now, the three models only differ in the boundary conditions: (i) Yabushita's gas spheres have envelopes which exert a constant pressure. (ii) In our model, the volume of the isothermal region is assumed to be fixed. (iii) Misner & Zapolsky introduce a more realistic envelope in which the equation of state has a physically admissible form. Our results suggest that imposing a fixed volume (Antonov point of view) is a better approximation than imposing a fixed pressure (Bonnor point of view) since it provides the same condition of stability as the reference model of Misner & Zapolsky (1964) and other works.

   
4.2 Perturbation profiles and secondary instabilities

The values of $\alpha$ at which a new mode of instability occurs are determined by Eq. (171). For large values of $\alpha$ (high order modes), we can introduce the asymptotic expansion (134) in the Milne variables. Then, after some algebra, we find that the values of $\alpha_{n}$ at which a new instability occurs are asymptotically given by

 

$${(7+42q-q^{2})^{1/2}\over 2(1+q)}\ln\alpha_{n}+\delta =\arctan\biggl ({q^{2}+20q+3\over (1+q)(7+42q-q^{2})^{1/2}}\biggr )+n\pi.$$ (172)

In continuity with the Newtonian limit, they follow the geometric progression

$$\alpha_{n}\sim \biggl\lbrack {\rm e}^{2\pi (1+q)\over (7+42q-... ...\biggr\rbrack^{n}\qquad (n\rightarrow +\infty, {\rm integer}).$$ (173)

For q=0,1/3,1 the ratio in the geometric progression is respectively 10.74, 6.25 and 6.13. It is minimum for q=7/11 with a value 5.93.

The profile of the energy density perturbation associated with each instability is obtained by substituting Eq. (169) in Eq. (163). Simplifying the derivative with the aid of Eq. (125) and introducing the Milne variable v, we obtain

$${\delta\epsilon\over\epsilon}={c_{1}\over 4\pi}(2-v).$$ (174)

The zeros of the profile $\delta\epsilon$ are determined by the intersections between the spiral in the (u,v) plane and the line v=2 (see Fig. 11). The first two zeros are represented in Fig. 16 for different values of q. We see that $\xi_{2}$ is always larger than $\alpha _{1}$ so it does not satisfy the requirement $0\le \xi\le \alpha$. Therefore, the first mode of instability has only one node at $\xi=\xi_{1}$. In other words, the perturbation profile that triggers the instability of an isothermal sphere in general relativity does not present a "core-halo'' structure (see Fig. 17) in continuity with the Newtonian study in the canonical ensemble (Chavanis 2001).

Figure 16: Critical parameter $\alpha _{1}$ as a function of q (full line) together with the first two zeros of the perturbation profile (174) satisfying $v(\xi _{i})=2$ (dashed lines). Since $\xi _{2}>\alpha _{1}$, this solution must be rejected.

Figure 17: Density perturbation profiles at the point of marginal stability for different values of q. The perturbation profiles do not present a core-halo structure.

The high order modes of instability have several nodes. Their positions (for large $\xi_{i}$) are obtained by expanding the Milne variable v with the aid of Eq. (134) and substituting the resulting expression in the condition v=2. This yields

 

$${(7+42q-q^{2})^{1/2}\over 2(1+q)}\ln\xi_{i}+\delta =-\arctan\biggl ({1+3q\over (7+42q-q^{2})^{1/2}}\biggr )+i\pi.$$ (175)

Therefore, the zeros of the density perturbation profiles also follow a geometric progression

$$\xi_{i}\sim \biggl\lbrack {\rm e}^{2\pi (1+q)\over (7+42q-q^{... ...\biggr\rbrack^{i}\qquad (i\rightarrow +\infty, {\rm integer}),$$ (176)

with the same ratio as the dimensionless radii $\alpha_{n}$.

We can also determine the velocity profile (165) associated with each mode of instability. Introducing the velocity of sound $c_{\rm s}=q^{1/2}c$ and expressing the result in terms of the dimensionless variables defined in Sect. 3.3, we get

$${\delta u^{1}\over c_{\rm s}}=-{c\sigma' \over 4\pi(q+1)\xi^{2}}{\rm e}^{\psi(\xi)\over 1+q}F(\xi),$$ (177)

where $F(\xi)$ is given by Eq. (169) at the critical points. In writing Eq. (177), it is implicitly understood that we are just at the onset of instability ( $\sigma'=0^{+}$), so that Eq. (177) is applicable with $\sigma'>0$ (the velocity profile at the critical points, i.e. $\sigma'=0$, is simply $\delta u^{1}=0$). In terms of the Milne variables, the velocity profile (177) can be written

 

$${\delta u^{1}\over c_{\rm s}}=-{c\sigma' \over 4\pi(q+1)}c_{1}\psi'{\rm e}^{\psi(\xi)\over 1+q}{u\over 1+pv}\biggl (pv-{1\over u}+q+1\biggr ).$$ (178)

For q=0, we recover the expression obtained in the Newtonian context from the Navier-Stokes equations (Chavanis 2001). Since $F(0)=F(\alpha)=0$, the velocity always vanishes at the center and at the boundary of the domain. The other zeros, for which $F(\xi_{i})=0$, correspond precisely to the values of $\alpha$ at which a new mode of instability occurs, i.e. $\xi_{i}=\alpha_{i}$. The nth mode of instability therefore has n+1 zeros $\xi=0$, $\xi=\alpha_{1}$, $\xi=\alpha_{2}$,..., $\xi=\alpha_{n}$. For high order modes of instability, the zeros follow asymptotically the geometric progression (173).

   
4.3 The baryon number: Another interpretation of the instability

It has often been suggested on the basis of phenomenological arguments that the baryon number N, or equivalently the rest mass M₀=Nm, is also maximum at the onset of instability. We wish to check this prediction explicitly in the case of an isothermal gas sphere.

The baryon number is obtained by multiplying the baryon number density n by the proper volume element ${\rm e}^{\lambda(r)/2}4\pi r^{2}{\rm d}r$ and integrating the resulting expression over the whole configuration. Using Eq. (106), we get (see, e.g., Weinberg 1972)

$$N=\int_{0}^{R}n\biggl\lbrack 1-{2GM(r)\over rc^{2}}\biggr\rbrack^{-1/2}4\pi r^{2}{\rm d}r.$$ (179)

Now, according to the equations of state (116)-(118), the baryon number density is related to the energy density by

$$n^{\gamma}={q\over K}\epsilon.$$ (180)

Introducing the dimensionless variables defined in Sect. 3.3, we find that the baryon number can be written

$${N\over N_{*}}={1\over\alpha^{3q+1\over 1+q}}\int_{0}^{\alpha... ...ck 1-p{M(\xi)\over \xi}\biggr \rbrack^{-1/2}\xi^{2}{\rm d}\xi,$$ (181)

with the normalization constant

$$N_{*}=4\pi R^{3}\biggl \lbrack {q^{2}c^{4}\over 4\pi GKR^{2}(1+q)}\biggr\rbrack^{1/\gamma}\cdot$$ (182)

Figure 18: Baryon number and total mass-energy as a function of the central density (through the parameter $\alpha$) for q=1/3. The peaks occur for the same values of the central density.

Figure 19: Total mass-energy versus baryon number for isothermal gaseous spheres (q=1/3).

The baryon number (181) is ploted versus the central density in Fig. 18 for q=1/3. It exhibits damped oscillations, like the total mass energy, and tends asymptotically to the value ${1+q\over 3q+1}Q^{1\over 1+q}(1-pQ)^{-1/2}$ as $\alpha\rightarrow +\infty$ (singular sphere). It is observed numerically that the peaks of mass-energy and of baryon number occur for the same values of the central density (we have checked this property for various values of q). This is also clear from Fig. 19 where the mass-energy M is ploted as a function of the baryon number N. From this figure, we can understand physically why the system becomes unstable after the first mass peak. Indeed, for a given rest mass Nmc² there exists several solutions with different energy Mc². Only the solution with the lowest energy ("fundamental'' state) is stable with respect to radial perturbations (see Weinberg 1972). We expect therefore some phase transitions to occur from a state of high energy to a state of lower energy (see the arrow in Fig. 19). Therefore, only the lower branch is stable and it corresponds to the configurations prior to the first mass peak (we come to the same conclusion by minimizing the binding energy E=(M-M₀)c² at fixed mass M). These results are consistent with those found by Tooper (1964) for polytropic gas spheres.

   
4.4 Characteristic length scales

We shall now reformulate the previous results in a different form by introducing appropriate length scales. This will provide a direct comparison with our discussion in the Newtonian context (Chavanis 2001). In general relativity, a gaseous sphere of mass M is unstable if its radius is smaller than a multiple of the Schwarzschild radius (Chandrasekhar 1964). For a confined isothermal gas, this condition is explicitly given by Eq. (150) which can be rewritten

$$R\le {1\over \chi_{\rm c}(q)}R_{\rm S},\qquad R_{\rm S}={2GM\over c^{2}}\cdot$$ (183)

On the other hand, in Newtonian gravity, the condition of gravitational collapse is that the size of the domain be larger than the Jeans length $L_{\rm J}$. For example, in the case of an isothermal gas confined within a box, the condition of instability reads (de Vega & Sanchez 2001; Chavanis 2001)

$$R\ge \sqrt{2.52...\over 3}L_{\rm J},\qquad L_{\rm J}=\biggl ({3R^{3}\over \beta GMm}\biggr )^{1/2}\cdot$$ (184)

It should be emphasized that these inequalities characterize the absence of hydrostatic equilibrium for an isothermal gas beyond these thresholds. This is more stringent than the instability of an equilibrium configuration (see below). Although the sense of these inequalities is opposite, there is no contradiction since the size of the domain enters explicitly in the definition of the Jeans length (through the average density). For $q\rightarrow 0$, the two criteria become equivalent as we now show. In the classical limit, the rest mass dominates the kinetic energy so that $\epsilon=\rho c^{2}$ where $\rho=n m$ is the classical mass density. Therefore, the equation of state (118) becomes $p=q\rho c^{2}$ and it must be compared with the equation of state of a classical isothermal gas $p={\rho\over m}kT$. This yields

$$q={kT\over mc^{2}}={1\over\beta mc^{2}}\qquad (\beta mc^{2}\rightarrow +\infty).$$ (185)

Now, in the classical limit $c\rightarrow +\infty$, the critical parameter $1/\chi_{\rm c}$, given by Eq. (149), diverges like $1/\chi_{\rm c}=1/2qv(\alpha_{1})\rightarrow +\infty$ while the Schwarzschild radius $R_{\rm S}$ goes to zero in such a way that the product $(1/\chi_{\rm c})R_{\rm S}$ remains finite. Therefore, the inequality (183) is reversed

$$R\le {\beta G M m\over v(\alpha_{1})}={3R^{3}\over v(\alpha_{... ...2}} \Longrightarrow R\ge \sqrt{v(\alpha_{1})\over 3}L_{\rm J},$$ (186)

and we recover the result (184) with $v(\alpha_{1})=2.52$. It can be noted that the same criterion can be written

$$M\ge M_{\rm c}\equiv 2.52 {RkT\over Gm}\cdot$$ (187)

Therefore, a classical gaseous sphere of radius R and temperature T can exist only below a limiting mass $M_{\rm c}$ given by formula (187). Obviously, the physical origin of this limiting mass is the same as for the limiting mass of neutron stars but with a completely different interpretation of the parameters.

It is also relevant to introduce another length scale

$$L_{\rm c}=\biggl ({9qc^{4}\over 4\pi G\epsilon_{0}(1+q)}\biggr )^{1/2},$$ (188)

which gives a good estimate of the "core radius'' of a relativistic isothermal gas sphere. In the classical limit, $L_{\rm c}$ reduces to the King's length (Binney & Tremaine 1987). In terms of the core radius, the parameter $\alpha$ can be written

$$\alpha=3{R\over L_{\rm c}}\cdot$$ (189)

The condition of instability $\alpha>\alpha_{1}$ is therefore equivalent to

$$R>{\alpha_{1}\over 3}L_{\rm c}.$$ (190)

Above this threshold, it is still possible to construct isothermal gas spheres in hydrostatic equilibrium but they are unstable. Secondary instabilities occur for radii that follow a geometric progression

$$R_{n}\sim\biggl\lbrack {\rm e}^{2\pi (1+q)\over (7+42q-q^{2})... ...ack^{n} L_{\rm c}\qquad (n\rightarrow +\infty, {\rm integer}).$$ (191)

In the Newtonian context, this hierarchy of scales was interpreted by Semelin et al. (1999) and Chavanis (2001) in relation with the fragmentation of the interstellar medium and the possible origin of its fractal structure.

The instability criteria (183) and (190), expressed in terms of distinct length scales, can also be expressed in terms of the redshift, with a similar distinction between the central redshift and the surface redshift. Indeed, the surface redshift z(R) is an increasing function of the mass energy $\chi$ while the central redshift z₀ is an increasing function of the central density $\alpha$. Therefore, the surface redshift must necessarily satisfy $z(R)<z(R)^{\rm c}$ since no hydrostatic equilibrium of isothermal spheres can exist above this threshold while the central redshift must satisfy $z_{0}<z_{0}^{\rm c}$ because configurations with $z_{0}>z_{0}^{\rm c}$ exist mathematically but are unstable.