2 Antonov instability for a relativistic gas sphere

   
2.1 The relativistic gas

We consider a system of N particles, each of mass m, in gravitational interaction. We allow the speed of the particles to be close to the velocity of light so that special relativity must be taken into account. However, in this first approach, we shall treat the gravitational field within the framework of Newtonian mechanics. This procedure is permissible if the typical size of the system is much larger than the Schwarzschild radius (see, e.g., Chandrasekhar 1942). Let $f(\vec{r},\vec{p},t)$ denote the distribution function of the system, i.e. $f(\vec{r},\vec{p},t){\rm d}^{3}\vec{r}{\rm d}^{3}\vec{p}$ gives the average mass of particles whose positions and momenta are in the cell $(\vec{r},\vec{p}; \vec{r}+{\rm d}^{3}\vec{r},\vec{p}+{\rm d}^{3}\vec{p})$ at time t. The integral of f over the momenta determines the spatial density

$$\rho=\int f {\rm d}^{3}\vec{p},$$ (1)

and the total mass is expressed as

$$M=\int \rho {\rm d}^{3}\vec{r}.$$ (2)

In a meanfield approximation, the total energy is given by

$$E=\int {f\over m} \epsilon \ {\rm d}^{3}\vec{r}{\rm d}^{3}\vec{p}+{1\over 2}\int\rho\Phi {\rm d}^{3}\vec{r}=K+W,$$ (3)

where K and W are the kinetic and potential energy respectively. According to the theory of special relativity, the energy of a particule reads

$$\epsilon=m c^{2}\biggl\lbrace \biggr (1+{p^{2}\over m^{2}c^{2}}\biggr )^{1/2}-1\biggr \rbrace\cdot$$ (4)

This expression does not include the rest mass, so it reduces to the usual kinetic energy ${p^{2}\over 2m}$ in the Newtonian limit. All the effects of gravity are incorporated in the potential energy which contains the gravitational potential $\Phi$ related to the density by the Newton-Poisson equation

$$\Delta\Phi=4\pi G\rho.$$ (5)

We now ask which configuration maximizes the Boltzmann entropy

$$S=-k\int {f\over m}\ln {f\over m} {\rm d}^{3}\vec{r}{\rm d}^{3}\vec{p},$$ (6)

subject to the conservation of mass M and energy E. To that purpose, we proceed in two steps. In this section, we maximize S[f] at fixed E, M and $\rho(\vec{r})$. This provides an optimal distribution function expressed in terms of $\rho$. Then, in Sect. 2.2, we maximize $S_{\max}[\rho]$ at fixed M and E. This optimization problem is not trivial and will be further discussed in Sects. 2.4-2.6.

Since the gravitational potential can be deduced from the density by solving the Poisson Eq. (5), maximizing S at fixed E, M and $\rho(\vec{r})$ is equivalent to maximizing S at fixed $\rho(\vec{r})$ and K. Writing the variational principle in the form

$$\delta S-k\beta\delta K-\int\lambda(\vec{r})\delta\rho {\rm d}^{3}\vec{r}=0,$$ (7)

we obtain the optimal distribution function

$$f(\vec{r},\vec{p})=A(\vec{r}){\rm e}^{-\beta \epsilon},$$ (8)

which is a global entropy maximum with the previous constraints. Equation (8) is the relativistic Maxwell-Boltzmann distribution with an inverse temperature

$$\beta={1\over kT}\cdot$$ (9)

The Lagrange multipliers $A(\vec{r})$ and $\beta$ must be related to the constraints $\rho(\vec{r})$ and K. As discussed in detail by Chandrasekhar (1942), these relations can be expressed in terms of the modified Bessel functions

$$K_{n}(z)=\int_{0}^{+\infty}{\rm e}^{-z\cosh\theta}\cosh(n\theta){\rm d}\theta.$$ (10)

Using Eqs. (1) (8) (4) and introducing the Juttner transformation ${p\over mc}=\sinh\theta$ in the integral, the density can be written

$$\rho(\vec{r})={4\pi m^{3} c^{3}\over x}A(\vec{r})K_{2}(x){\rm e}^{x},$$ (11)

where we have introduced the dimensionless parameter

$$x=\beta mc^{2}={mc^{2}\over kT}\cdot$$ (12)

This parameter quantifies the importance of relativistic effects. The classical limit corresponds to $x\rightarrow +\infty$ ( $kT\ll mc^{2}$) and the ultra-relativistic limit to $x\rightarrow 0$ ( $kT\gg mc^{2}$). The distribution function (8) can now be expressed in terms of the density as

$$f(\vec{r},\vec{p})={x \over 4\pi m^{3}c^{3}} {{\rm e}^{-x}\over K_{2}(x)}\rho(\vec{r}) {\rm e}^{-\beta\epsilon}.$$ (13)

In the classical limit $(x\rightarrow +\infty)$ we recover the standard formula

$$f(\vec{r},\vec{p})={1\over (2\pi mkT)^{3/2}}\rho(\vec{r}){\rm e}^{-{p^{2}\over 2mkT}},$$ (14)

and in the ultra-relativistic limit ( $x\rightarrow 0$), we get

$$f(\vec{r},\vec{p})={c^{3}\over 8\pi k^{3}T^{3}}\rho(\vec{r}){\rm e}^{-{pc\over kT}}.$$ (15)

Similarly, after some elementary transformations, the kinetic energy can be expressed in terms of the normalized inverse temperature x by

$$K={\cal F}(x)Mc^{2},\qquad {\cal F}(x)={3K_{3}(x)+K_{1}(x)\over 4K_{2}(x)}-1.$$ (16)

Using the recursion formula

$$K_{n-1}(x)-K_{n+1}(x)=-{2n\over x}K_{n}(x),$$ (17)

the function ${\cal F}(x)$ can be written in the equivalent form

$${\cal F}(x)={K_{1}(x)\over K_{2}(x)}+{3\over x}-1.$$ (18)

It has the asymptotic behaviors

$${\cal F}(x)\sim {3\over 2x}\qquad (x\rightarrow +\infty),$$ (19)

$${\cal F}(x)\sim {3\over x}\qquad (x\rightarrow 0).$$ (20)

When substituted in Eq. (16), we recover the usual expressions of the kinetic energy $K={3\over 2}NkT$ in the classical limit and K=3NkT in the ultra-relativistic limit.

We can now express the entropy in terms of the density $\rho$ and the inverse temperature x. Substituting the optimal distribution function (13) in Eq. (6), we get, up to an additional constant

$$S=kN{\cal G}(x)-k\int{\rho\over m}\ln {\rho\over m} {\rm d}^{3}\vec{r},$$ (21)

where

$${\cal G}(x)=x\lbrack {\cal F}(x)+1\rbrack +\ln K_{2}(x) -\ln x.$$ (22)

The function ${\cal G}(x)$ has the asymptotic behaviors

$${\cal G}(x)\sim -{3\over 2}\ln x \qquad (x\rightarrow +\infty),$$ (23)

$${\cal G}(x)\sim -3\ln x \qquad (x\rightarrow 0).$$ (24)

The thermal contribution to the entropy in the classical limit is $S_{\rm th}={3\over 2}kN\ln T$ and in the ultra-relativistic limit $S_{\rm th}={3}kN\ln T$.

There exists a general relation between the derivatives of ${\cal F}$ and ${\cal G}$ that we shall need in the following. Differentiating Eq. (22) with respect to x and using the identity

$$K'_{n}(x)=-K_{n-1}(x)-{n\over x}K_{n}(x),$$ (25)

for n=2, we find that

$${\cal G}'(x)=x{\cal F}'(x).$$ (26)

   
2.2 First and second order variations of entropy

In the preceding section, we have expressed the entropy and the kinetic energy in terms of the density $\rho(\vec{r})$ and the temperature T (through the variable x). We now wish to maximize the entropy $S[\rho]$ at fixed E and M. For convenience, we shall introduce a new variable $y={\cal F}(x)$. In terms of this variable, the total energy and the entropy can be written

$$E=yMc^{2}+{1\over 2}\int\rho\Phi {\rm d}^{3}\vec{r},$$ (27)

$$S=kN{\cal G}(x(y))-k\int{\rho\over m}\ln {\rho\over m} {\rm d}^{3}\vec{r}.$$ (28)

We can now determine the variations of S around a given density profile $\rho(\vec{r})$. To second order in the expansion, we get

 

$$\delta S=kN{{\rm d}{\cal G}\over {\rm d}y}\delta y+kN {{\rm d}^{2... ... d}^{3}\vec{r} -{k\over m}\int {(\delta\rho)^{2}\over 2\rho}{\rm d}^{3}\vec{r}.$$ (29)

Using the identity (26), we find that

$${{\rm d}{\cal G}\over {\rm d}y}={{\rm d}{\cal G}\over {\rm d}x}{{\rm d}x\over {\rm d}y}={{\cal G}'(x)\over {\cal F}'(x)}=x.$$ (30)

Differentiating one more time with respect to y, we obtain

$${{\rm d}^{2}{\cal G}\over {\rm d}y^{2}}={{\rm d}x\over {\rm d}y}={1\over {\cal F}'(x)}\cdot$$ (31)

Substituting the above results in Eq. (29), we get

 

$$\delta S=kN x\delta y+kN{1\over {\cal F}'(x)} {(\delta y)^{2}\ove... ...m d}^{3}\vec{r}-{k\over m}\int {(\delta\rho)^{2}\over 2\rho}{\rm d}^{3}\vec{r}.$$ (32)

We now need to express the variation $\delta y$ in terms of $\delta \rho$. From the conservation of energy, we have the exact identity

$$0=\delta E=Mc^{2}\delta y+\int \Phi\delta \rho {\rm d}^{3}\vec{r}+{1\over 2}\int\delta\rho\delta \Phi {\rm d}^{3}\vec{r}.$$ (33)

Substituting the foregoing expression for $\delta y$ from Eq. (33) in Eq. (32), we obtain

 

$$\delta S=-{1\over T}\int \Phi\delta\rho {\rm d}^{3}\vec{r}-{1\ove... ... d}^{3}\vec{r} -{k\over m}\int {(\delta\rho)^{2}\over 2\rho}{\rm d}^{3}\vec{r}.$$ (34)

Introducing a Lagrange multiplier $\alpha$ to satisfy the conservation of mass, the condition that S is an extremum is written (to first order)

 

$$0=\delta S-\alpha\delta M=-\int\biggl\lbrack {\Phi\over T} +{k\ov... ...1+\ln {\rho\over m}\biggr )+\alpha \biggr\rbrack \delta\rho {\rm d}^{3}\vec{r}.$$ (35)

This condition must be satisfied for any variations $\delta \rho$. This yields the Boltzmann distribution

 

$$\rho=A{\rm e}^{-{m\Phi\over kT}},$$ (36)

like for a classical gas. The condition that the critical point (36) is an entropy maximum requires that

 

$$\delta^{2} S=-{1\over 2T}\int\delta\rho\delta\Phi {\rm d}^{3}\vec... ...\biggr )^{2} -{k\over m}\int {(\delta\rho)^{2}\over 2\rho}{\rm d}^{3}\vec{r}<0,$$ (37)

for any variation $\delta \rho$ that conserves mass to first order.

   
2.3 The Virial theorem for a relativistic gas

In this section, we derive the form of the Virial theorem appropriate to an isothermal gas described in the framework of special relativity. Quite generally, the potential energy of a self-gravitating system can be expressed in the form (Binney & Tremaine 1987)

 

$$W=-\int\rho \vec{r}\nabla\Phi {\rm d}^{3}\vec{r}.$$ (38)

If the system is in hydrostatic equilibrium, then

 

$$\nabla p=-\rho\nabla\Phi.$$ (39)

Substituting this identity in Eq. (38) and integrating by parts, we get

 

$$W=\int \vec{r}\nabla p {\rm d}^{3}\vec{r}=\oint p\vec{r}{\rm d}\vec{S}-3\int p {\rm d}^{3}\vec{r}.$$ (40)

If the pressure $p_{\rm b}$ on the boundary of the system is uniform, we can write

 

$$\oint p\vec{r}{\rm d}\vec{S}=p_{\rm b}\oint \vec{r}{\rm d}\vec{S}=p_{\rm b}\int \nabla .\vec{r} {\rm d}^{3} \vec{r}=3p_{\rm b}V,$$ (41)

where V is the total volume of the system. Therefore, for any system in hydrostatic equilibrium, one has

 

$$W=3p_{\rm b}V-3\int p{\rm d}^{3}\vec{r},$$ (42)

which can be considered as the general form of the Virial theorem for self-gravitating systems.

Now, the pressure of an ideal gas can be expressed as (Chandrasekhar 1942)

 

$$p={1\over 3}\int {f\over m}p{\partial\epsilon\over\partial p}{\rm d}^{3}\vec{p}.$$ (43)

For a relativistic gas described by the distribution function (8), we get

 

$$p=-{4\pi\over 3 m\beta} A(\vec{r})\int_{0}^{+\infty}{\partial\over\partial p}\left({\rm e}^{-\beta\epsilon}\right)p^{3}{\rm d}p.$$ (44)

Integrating by parts, we obtain

 

$$p={\rho\over m\beta}={\rho\over m}{k}T.$$ (45)

Therefore, the equation of state for a (non quantum) relativistic gas is the same as for a classical gas. Written in the form

 

$$W+3NkT=3p_{\rm b}V,$$ (46)

the Virial theorem (42) also has the same form as for an isothermal classical gas. However, the Virial theorem is usually expressed in terms of the kinetic energy instead of the temperature. Therefore, the appropriate form of the relativistic Virial theorem reads

 

$$W+{3Mc^{2}\over {\cal F}^{-1}({K\over Mc^{2}})}=3p_{\rm b}V.$$ (47)

In the classical limit ( $x\rightarrow +\infty$), it reduces to the well-known formula

 

$$W+2K=3p_{\rm b}V,$$ (48)

and in the ultra-relativistic limit ( $x\rightarrow 0$), we get

 

$$W+K=3p_{\rm b}V.$$ (49)

It should be emphasized that Eq. (47) is only valid for a relativistic gas in thermal equilibrium.

   
2.4 The equilibrium phase diagram

Since the equation of state for a relativistic gas is the same as for a classical gas, the equilibrium configurations of such systems correspond to the isothermal gas spheres extensively described in the monograph of Chandrasekhar (1942). Only the onset of the gravitational instability will be modified by relativistic effects.

For non-rotating systems, the equilibrium states are expected to be spherically symmetric. In that case, the Poisson Eq. (5) together with the Boltzmann distribution (36) yield the second order differential equation

 

$${1\over r^{2}}{{\rm d}\over {\rm d}r}\biggl (r^{2}{{\rm d}\Phi\over{\rm d}r}\biggr )=4\pi GA{\rm e}^{-\beta m\Phi}.$$ (50)

This equation can also be deduced from the condition of hydrostatic equilibrium (39) when the pressure is related to the density according to the equation of state (45). It is well-known that the density profile of such isothermal configuations behaves like $\rho\sim r^{-2}$ at large distances so that their total mass is infinite. Following Antonov (1962), we shall avoid this infinite mass problem by confining artificially the system within a spherical box of radius R. It is only under this simplifying assumption that a rigorous thermodynamics of self-gravitating systems can be carried out (see, e.g., Padmanabhan 1990; Chavanis et al. 2001). This procedure is justified physically by the realization that a distribution of matter never extends to infinity so R represents an upper cut-off at which other processes intervene to limitate the spatial extent of the system. Of course, different cut-offs are possible but fixing the volume is consistent with the traditional viewpoint of statistical mechanics and it is sufficient to capture the essential physics of the problem (see the different comparisons of truncated models performed by Katz 1980).

We now wish to determine the equilibrium phase diagram of a relativistic isothermal gas. To that purpose, we introduce the function $\psi=\beta m (\Phi-\Phi_{0})$ where $\Phi_{0}$ is the gravitational potential at r=0. Then, the density field (36) can be written

$$\rho=\rho_{0}{\rm e}^{-\psi},$$ (51)

where $\rho_{0}$ is the central density. Introducing the notation $\xi=(4\pi G\beta m \rho_{0})^{1/2}r$, the Boltzmann-Poisson Eq. (50) reduces to the standard Emden form

$${1\over\xi^{2}}{{\rm d}\over {\rm d}\xi}\biggl (\xi^{2}{{\rm d}\psi\over {\rm d}\xi}\biggr )={\rm e}^{-\psi}.$$ (52)

Equation (52) has a simple analytic solution, the singular sphere

$${\rm e}^{-\psi_{\rm s}}={2\over\xi^{2}},$$ (53)

whose central density is infinite. The regular solution of Eq. (52) satisfying the boundary conditions

$$\psi(0)=\psi'(0)=0,$$ (54)

at the center of the sphere must be computed numerically. In the case of bounded isothermal systems, we must stop the integration at the normalized box radius

$$\alpha=(4\pi G\beta m\rho_{0})^{1/2}R.$$ (55)

We shall now relate the parameter $\alpha$ to the temperature and the energy. According to the Poisson Eq. (5), we have

 

$$GM=\int_{0}^{R} 4\pi G\rho r^{2}{\rm d}r= \int_{0}^{R}{{\rm d}\ov... ...rm d}r}\biggr ){\rm d}r=\biggl (r^{2}{{\rm d}\Phi\over {\rm d}r}\biggr )_{r=R},$$ (56)

which is just a particular case of the Gauss theorem. Introducing the dimensionless variables defined previously, we get

$$\eta\equiv {\beta GMm\over R}=\alpha\psi'(\alpha).$$ (57)

The relation between $\alpha$ and the normalized temperature $\eta$ is not affected by special relativity.

For the total energy, using the Virial theorem (46) and the expression (16) for the kinetic term, we have

$$E=K+W={\cal F}(x)Mc^{2}-{3N\over\beta}+3p(R)V.$$ (58)

Now, the pressure at the boundary of the domain can be written

$$p(R)={\rho(R)\over m\beta}={\rho_{0}{\rm e}^{-\psi(\alpha)}\over m\beta}\cdot$$ (59)

Expressing the central density in terms of $\alpha$, using Eq. (55), we have equivalently

$$p(R)={\alpha^{2}\over 4\pi G\beta^{2}m^{2}R^{2}}{\rm e}^{-\psi(\alpha)}.$$ (60)

The total energy therefore reads

$$\Lambda\equiv -{ER\over GM^{2}}=-{Rc^{2}\over GM}{\cal F}(x)+... ...psi'(\alpha)}-{{\rm e}^{-\psi(\alpha)}\over\psi'(\alpha)^{2}},$$ (61)

where we have used Eq. (57) to eliminate the temperature in the last two terms.

It will be convenient in the following to introduce the parameter

$$\mu={Rc^{2}\over GM}\equiv {2R\over R_{\rm S}^{\rm class}},$$ (62)

which is twice the ratio between the system radius R and the "classical'' Schwarzschild radius

$$R_{\rm S}^{\rm class}={2GM\over c^{2}},$$ (63)

constructed with the total mass M of the system. Clearly, our semi-relativistic treatment, which ignores general relativity, is only valid for $\mu\gg 1$. However, in our rather formal analysis, we shall treat $\mu$ as a free parameter varying in the range $0\le\mu<+\infty$. The relativistic parameters $\mu$ and x are related to each other by

$$x=\mu\eta.$$ (64)

Therefore, in accordance with Eqs. (61) and (64), the relation between the parameter $\alpha$ and the normalized energy $\Lambda$ takes the form

$$\Lambda\equiv -{ER\over GM^{2}}=-\mu {\cal F}\lbrack \mu\alph... ...(\alpha)}-{{\rm e}^{-\psi(\alpha)}\over\psi'(\alpha)^{2}}\cdot$$ (65)

In the classical limit ( $\mu\rightarrow +\infty$), we recover the result of Lynden-Bell & Wood (1968)

$$\Lambda={3\over 2\alpha\psi'(\alpha)}-{{\rm e}^{-\psi(\alpha)}\over\psi'(\alpha)^{2}},$$ (66)

and in the formal limit $\mu\rightarrow 0$, we get

$$\Lambda=-{{\rm e}^{-\psi(\alpha)}\over\psi'(\alpha)^{2}}\cdot$$ (67)

Figure 1: Equilibrium phase diagram for isothermal gas spheres described in the framework of special relativity. Relativistic effects shift the onset of instability to larger energies.

In Fig. 1, we have represented the equilibrium phase diagram $\beta-E$ for different values of the relativistic parameter $\mu$. We see that the effect of relativity is to shift the spiral to the left. Therefore, the "gravothermal catastrophe'' (corresponding to the absence of equilibrium below a critical energy) occurs sooner than in the Newtonian case. The critical energy is ploted as a function of the relativistic parameter $\mu$ in Fig. 2. For $\mu\rightarrow +\infty$, we recover the classical result of Antonov $\Lambda_{\rm c}(+\infty)=0.335$ and, in the formal limit $\mu\rightarrow 0$, we get $\Lambda_{\rm c}(0)=-0.345$. Clearly, the spiral is not destroyed by relativistic effects. Only is its shape slightly modified: the relativistic spirals are more "stretched'' than the classical one. Note that the maximum value of $\eta$ is independant on $\mu$ and is equal to its classical value $\eta_{\rm c}=2.52$. On the other hand, substituting the expression (53) for the singular sphere in Eqs. (57) and (65), we find that the center of the spiral is determined by the equations

$$\eta_{\rm s}=2,$$ (68)

$$\Lambda_{\rm s}(\mu)=1-\mu{\cal F}(2\mu).$$ (69)

In the limit $\mu\rightarrow +\infty$, $\Lambda_{\rm s}(+\infty)={1\over 4}$ and in the limit $\mu\rightarrow 0$, $\Lambda_{\rm s}(0)=-{1\over 2}$.

The spiral is parametrized by the normalized box radius $\alpha$that goes from 0 (ordinary gas) to $+\infty$ (singular sphere) when we spiral inwards. If one prefers, we can use a parametrization in terms of the density contrast

$${\cal R}={\rho_{0}\over \rho(R)}={\rm e}^{\psi(\alpha)},$$ (70)

that goes from 1 to $+\infty$. In Fig. 3, we plot the critical density contrast (corresponding to $\Lambda _{\rm c}$) as a function of the relativistic parameter $\mu$. For $\mu\rightarrow +\infty$, we recover the classical value ${\cal R}_{\rm c}(+\infty)=709$(and $\alpha_{\rm c}(+\infty)=34.4$). For $\mu =0$, we get ${\cal R}_{\rm c}(0)=132$ (and $\alpha_{\rm c}(0)=16.0$). It is found that instability occurs for smaller density contrasts when relativity is accounted for.

Figure 2: Critical energy $\Lambda _{\rm c}$ as a function of the relativistic parameter $\mu$.

Figure 3: Critical density contrast ${\cal R}_{\rm c}$ as a function of the relativistic parameter $\mu$.

   
2.5 The Milne variables

It will be convenient in the following to introduce the Milne variables (u,v) defined by (Chandrasekhar 1942)

$$u={\xi {\rm e}^{-\psi}\over\psi'},\qquad {\rm and}\qquad v=\xi\psi'.$$ (71)

Taking the logarithmic derivative of u and v with respect to $\xi$ and using Eq. (52), we get

$${1\over u}{{\rm d}u\over {\rm d}\xi}={1\over\xi}(3-v-u),$$ (72)

$${1\over v}{{\rm d}v\over {\rm d}\xi}={1\over\xi}(u-1).$$ (73)

Taking the ratio of these equations, we find that the variables u and v are related to each other by a first order differential equation

$${u\over v}{{\rm d}v\over {\rm d}u}=-{u-1\over u+v-3}\cdot$$ (74)

The solution curve in the (u,v) plane is well-known and is represented in Fig. 4. Its striking oscillating behavior has been described by a number of authors (see in particular Chandrasekhar 1942). We refer to Padmanabhan (1989) and Chavanis (2001) for the description of its main characteristics in connexion with the present work.

It turns out that the normalized temperature and the normalized energy can be expressed very simply in terms of the values of u and v at the normalized box radius $\alpha$. Indeed, writing $u_{0}=u(\alpha)$ and $v_{0}=v(\alpha)$ and using Eqs. (57), (65), we get

$$\eta=v_{0},$$ (75)

$$\Lambda=-\mu{\cal F}(\mu v_{0})+{3\over v_{0}}-{u_{0}\over v_{0}}\cdot$$ (76)

The intersections between the curves defined by Eqs. (75), (76) and the spiral in the (u,v) plane determine the values of $\alpha$ corresponding to a given temperature or energy. Considering Eq. (75), we find that there is no intersection for $\eta={\beta GM\over R}>v_{\max}=2.52$. Therefore, below a critical temperature $kT_{\rm c}={GmM\over 2.52 R}$, an isothermal sphere is expected to collapse. This classical result is not altered by special relativity. Considering now the microcanonical ensemble, we first note that Eq. (76) can be rewritten

$$u_{0}=3-\mu v_{0}{\cal F}(\mu v_{0})-\Lambda v_{0}.$$ (77)

In the classical limit ( $\mu\rightarrow +\infty$) it reduces to the straight line found by Padmanabhan (1989)

$$u_{0}={3\over 2}-\Lambda v_{0},$$ (78)

and in the limit $\mu\rightarrow 0$, we find another straight line

$$u_{0}=-\Lambda v_{0}.$$ (79)

The curve (77) is ploted in Fig. 4 for a fixed value of $\mu$ and for different values of $\Lambda$. For $\Lambda>\Lambda_{\rm c}(\mu)$ there is no intersection, for $\Lambda=\Lambda_{\rm c}(\mu)$ the curve (77) is tangent to the spiral and for $\Lambda<\Lambda_{\rm c}(\mu)$ there are one or several intersections. We recover therefore by this graphical construction the existence of a critical energy below which no hydrostatic equilibrium can exist for isothermal spheres. In Fig. 5, we plot the same diagram as Fig. 4 but for different values of $\mu$and, in each case, for the critical energy $\Lambda _{\rm c}(\mu )$. The intersection with the spiral determines the value of $\alpha$ at the critical point. This figure confirms that instability occurs sooner (i.e., for smaller values of $\alpha$ or smaller density contrasts) when relativistic effects are taken into account.

Figure 4: The (u,v) plane. All isothermal spheres must necessarily lie on the spiral. There exists solutions in the canonical ensemble only for $\eta <2.52$. In the microcanonical ensemble, the critical energy depends on the relativistic parameter $\mu$. For the value $\mu =1$ adopted in the figure, there exists solutions only for $\Lambda <\Lambda _{\rm c}=0.06$.

Figure 5: Same as Fig. 4 but for different values of $\mu$ and, in each case, for the critical parameter $\Lambda _{\rm c}(\mu )$ above which there is no equilibrium solution.

   
2.6 The condition of thermodynamical stability

We now address the question of thermodynamical stability. Let us recall that the isothermal spheres lying on the spirals of Fig. 1 are critical points of entropy but they are not necessarily entropy maxima. To determine whether they are local entropy maxima or saddle points, we must examine the sign of the second variations of entropy given by Eq. (37). In fact, this condition of stability has the same form as in the classical case provided that we make the substitution

$$-{1\over 3NkT^{2}}\ \rightarrow \ {k\over 2Mmc^{4}}{1\over {\cal F}'(x)}\cdot$$ (80)

Therefore, the analysis of Padmanabhan (1989) for the classical Antonov instability can be extended straightforwardly. Introducing the mass perturbation $q(r)\equiv \delta M(r)=\int_{0}^{r}4\pi r^{'2}\delta\rho(r'){\rm d}r'$ within the sphere of radius r such that

$$\delta\rho={1\over 4\pi r^{2}}{{\rm d}q\over {\rm d}r},$$ (81)

the second variations of entropy can be put in a quadratic form

$$\delta^{2}S=\int_{0}^{R}\int_{0}^{R}{\rm d}r{\rm d}r' q(r)K(r,r')q(r'),$$ (82)

with

 

$$K(r,r')={k\over 2Mmc^{4}}{1\over {\cal F}'(x)}{{\rm d}\Phi\over {... ...biggl ({1\over 4\pi\rho r^{2}}{{\rm d}\over {\rm d}r}\biggr )\biggr\rbrack\cdot$$ (83)

Clearly, the conservation of mass imposes the boundary conditions q(0)=q(R)=0. The problem of stability can therefore be reduced to the study of the eigenvalue equation

$$\int_{0}^{R}K(r,r')F_{\lambda}(r'){\rm d}r'=\lambda F_{\lambda}(r),$$ (84)

with the boundary conditions $F_{\lambda}(0)=F_{\lambda}(R)=0$. If all the eigenvalues are negative, then $\delta^{2}S<0$ and the critical point is a local entropy maximum. If one eigenvalue is positive, the critical point is an unstable saddle point. The point of marginal stability is determined by the condition that the largest eigenvalue is equal to zero ($\lambda=0$). We thus have to solve the differential equation

 

$$\biggl \lbrack {k\over m}{{\rm d}\over {\rm d}r}\biggl ({1\over 4... ...ck F(r) =-{k\over Mmc^{4}}{1\over {\cal F}'(x)}V{{\rm d}\Phi\over {\rm d}r}(r),$$ (85)

with

$$V=\int_{0}^{R}{{\rm d}\Phi\over {\rm d}r}(r')F(r'){\rm d}r',$$ (86)

and F(0)=F(R)=0. Introducing the dimensionless variables defined in Sect. 2.4, it can be rewritten

$$\biggl\lbrack {{\rm d}\over {\rm d}\xi}\biggl ({{\rm e}^{\psi... ...i^{2}}\biggr\rbrack F(\xi)=\chi {{\rm d}\psi\over {\rm d}\xi},$$ (87)

with

$$\chi=-{1\over x^{2}{\cal F}'(x)} {1\over \alpha^{2}\psi'(\alp... ...{\alpha}{{\rm d}\psi\over {\rm d}\xi}(\xi')F(\xi'){\rm d}\xi',$$ (88)

and $F(0)=F(\alpha)=0$. As shown by Padmanabhan (1989), the solutions of the differential equation (87) can be expressed in terms of the solutions of the Emden Eq. (52) as

$$F(\xi)=\chi\biggl\lbrack {1\over 1-u_{0}}(\xi^{3}{\rm e}^{-\psi}-\xi^{2}\psi')+\xi^{2}\psi'\biggr\rbrack.$$ (89)

We can check that this function satisfies the boundary conditions $F(0)=F(\alpha)=0$. The point of marginal stability is obtained by substituting the solution (89) in Eq. (88). The integrations can be carried out and the solutions expressed in terms of the Milne variables u₀ and v₀ (see Padmanabhan 1989 for more details). In our semi-relativistic treatment, we obtain

$$2u_{0}^{2}+u_{0}v_{0}-7u_{0}+3-\mu^{2}v_{0}^{2}{\cal F}'(\mu v_{0})(u_{0}-1)=0,$$ (90)

where we have used Eqs. (64) and (75) to simplify the last term.

In the classical limit $\mu\rightarrow +\infty$ we recover the result of Padmanabhan (1989)

 

4u₀²+2u₀v₀-11u₀+3=0, (91)

and in the formal limit $\mu\rightarrow 0$, we find

 

2u₀+v₀-4=0. (92)

The intersections between the curve (90) and the spiral in the (u,v) plane determine the values of $\alpha$ for which a new mode of stability is lost (i.e., a new eigenvalue $\lambda$ becomes positive). Since the curve (90) passes through the singular sphere $(u_{\rm s},v_{\rm s})=(1,2)$ at the center of the spiral, there is an infinity of intersections. The first intersection (for which $\alpha$ is minimum) corresponds to the point of marginal stability denoted $\alpha _{1}$. We can show that the points determined by Eq. (90) are precisely those for which $\Lambda$ is extremum in agreement with the turning point analysis of Katz (1978). Indeed, differentiating the expression (76) for $\Lambda$ with respect to $\alpha$, we get

$${{\rm d}\Lambda\over {\rm d}\alpha}=-\mu^{2}{\cal F}'(\mu v_{... ...}+{u_{0}\over v_{0}^{2}}{{\rm d}v_{0}\over {\rm d}\alpha}\cdot$$ (93)

Using Eqs. (72), (73), we obtain

$${{\rm d}\Lambda\over {\rm d}\alpha}={1\over \alpha v_{0}}\big... ...-7u_0+3-\mu^{2}v_0^{2}{\cal F}'(\mu v_0)(u_0 -1)\biggr \rbrack$$ (94)

and we check that the condition ${{\rm d}\Lambda\over {\rm d}\alpha}=0$ is equivalent to Eq. (90).

Figure 6: Graphical construction to determine the nodes of the perturbation profile $\delta \rho$ at the point of marginal stability. The construction is done explicitly for $\mu =0$ (formally), for which $\alpha _{1}=16.0$ and u0=0.817. There is only one zero satisfying $\xi _{1}<\alpha _{1}$ so that the perturbation profile does not present a "core-halo'' structure. This property is maintained until $\mu >1.61$.

Figure 7: Density perturbation profile at the point of marginal stability as a function of the relativistic parameter.

It is also easy to determine the form of the perturbation that triggers the instability at the critical point $\Lambda _{\rm c}(\mu )$. According to Eq. (81), the eigenfunction associated with the eigenvalue $\lambda=0$ can be written

$${\delta\rho\over\rho_{0}}={1\over 4\pi\xi^{2}}{{\rm d}F\over {\rm d}\xi},$$ (95)

where $F(\xi)$ is given by Eq. (89). Simplifying the derivative with the aid of Eq. (52), we can express the perturbation profile in terms of the Milne variables (71) as

$${\delta\rho\over\rho}={\chi\over 4\pi}{1\over 1-u_{0}}(3-v-u_{0}).$$ (96)

The qualitative behavior of the perturbation profile can be studied without numerical integration by a graphical construction (see Padmanabhan 1989). The density perturbation $\delta \rho$ becomes zero at the point(s) $\xi_{i}$ such that $u_{0}=3-v(\xi_{i})$. In Fig. 6, we first draw the line u=3-v. This line passes through the singular sphere $(u_{\rm s}=1$, $v_{\rm s}=2)$ and also cuts the spiral at the points of vertical tangent (see Eq. (72)). In particular, the first intersection corresponds to $\alpha_{*}=22.5$ and (u_*,v_*)=(0.793,2.208). Then, we draw the line $u=u_{0}=u(\alpha_{1})$. The intersection between these two lines determines $v(\xi_{i})$. The intersection between $v=v(\xi_{i})$ and the spiral determines the zeros of $\delta \rho$. For $\alpha_{1}>\alpha_{*}$, there are two intersections satisfying $\xi_{i}<\alpha_{1}$ so that the perturbation profile presents a "core-halo'' structure. This is the case in particular in the classical limit $\mu\rightarrow +\infty$ for which $\alpha_{1}=34.4$ (see Padmanabhan 1989). By contrast, for $\alpha_{1}<\alpha_{*}$, there is only one intersection satisfying $\xi _{1}<\alpha _{1}$ so that the perturbation profile does not present a "core-halo'' structure. This is the case in particular in the (formal) limit $\mu =0$ for which $\alpha _{1}=16.0$. The "core-halo'' structure disappears for $\alpha_{1}=\alpha_{*}=22.5$ corresponding to a relativistic parameter $\mu_{*}=1.61$. Since our study is valid for $\mu\gg 1$, we deduce that the density perturbation profile always presents a "core-halo'' structure in the cases of physical interest. However, relativistic effects have the tendency to reduce the extent of the halo (see Fig. 7).

The previous results are valid in the microcanonical ensemble in which the energy is fixed. In the canonical ensemble, we must consider maxima of the free energy $J=S-\beta E$ at fixed temperature. In that case, the condition of stability is given by Eq. (85) with V=0. Since the relativistic function ${\cal F}(x)$ does not appear anymore in the equations, we conclude that special relativity does not change the classical results in the canonical ensemble. In particular, the perturbation profile does not present a "core-halo'' structure at the critical point $\eta_{\rm c}$ (Chavanis 2001).