text:
# Thermodynamics and Self-Gravitating Systems

## 1 Introduction

A strident aspect that puts self-gravitating systems in a class of their own is that they cannot be small sub-systems of a large ensemble except in very imaginary situations. A small sub-system of a large self-gravitating system has its energy almost entirely determined by the large system. Moreover, the sum of energies of all the small sub-systems is not equal to the total energy of the system in contrast with a basic tenet of classical thermodynamics. This has far reaching consequences. Much of the rich harvest of classical and quantum thermodynamics as exposed in Landau and Lifshitz’s Statistical Physics  is useless in gravitational thermodynamics. For instance stable and unstable isolated self-gravitating systems may have heat capacities of both signs. This is also the case for a system in a heat bath. Thus $`C_V>0`$ is no more a signature of stable thermal states like in classical thermodynamic ensembles. In this respect gravitational thermodynamics appears slightly exotic though not to astronomers as pointed out by Lynden-Bell  in his review on progress in understanding of gravitational thermodynamics. It is also what makes thermodynamics of self-gravitating systems interesting.

The paper is based on the Gibbs density of states function in the meanfield approximation. Astronomers and astrophysicists are more familiar with Boltzmann entropy. Gibbs definition is more general in that it provides Boltzmann entropy as a first approximation as well as stability conditions and fluctuations. It has a wider range of applications and extends to meanfield theory. The Gibbs approach is in our view the least frustrating way to understand statistical thermodynamics and provides the best logically connected approach to the subject.

Section 1 begins with a review of statistical equilibrium theory in the simple case of $`N`$ identical point masses in an external field first and gravitationally interacting particles in a finite volume next. The case $`N=2`$ is interesting and has been studied in detail by Padmanabhan . Here we consider $`N>>2`$ and situations in which particle pair formation has a negligible role on relevant time scales. In this case meanfield theory is a very good approximation. A steepest descent calculation gives equilibrium configurations of stable as well as metastable states. Metastable states are particularly interesting in gravitational thermodynamics222Most globular clusters  and rich clusters of galaxies  have their cores in a quasi thermalized metastable state.. The method provides also stability conditions.

In section 2 we deal with thermodynamic functions and stability conditions in general. The latter are first examined with the help of Gibbs’ density of states function but are subsequently translated into conditions involving only thermodynamic functions. The theory applies to isolated non-extensive systems, the type of system considered in section 1, but is of much broader applicability. It predicts in particular that under very general assumptions stable thermal equilibrium configurations of isolated systems near instability have always negative heat capacities which turn to be positive when systems become unstable! In classical thermodynamics a small thermally stable subsystem of an ensemble has always a positive heat capacity. These results imply that stable isolated systems near instability are necessarily unstable if put in a heat bath. This well known “nonequivalence of ensembles” has no analogue in classical thermodynamics.

Section 3 develops the theory of fluctuations. We deal with fluctuations of temperature as a generic example. Fluctuations of other thermodynamic functions can be dealt with in a similar way. The theory has also broad applicability; again results are expressed in terms of thermodynamic functions and hold only near instability. For nonequivalent ensembles the behavior of equilibrium near instability departs considerably from classical systems and the effect of fluctuations is consequently quite different.

Section 4 gives two examples of application with some detail: isothermal spheres in various ensembles and liquid ellipsoids. The purpose here is to illustrate the power and limitations of thermodynamics in self-gravitating systems. The last sub-section reviews briefly various other applications in astronomy where thermodynamics has been used or might have been with considerable benefit.

Each section has its own summary, conclusion and/or comments.

## 2 Statistical thermodynamics

It may interest the reader to see the explicit connection between the entropies of Gibbs and Boltzmann by starting with the simpler case of a system in an external field first. Here some of the mathematical techniques appear in a more transparent form. We go then over to self-gravitating systems in the meanfield approximation which is the main object of this section.

(i) Gibbs’s density of states function and Boltzmann entropy.

Consider first the simple case of $`N`$ identical point particles of mass $`m`$ in a volume $`V`$, with perfectly reflecting massless walls in an external gravitational potential333To avoid any confusion we use here a tilde over all quantities with the same symbols and the same meanings as in self-gravitating systems. $`\stackrel{~}{U}(𝐫)`$. The total energy $`\stackrel{~}{H}`$, the sum of kinetic $`\stackrel{~}{H}_K`$ and potential energy $`\stackrel{~}{H}_P`$ of the system is a constant of motion. Let $`\stackrel{~}{E}`$ be its value. The total energy $`\stackrel{~}{H}`$ is a function of the coordinates of each particle in phase space $`(𝐫_i,𝐩_i)`$, $`(i=1,2,\mathrm{},N)`$:

$$\stackrel{~}{H}=\stackrel{~}{H}_K+\stackrel{~}{H}_P=\underset{i}{}\frac{1}{2m}p_i^2+\underset{i}{}m\stackrel{~}{U}(𝐫_i)=\underset{i}{}\stackrel{~}{E}_i=\stackrel{~}{E}.$$
(2.1)
The total energy $`\stackrel{~}{H}`$ is here equal to the sum of the energies of individual particles $`\stackrel{~}{E}_i`$. The density of states function $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E},V,N,m)`$ or in short $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})`$ is the sum of all possible states444The Gibbs density fo states function has dimensions $`[\stackrel{~}{\mathrm{\Omega }}]=erg^{3N1}sec^{3N}`$. Thus $`\stackrel{~}{\mathrm{\Omega }}`$ like later $`\mathrm{\Omega }`$ is defined up to a constant that depends on the units chosen. of the system with energy $`\stackrel{~}{E}`$ divided by $`N!`$:

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})=\frac{1}{N!}\delta (\stackrel{~}{E}\stackrel{~}{H})\underset{i}{}d\omega _i\mathrm{where}d\omega _i=d^3r_id^3p_i(\mathrm{}<𝐩_i<+\mathrm{});$$
(2.2)
$`\delta (\stackrel{~}{E}\stackrel{~}{H})\underset{i}{}d\omega _i`$ is known as the microcanonical distribution  and the basic justification for this distribution is the ergodic hypothesis. The Gibbs entropy $`\stackrel{~}{S}_g`$ is the logarithm of $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})`$ divided by Boltzmann’s constant $`k`$. Here we take $`k=1`$ which amounts mainly to measure the temperature in units of energy. Thus

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})=e^{\stackrel{~}{S}_g}.$$
(2.3)
To evaluate the Gibbs entropy we first replace $`\delta (\stackrel{~}{E}\stackrel{~}{H})`$ by a Bromwich integral of $`1`$ (its inverse Laplace transform) in brackets:

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})=\frac{1}{N!}\left[\frac{1}{2\pi i}_{bi\mathrm{}}^{b+i\mathrm{}}e^{\beta (\stackrel{~}{E}\underset{i}{}\stackrel{~}{E}_i)}𝑑\beta \right]\underset{i}{}d\omega _i.$$
(2.4)
This expression for $`\stackrel{~}{\mathrm{\Omega }}`$ can be rewritten after a slight rearrangement of terms as

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})=\frac{1}{2\pi i}_{bi\mathrm{}}^{b+i\mathrm{}}\left[\frac{1}{N!}\underset{i}{}e^{\beta \stackrel{~}{E}_i}d\omega _i\right]e^{\beta \stackrel{~}{E}}𝑑\beta \frac{1}{2\pi i}_{bi\mathrm{}}^{b+i\mathrm{}}\stackrel{~}{\mathrm{\Psi }}e^{\beta \stackrel{~}{E}}𝑑\beta ,$$
(2.5)
$`\stackrel{~}{\mathrm{\Psi }}`$, the quantity between brackets, is easily seen to reduce to

$$\stackrel{~}{\mathrm{\Psi }}=\frac{1}{N!}\left(e^{\beta \stackrel{~}{E}_0}𝑑\omega \right)^N\mathrm{where}\stackrel{~}{E}_0=\frac{1}{2m}p^2+m\stackrel{~}{U}(𝐫)\mathrm{and}d\omega =d^3rd^3p.$$
(2.6)
$`\stackrel{~}{E}_0`$ is the energy of of one particle with coordinates $`(𝐫,𝐩)`$. We can of course also write $`\stackrel{~}{\mathrm{\Psi }}`$ as follows

$$\stackrel{~}{\mathrm{\Psi }}=\underset{N^{}=0}{\overset{N^{}=\mathrm{}}{}}\frac{1}{N^{}!}\left(e^{\beta \stackrel{~}{E}_0}𝑑\omega \right)^N^{}\delta _{N^{}N}.$$
(2.7)
We then approximate $`\delta _{N^{}N}`$ by $`\delta (N^{}N)`$ since $`N>>1`$ and replace immediately $`\delta (N^{}N)`$ by a Bromwich integral like we did with $`\delta (\stackrel{~}{E}\stackrel{~}{H})`$:

$$\stackrel{~}{\mathrm{\Psi }}\frac{1}{2\pi i}_{ai\mathrm{}}^{a+i\mathrm{}}\underset{N^{}=0}{\overset{N^{}=\mathrm{}}{}}\frac{1}{N^{}!}\left(e^{\beta \stackrel{~}{E}_0}𝑑\omega \right)^N^{}e^{\alpha (N^{}N)}d\alpha .$$
(2.8)
$`\stackrel{~}{\mathrm{\Psi }}`$ can now be rewritten in this form

$$\stackrel{~}{\mathrm{\Psi }}\frac{1}{2\pi i}_{ai\mathrm{}}^{a+i\mathrm{}}e^{\left(\alpha N+{\scriptscriptstyle \stackrel{~}{f}𝑑\omega }\right)}𝑑\alpha \stackrel{~}{f}=e^{\alpha \beta \stackrel{~}{E}_0}.$$
(2.9)
and if we substitute this expression for $`\stackrel{~}{\mathrm{\Psi }}`$ back into (2.5) we obtain

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})\frac{1}{(2\pi i)^2}_{ai\mathrm{}}^{a+i\mathrm{}}_{bi\mathrm{}}^{b+i\mathrm{}}e^{\stackrel{~}{\sigma }}𝑑\alpha 𝑑\beta \stackrel{~}{\sigma }=\alpha N+\beta \stackrel{~}{E}+\stackrel{~}{f}𝑑\omega .$$
(2.10)
We shall now evaluate $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})`$ by a steepest descent technique  as follows. We first look for the extremum of $`\stackrel{~}{\sigma }(\alpha ,\beta )`$; this is obtained from the following equalities:

$$\frac{\stackrel{~}{\sigma }}{\alpha }=N+\stackrel{~}{f}𝑑\omega =0,\frac{\stackrel{~}{\sigma }}{\beta }=\stackrel{~}{E}\stackrel{~}{E}_0\stackrel{~}{f}𝑑\omega =0.$$
(2.11)
The equations define a point $`(\stackrel{~}{\alpha }_e,\stackrel{~}{\beta }_e)`$ in $`(\alpha ,\beta )`$ space in terms of $`\stackrel{~}{E}`$ and $`N`$ and a corresponding equilibrium value $`\stackrel{~}{f}_e`$ for $`\stackrel{~}{f}`$:

$$\stackrel{~}{f}_e=e^{\stackrel{~}{\alpha }_e\stackrel{~}{\beta }_e\stackrel{~}{E}_0}=e^{\stackrel{~}{\alpha }_e\stackrel{~}{\beta }_e[\frac{1}{2m}p^2+m\stackrel{~}{U}(𝐫)]}.$$
(2.12)
$`\stackrel{~}{f}_e`$ is the Boltzmann distribution of energies in the external field $`\stackrel{~}{U}(𝐫)`$ calculated at the point of extremum. The extremal value $`\stackrel{~}{\sigma }_e`$ of $`\stackrel{~}{\sigma }`$ is the Boltzmann entropy $`\stackrel{~}{S}`$ and it is easily seen that

$$\stackrel{~}{S}=\stackrel{~}{\sigma }_e=\stackrel{~}{f}_e\mathrm{ln}(\stackrel{~}{f}_e)𝑑\omega +N.$$
(2.13)
We may thus write (2.10) as follows introducing $`\delta ^2\stackrel{~}{\sigma }`$ which represents the sum of all the terms of order higher than one in a Taylor expansion of $`\stackrel{~}{\sigma }`$ near the extremum $`\stackrel{~}{\sigma }_e`$:

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})=e^{\stackrel{~}{S}_g}e^{\stackrel{~}{S}}\frac{1}{(2\pi i)^2}_{ai\mathrm{}}^{a+i\mathrm{}}_{bi\mathrm{}}^{b+i\mathrm{}}e^{\delta ^2\stackrel{~}{\sigma }}𝑑\alpha 𝑑\beta .$$
(2.14)
Since the exponent $`\delta ^2\stackrel{~}{\sigma }`$ is of order $`N`$ the integrant in (2.14) is very steep and a good approximation is obtained by limiting the integrations to terms of order two. The terms of order two are readily found; with $`\delta \alpha =\alpha \stackrel{~}{\alpha }_e`$ and $`\delta \beta =\beta \stackrel{~}{\beta }_e`$ we may write

$`\delta ^2\stackrel{~}{\sigma }`$ $``$ $`{\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{^2\stackrel{~}{\sigma }}{\alpha ^2}}\right)_e(\delta \alpha )^2+2\left({\displaystyle \frac{^2\stackrel{~}{\sigma }}{\alpha \beta }}\right)_e\delta \alpha \delta \beta +\left({\displaystyle \frac{^2\stackrel{~}{\sigma }}{\beta ^2}}\right)_e(\delta \beta )^2\right],`$ (2.15)
$`=`$ $`{\displaystyle \frac{1}{2}}\left[N(\delta \alpha )^22E\delta \alpha \delta \beta +N\overline{\stackrel{~}{E}_0^2}(\delta \beta )^2\right]`$
$`=`$ $`{\displaystyle \frac{1}{2}}N\left[(\delta \alpha \overline{\stackrel{~}{E}_0}\delta \beta )^2+\overline{(\stackrel{~}{E}_0\overline{\stackrel{~}{E}_0})^2}\delta \beta ^2\right],`$
in this formula mean values calculated in the phase space of a particle are denote with an overbar; for instance:

$$\overline{\stackrel{~}{E}_0}=\frac{1}{N}\stackrel{~}{E}_0\stackrel{~}{f}_e𝑑\omega $$
(2.16)
We can now fix the position of the imaginary axis for the complex variables $`(\alpha ,\beta )`$ by introducing a new set of variables which vary between $`\pm \mathrm{}`$:

$$\stackrel{~}{\alpha }^{}=i(\delta \alpha \overline{\stackrel{~}{E}_0}\delta \beta )\mathrm{and}\stackrel{~}{\beta }^{}=i\delta \beta .$$
(2.17)
This gives to (2.10) the following form:

$$\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})=e^{\stackrel{~}{S}_g}e^{\stackrel{~}{S}}\frac{1}{(2\pi )^2}_{\mathrm{}}^+\mathrm{}e^{\frac{N}{2}\stackrel{~}{\alpha }^2}𝑑\stackrel{~}{\alpha }^{}_{\mathrm{}}^+\mathrm{}e^{\frac{N}{2}\overline{(\stackrel{~}{E}_0\overline{\stackrel{~}{E}_0})^2}\stackrel{~}{\beta }^2}𝑑\stackrel{~}{\beta }^{}.$$
(2.18)
The integrals are easily found and the end product is this:

$$\stackrel{~}{S}_g\stackrel{~}{S}\mathrm{ln}(N)\mathrm{ln}(2\pi )\frac{1}{2}\mathrm{ln}\left[\overline{(\stackrel{~}{E}_0\overline{\stackrel{~}{E}_0})^2}\right].$$
(2.19)
This shows incidentally that $`\stackrel{~}{S}`$ is a maximum of $`\stackrel{~}{\sigma }`$ for any $`\stackrel{~}{U}`$; in other words the system is stable and $`\stackrel{~}{\mathrm{\Omega }}(\stackrel{~}{E})`$ is convergent. Moreover, since $`\stackrel{~}{S}`$ is of order $`N`$ and $`N>>1`$ the last three terms are clearly negligible compared to $`\stackrel{~}{S}`$ so that

$$\stackrel{~}{S}_g\stackrel{~}{S}.$$
(2.20)
This is the main result of this first subsection.

(ii) The Gibbs density of states function for self-gravitating systems.

If the same system of $`N`$ particles of mass $`m`$ in a volume $`V`$ are interacting gravitationally, the total conserved energy is

$$H=H_K+H_P=\underset{i}{}\frac{1}{2m}p_i^2\frac{1}{2}\underset{ij}{}\frac{Gm^2}{|𝐫_i𝐫_j|}=E.$$
(2.21)
$`G`$ is the gravitational constant. The density of states function $`\mathrm{\Omega }(E)`$ is the same as (2.2) without tildes:

$$\mathrm{\Omega }(E)=\frac{1}{N!}\delta (EH)\underset{i}{}d\omega _i.$$
(2.22)
As is well known, $`\mathrm{\Omega }`$ diverges. This can easily be seen as follows555Here we follow Padmanabhan .. Take equation (2.22), replace again $`\delta (EH)`$ by a Bromwich integral of $`1`$ with variable $`EH=E(H_K+H_P)`$; in this way we obtain an expression similar to (2.4) :

$$\mathrm{\Omega }(E)=\frac{1}{N!}\left[\frac{1}{2\pi i}_{bi\mathrm{}}^{b+i\mathrm{}}e^{\beta (EH_KH_P)}𝑑\beta \right]\underset{i}{}d\omega _i.$$
(2.23)
Integrate $`e^{\beta H_K}`$ over momentum space which is easy to do; this gives

$$\mathrm{\Omega }(E)=\frac{(2\pi m)^{\frac{3}{2}N}}{N!}_V\left[\frac{1}{2\pi i}_{bi\mathrm{}}^{b+i\mathrm{}}\frac{e^{\beta (EH_P)}}{\beta ^{\frac{3}{2}N}}𝑑\beta \right]\underset{i}{}d^3r_i.$$
(2.24)
Then use Table 15.2 in Arfkin  to calculate the Bromwich integral in the square brackets and replace $`H_P`$ by its definition in (2.21):

$$\mathrm{\Omega }(E)=\frac{(2\pi m)^{\frac{3}{2}N}}{N!(\frac{3}{2}N1)!}_V\left(E+\frac{1}{2}\underset{ij}{}\frac{Gm^2}{|𝐫_i𝐫_j|}\right)^{\frac{3N}{2}1}\underset{i}{}d^3r_i,$$
(2.25)
and finally change to new variables $`𝐫=𝐫_1𝐫_2`$ and $`𝐫_u`$ with $`u,v=2,3,\mathrm{},N`$ in terms of which

$$\mathrm{\Omega }(E)=\frac{(2\pi m)^{\frac{3}{2}N}}{N!(\frac{3}{2}N1)!}_V\underset{u}{}d^3r_u,$$
(2.26)
in this equality,

$$=_V\left(E+\frac{Gm^2}{r}+\underset{u=3}{\overset{u=N}{}}\frac{Gm^2}{|𝐫+𝐫_2𝐫_u|}+\frac{1}{2}\underset{uv}{}\frac{Gm^2}{|𝐫_u𝐫_v|}\right)^{\frac{3N}{2}1}d^3r.$$
(2.27)
We now see that for $`r0`$, $`r^{3\frac{3N}{2}}drr^{4\frac{3N}{2}}`$. So for $`N3`$ and $`r0`$, $`\mathrm{\Omega }\mathrm{}`$. Considerations of simple models show that the system with a subset of particles closely bound together by gravity and the remainder banging around with the high energy released has a large phase space volume associated with it . Such a state can presumably be reached given sufficient time for the system to evolve. A system of point particles is however not realistic. In more realistic models, one make use of small hard spheres or fermions or remove the divergence of the potential with artificial cutoffs. This makes $`\mathrm{\Omega }(E)`$ converge.

It is important to realize that once $`\mathrm{\Omega }`$ converges by say taking hard spheres, the stable state is not necessarily a dense core with a dilute halo; other configurations may be more stable. For a sufficiently small cut off at least two particles will be very close together since such states would dominate were the hard spheres reduced to points. Other configurations may only be more stable if the cut offs are big enough.

Models like these have been reanalyzed recently by Chavanis  who reviews earlier interesting works with various cutoffs. Neither of the models give however a direct evaluation of $`\mathrm{\Omega }(E)`$.

(iii) The meanfield approximation.

A meanfield theory is a most common approximation when short distance effects are negligible on the timescales considered. A steepest descent calculation gives a meaningful statistical thermodynamics, equilibrium configurations and stability conditions.

The following calculation is kept simple by not introducing any sort of cut off while embedding divergences into terms that would normally converge if there was a cutoff. These quantities when finite have no role in the results that interest us.

One should however be aware that the maximum of the new integrand $`e^\sigma `$ is not unique nor is it necessarily the highest. A second maximum may be much higher if the hard spheres are small enough. However, the calculation of that maximum is precisely brushed under the rug. A correct calculation is far too difficult. Our entropy will be associated with what is often referred to as a “local” maximum. This local maximum may not be the highest but it may be the more important one on time scales relevant in astronomy.

The calculation of $`\mathrm{\Omega }(E)`$ in the meanfield approximation is slightly intricate. We give here a heuristic derivation that will make the answer plausible. A complete and rigorous derivation that follows paper  is given in Appendix A.

Let us ask how (2.10) should be modified when $`\stackrel{~}{U}(𝐫)`$ is replaced by the gravitational field of the particles themselves? The field is now variable with $`\mathrm{}^3`$ degrees of liberty. The correct $`\mathrm{\Omega }`$ will have a functional integration over the space of all fields instead of a discrete integration over the $`𝐫_i`$-space. $`\stackrel{~}{\sigma }`$ needs a correction also that must depend on the field because the mean value of the energy in equation (2.11) is not the correct expression for a a self-gravitating system: the potential energy is counted twice. On the other hand the correction of $`\stackrel{~}{\sigma }`$ should not depend on $`\alpha `$ or $`\beta `$ since $`\stackrel{~}{S}`$ must also be the correct expression for the Boltzmann entropy whether the field is a self-field or an external one; the local distribution does not know what the origin of the field is.

With these remarks in mind, we now write what is the expression for the Gibbs density of states function $`\mathrm{\Omega }`$ in the “meanfield approximation” derived from (2.23) in Appendix A. First we must use a field variable $`W`$ slightly different from and more convenient than $`\stackrel{~}{U}`$. We replace $`\stackrel{~}{U}`$ by $`\beta ^{1/2}W(𝐫)`$. In terms of $`W(𝐫)`$

$$\mathrm{\Omega }(E)=\frac{1}{(2\pi i)^2}_{\mathrm{}}^+\mathrm{}𝒟W_{ai\mathrm{}}^{a+i\mathrm{}}_{bi\mathrm{}}^{b+i\mathrm{}}e^\sigma 𝑑\alpha 𝑑\beta ,$$
(2.28)
in this expression $``$ is a constant which actually diverges if there is no short distant cutoff (see appendix A). $`𝒟W`$ is the element of volume in function space (see also appendix A) while

$$\sigma [\alpha ,\beta ,W(𝐫)]=\alpha N+\beta E+e^{\alpha \beta [\frac{1}{2m}p^2+m\beta ^{1/2}W(𝐫)]}𝑑\omega +\frac{1}{8\pi G}_{\mathrm{}}^+\mathrm{}W\mathrm{\Delta }Wd^3r.$$
(2.29)
$`\mathrm{\Delta }`$ is the Laplacian. $`\sigma `$ differs from $`\stackrel{~}{\sigma }`$ in (2.10) mainly by the term $`W\mathrm{\Delta }W`$. This addition provides the correct mean value of the total energy as well as Poisson’s equations for equilibrium configurations as we shall now see. And it does not depend on $`\alpha ,\beta `$.

(iv) A steepest descent evaluation.

We start as in the case with an external field by calculating the extremal values of $`\sigma [\alpha ,\beta ,W(𝐫)]`$; these are associated with equilibrium configurations.

The equilibrium value $`W_e(𝐫)`$ of $`W(𝐫)`$ or rather $`U(𝐫)=\beta _e^{1/2}W_e(𝐫)`$ is defined by

$$\frac{\delta \sigma }{\delta W(𝐫)}=0\mathrm{\Delta }U=4\pi G\rho \mathrm{with}\rho (𝐫)=m_{\mathrm{}}^+\mathrm{}e^{\alpha _e\beta _e[\frac{1}{2m}p^2+mU(𝐫)]}d^3p,$$
(2.30)
$`\rho `$ is the mean mass density of the particles in a meanfield $`U`$ with a Boltzmann distribution

$$f_e=e^{\alpha _e\beta _e[\frac{1}{2m}p^2+mU(𝐫)]}=e^{\alpha _e\beta _eE_0},$$
(2.31)
$`E_0`$ is here the energy of one particle in the meanfield $`U`$ itself, a solution of Poisson’s equation (2.30). The equilibrium values $`(\alpha _e,\beta _e)`$ of $`(\alpha ,\beta )`$ are fixed by equations similar to (2.11) but $`E`$ is not quite the same as $`\stackrel{~}{E}`$:

$$\frac{\sigma }{\alpha }=0N=f_e𝑑\omega \mathrm{but}\frac{\sigma }{\beta }=0E=\left[\frac{1}{2m}p^2+\frac{1}{2}mU(𝐫)\right]f_e𝑑\omega .$$
(2.32)
The maximum $`\sigma _e`$ of $`\sigma `$, which is approximately the maximum of the integral when $`N>>1`$ is equal to the Boltzmann entropy. It is indeed easy to see that (2.31)(2.32) provide the same expression for the extremal value of $`\sigma `$ as (2.29):

$$S(E,N,V)\sigma _e=f_e\mathrm{ln}(f_e)𝑑\omega +N.$$
(2.33)
We may thus write (2.28) as follows introducing $`\delta ^2\sigma `$ which represents the sum of all the terms of order higher than one in a Taylor expansion near the extremum:

$$\mathrm{\Omega }(E)=e^{S_g}e^S\frac{1}{(2\pi i)^2}_{\mathrm{}}^+\mathrm{}𝒟W_{ai\mathrm{}}^{a+i\mathrm{}}_{bi\mathrm{}}^{b+i\mathrm{}}e^{\delta ^2\sigma }𝑑\alpha 𝑑\beta .$$
(2.34)
Since, as before, the exponent is of order $`N`$ the exponent in the integral is very steep and a good approximation will be obtained by limiting the integration to terms of order two. The reader interested in the derivation of $`\delta ^2\sigma `$ and its explicit expression will find them both in Appendix B. There we show that $`\mathrm{\Omega }(E)`$ is approximately given by

$$\mathrm{\Omega }(E)e^{S\mathrm{ln}(N)\mathrm{ln}(2\pi b/\beta _e)}_{\mathrm{}}^+\mathrm{}e^{\delta ^2\sigma ^{}}𝒟(\delta W),$$
(2.35)
$`\delta ^2\sigma ^{}`$ is a non-local quadratic functional in $`\delta W`$ which must be negative for the integral to converge, something of the form:

$$\delta ^2\sigma ^{}=_{\mathrm{}}^+\mathrm{}\delta W𝒪(𝐫,𝐫^{})\delta W^{}d^3rd^3r^{}<0.$$
(2.36)
$`\delta W`$ can in principle be expanded in terms of a complete set of orthonormal eigenfunctions say $`\xi _a(𝐫)`$ of the linear operator $`𝒪(𝐫,𝐫^{})`$ with arbitrary discrete variables $`\delta Y^a`$. The quadratic form $`\delta ^2\sigma ^{}`$ is then replaced by an infinite sum of squares $`(\delta Y^a)^2`$ with coefficients that for convenience we write $`\frac{1}{2}\lambda _a^{}`$; thus

$$\delta ^2\sigma ^{}=\frac{1}{2}\underset{a=1}{\overset{a=\mathrm{}}{}}\lambda _a^{}(\delta Y^a)^2.$$
(2.37)
The $`\lambda _a^{}`$’s can be arranged in increasing order $`\lambda _1^{}\lambda _2^{}\lambda _3^{}\mathrm{}`$. Clearly the integral in (2.35) will be convergent if $`\lambda _1^{}>0`$; $`S`$ is then a maximum and $`\lambda _1^{}>0`$ is a condition of stability. If $`N>>1`$ and $`\delta ^2\sigma ^{}`$ is sufficiently small, $`S_gS`$ as is usually expected from a steepest descent calculation. We can in principle evaluate the limit of validity of this quasi-equality by calculating $`\delta ^2\sigma ^{}`$.

(v) Convergence of $`\mathrm{\Omega }(E)`$ in meanfield theory.

There are two factors that cause $`\mathrm{\Omega }(E)`$ to be zero if it converges: $``$ as we said diverges and the spectrum of eigenvalues $`\lambda _a^{}`$ is unbounded if the domain of definition of eigenfunctions is bounded. This is true under very general conditions that do not depend on the form of a finite domain . It is true in particular for isothermal spheres  as shown in appendix B. In other volumes we may restrict $`\delta W`$’s to functions that satisfy a Poisson equation $`\mathrm{\Delta }\delta \rho =4\pi G\delta W`$ and use these $`\delta \rho `$’s as variables; they are zero outside $`V`$. In that case, either $`\delta ^2\sigma ^{}>0`$, the Boltzmann entropy is not a maximum of $`\sigma `$ which tends to plus infinity and $`\mathrm{\Omega }(E)`$ is undefined in this approximation because $``$ in the denominator of (2.34) divergence. Or $`\delta ^2\sigma ^{}<0`$ and the Gibbs density of states function tends to zero! There is no way out of this difficulty in a meanfield steepest descent calculation666see also the end of Appendix B in this respect..

The culprit of the trouble is easily traced to the Laplacian operator whose inverse is the divergent $`1/r`$ potential . Some form of short distance cutoff should have the effect of replacing the Laplacian by a non local operator with a convergent set of negative eigenvalues $`\lambda _a^{}`$ and a convergent $``$.

Notice one important point: that the number of nodes of the eigenfunctions grows steadily as $`\lambda _a^{}`$ increases. When the spacing between the nodes in the finite domain is of the order of the short distance cutoff the Laplacian in $`\delta ^2\sigma ^{}`$ is no more a reasonable representation for a convergent $`\mathrm{\Omega }(E)`$. This means however that a short distance cutoff will changes the higher end of the spectrum and that $`\lambda _1^{}>0`$ is a perfectly good stability condition for the meanfield configurations with a Boltzmann distribution.

(vi) Local and global maxima.

Had we made a correct calculation with some specific cutoff not too big, we would most likely have found at least two maxima of $`\sigma `$, the one we have just found and another one associated with a highly concentrated core and entropy $`S_0`$ say. This $`S_0`$ would tend to infinity if the cutoff went to zero. A steepest descent calculation would have given

$$\mathrm{\Omega }(E)=e^{S_g}e^{S_0+S}$$
(2.38)
Since both $`S`$ and $`S_0`$ are of order $`N>>1`$ the Gibbs entropy, for most values of the energy $`E`$ and volume $`V`$ would be approximately equal to either $`S`$ or $`S_0`$. Such a transition appears in a Monte-Carlo calculation of $`\mathrm{\Omega }(E)`$ by de Vega and Sanchez .

In principle the global maximum is the equilibrium configuration after an infinite time. However, the system may find itself in the local maximum for a very long time indeed. This makes local maxima interesting.

(vii) Remarks about relativistic systems and other ensembles.

(a) Relativistic systems.

The Gibbs density of states function can and has been calculated in general relativity . The advantage of general relativity is that it is a mean field theory to start with. The disadvantage is that the formalism is far more complicated, arbitrariness of coordinates creates additional difficulties in functional integrations and moreover a density of energy is not defined. With so much trouble, one finds however that the total entropy is again given by (2.33) but of course $`f_e`$ is different from (2.31).

Other ensembles have not found wide applications in astronomy but they are of theoretical interest.

(b) The grand canonical ensemble.

In this ensemble $`N`$ and $`E`$ are not kept fixed. What is fixed is $`\alpha `$ and $`\beta `$. It follows that the Gibbs equivalent for this ensemble $`\mathrm{\Omega }_{GC}`$, called the grand partition function, is in the meanfield approximation given by expression (2.28) in which $`\frac{1}{2\pi i}e^{\alpha N}d\alpha `$ and $`\frac{1}{2\pi i}e^{\beta E}d\beta `$ integrations are removed from $`e^\sigma `$ because $`N`$ and $`E`$ are not fixed anymore; thus

$$\mathrm{\Omega }_{GC}=\frac{1}{}_{\mathrm{}}^+\mathrm{}e^{\sigma _{GC}}𝒟W\sigma _{GC}=\sigma +\alpha N\beta E.$$
(2.39)
$`\mathrm{\Omega }_{GC}`$ is thus a double Laplace transform of $`\mathrm{\Omega }(E)`$. Thermodynamic equilibrium is associated with the extremum $`\sigma _{GCe}`$ of the exponent $`\sigma _{GC}`$. This extremum is simply related to the entropy:

$$\sigma _{GCe}=S+\alpha N\beta E.$$
(2.40)
The thermodynamic potential has no conventional name. $`N`$ and $`E`$ are now mean values which are however related to $`\alpha ,\beta `$ by the same equations (2.32). In this ensemble like in the following ones equations (2.30), (2.31) and (2.32) always hold at the point where the exponent is extremum. But equations (2.32) have slightly different interpretations.

(c) The canonical ensemble.

In this ensemble the system has a fixed number of particles and the temperature or $`\beta `$ is kept fixed rather than the total energy. The Gibbs equivalent for this ensemble $`\mathrm{\Omega }_C`$ is called the partition function and it is obtained from (2.28) by removing the $`\frac{1}{2\pi i}e^{\beta E}d\beta `$ integration from $`e^\sigma `$. Thus

$$\mathrm{\Omega }_C=\frac{1}{2\pi i}_{ai\mathrm{}}^{a+i\mathrm{}}_{\mathrm{}}^+\mathrm{}e^{\sigma _C}𝑑\alpha 𝒟W\sigma _C=\sigma \beta E.$$
(2.41)
The corresponding thermodynamic equilibrium function which plays the role of $`S`$ and associated with the extremum $`\sigma _{Ce}`$ of the exponent $`\sigma _C`$ is

$$\sigma _{Ce}=S\beta E=\beta F.$$
(2.42)
The thermodynamic potential $`F`$ is known as the free energy.

(d) The grand microcanonical ensemble.

This fourth ensemble not often considered  is one in which the total energy is fixed but the particle number is allowed to fluctuate. This ensemble is physically plausible if there are pair creations of particles. Thus $`\alpha `$ and the total energy are fixed. The Gibbs equivalent $`\mathrm{\Omega }_{GMC}`$ is obtained by removing the $`\frac{1}{2\pi i}e^{\alpha N}d\alpha `$ integration from $`e^\sigma `$:

$$\mathrm{\Omega }_{GMC}=\frac{1}{2\pi i}_{bi\mathrm{}}^{b+i\mathrm{}}_{\mathrm{}}^+\mathrm{}e^{\sigma _{GMC}}𝑑\beta 𝒟W\sigma _{GMC}=\sigma +\alpha N.$$
(2.43)
Other ensembles have been studied like isothermal spheres under constant pressure.

(viii) Summary and conclusions.

Here now are the main points of this section.

(1) We started from the Gibbs density of states function $`\mathrm{\Omega }`$, equation (2.22), which we wrote in a meanfield approximation as equations (2.28) (2.29).

(2) We used a steepest descent method to calculate equilibrium configurations which are given by Poisson’s equation (2.30) with a Boltzmann distribution of particles in phase space (2.31). We also calculated the entropy (2.33) and got a meaningful expression if (2.34) is convergent. The condition of convergence is a condition of stability. The results are only valid if some sort of short distance cutoff is assumed, the details of which are however unimportant. The limit of validity of the whole approximation scheme can be evaluated if $`V`$ is specified. It is not clear however if the Boltzmann entropy is a global or a local maximum. That depends on the cutoff and it would show up in a correct calculation.

Lynden-Bell and Wood, in dealing with isothermal spheres proceeded in a simpler way. They started from the entropy defined by (2.33) with an unknown distribution function $`f`$ restricted by conditions (2.32) in which $`\beta ^{1/2}W=U`$ is a solution of Poisson’s equation (2.30) with $`\rho =mf𝑑\omega `$. This out of equilibrium entropy is then extremized ($`\delta S=0`$) with arbitrary $`\delta f`$’s subject to two restrictions $`\delta N=\delta E=0`$. The extremum of $`f`$ is the Boltzmann distribution function (2.31). Padmanabhan  used from the same method to find stability conditions. This method is simpler and well adapted to gaseous systems. It short-circuits many of the difficulties encountered above. Thus the statistical thermodynamics of Gibbs ends up with the same mathematics but it must be said again that “it has a particular beauty of its own \[and\] is applicable quite generally to every physical system” .

Section 1 dealt with the thermodynamic equilibrium of $`N`$ point particles of mass $`m`$ isolated in a volume $`V`$ with energy $`E`$ and entropy $`S`$. Section 2 deals with thermodynamic functions which give statistical thermodynamics its physical content. We also develop the thermodynamic stability theory of equilibrium configurations. The theory of stability in section 2 applies to a much wider class of systems than the one studied in this section.

## 3 Thermodynamic equilibrium and stability

The first subsection deals with the thermodynamics of the system of $`N`$ particles with gravitational interactions in a volume $`V`$. This is our generic case. The rest of section 2 is more general. Please notice that we shall be concerned by equilibrium configurations only; we drop therefore the index e from $`\alpha _e,\beta _e`$ and $`f_e`$ which is we have no reason to use anymore and write simply $`\alpha ,\beta `$ and $`f`$.

(i) Thermodynamic functions of self gravitating particles.

Besides $`E,N,V`$ and $`S(E,N,V)`$ there are other thermodynamic functions that give a physical content to the results. Consider in particular the derivatives of $`S`$ with respect to $`E,N,V`$. The derivatives are easily calculated by considering $`\sigma `$ defined in (2.29). This is a function of $`E,N,V`$ as well as of $`\alpha ,\beta ,W(𝐫)`$. But the entropy is an extremum of $`\sigma `$ for which $`\sigma /\alpha =\sigma /\beta =\delta \sigma /\delta W(𝐫)=0`$. Therefore the partial derivatives of S with respect to $`E,N,V`$ are the same as the partial derivatives of $`\sigma `$ with respect to $`E,N,V`$ keeping $`\alpha ,\beta ,W(𝐫)`$ fixed. The partial derivatives must of course be evaluated at the extremum. A quick look at (2.29) shows immediately that:

$$\frac{S}{E}=\beta ,\frac{S}{N}=\alpha ,\frac{S}{V}=\left(_{\mathrm{}}^+\mathrm{}fd^3p\right)_b=\frac{\rho _b}{m};$$
(3.1)
the index b means “on the boundary”777 There would be no mechanical equilibrium otherwise as the boundary must be an equipotential. Indeed, $`\rho _b`$ must be the same everywhere on the boundary as mechanical equilibrium requires: $`P\rho U=0`$ which implies, see (3.2), $`U\mathrm{ln}\rho `$ and $`U_e`$ must be the same everywhere on the boundary because the pressure force is normal to the surface.. It is worth noticing that the local pressure

$$P(𝐫)=_{\mathrm{}}^+\mathrm{}\frac{1}{3}m\dot{𝐫}^2fd^3p=\frac{1}{\beta }_{\mathrm{}}^+\mathrm{}fd^3p=\frac{\rho (𝐫)}{m\beta };\mathrm{thus}\frac{S}{V}=\beta P_b,$$
(3.2)
where $`P_b`$ is the pressure on the boundary. It follows from (3.1) and (3.2) that

$$dS=\beta dE\alpha dN+\beta P_bdV\mathrm{or}\frac{1}{\beta }dS=dE\frac{\alpha }{\beta }dN+P_bdV.$$
(3.3)
This differential expression, the most important one in statistical thermodynamics, shows that

$$\frac{1}{\beta }=T\mathrm{is}\mathrm{the}\mathrm{temperature},\mathrm{and}\frac{\alpha }{\beta }=\mu \mathrm{is}\mathrm{the}\mathrm{Gibbs}\mathrm{chemical}\mathrm{potential}.$$
(3.4)
$`T,\mu .P_b`$ are global thermodynamic quantities; they concern the whole system.

(ii) Conjugate thermodynamic functions.

Of particular interest in stability analysis are pairs of conjugate thermodynamic functions with respect to the entropy like $`(E,\beta =S/E)`$, $`(N,\alpha =S/N)`$ and $`(V,\beta P_b=S/V)`$. Pairs of conjugate parameters in other ensembles are easily obtained from (3.3). For instance in the grand canonical ensemble the pairs of conjugate parameters with respect to $`\sigma _{GCe}`$ are $`(\beta ,E)`$, $`(\alpha ,N)`$ and $`(V,\beta P_b)`$ since from (3.3) and (2.40) we see that

$$d\sigma _{GCe}=Ed\beta +Nd\alpha +\beta P_bdV.$$
(3.5)
In the canonical ensemble pairs of conjugate parameters associated with the free energy or rather $`\beta F`$ are $`(\beta ,E)`$ and $`(N,\alpha )`$ and $`(V,\beta P_b)`$ and for the grand microcanonical ensemble the pairs of conjugate thermodynamic functions with respect to $`\sigma _{GMCe}`$ are $`(E,\beta ),(\alpha ,N),(V,\beta P_b)`$. From here on the thermodynamics of self-gravitating systems can be developed along the royal path taken by Lynden-Bell and Wood  for studying isothermal spheres. Now we continue at a more general level.

(iii) Stability conditions for equilibrium states in general.

If the Gibbs density of states function in the mean field approximation is convergent and the steepest descent evaluation of $`\mathrm{\Omega }(E)`$ to order two is meaningful, equation (2.34) together with equation (2.37) suggests that with a proper choice of variables in $`\mathrm{\Omega }(E)`$, say $`x^a`$, in (2.28) the stability conditions would appear directly in their simplest form. Thus starting from any Gibbs density of states function of any ensemble in these appropriate variables we would write $`\mathrm{\Omega }(E)`$ in a form like this:

$$\mathrm{\Omega }(E)=\frac{1}{𝒞}e^{w(E;x^a)}\underset{a}{}dx^a,$$
(3.6)
where $`𝒞`$ is a constant that depends on the choice of the $`x^a`$’s. The extremum of $`w`$ which defines equilibrium configurations is then of the form

$$\frac{w}{x^a}=0\mathrm{whose}\mathrm{solutions}x^a=X^a(E);$$
(3.7)
there may be more than one solution for given $`E`$ and whatever other parameters there may be, say $`\zeta _q,(q=1,2,\mathrm{})`$. Near an extremum where $`x^a=X^a(E)`$, the expansion of $`w`$ would be of the following form:

$$w=w_e\frac{1}{2}\underset{a}{}\lambda _a[x^aX^a(E)]^2+O_3.$$
(3.8)
Equations (3.7) and (3.8) are similar to those obtained for testing equilibrium and stability of a mechanical system with potential energy $`w`$, variables $`x^a`$ and parameters $`(E,\zeta _q)`$ though the number of variables in this case is usually finite. The whole analysis as can be seen is in fact valid for self-gravitating systems or other systems with short range forces in thermal equilibrium.

The calculation of the stability limit $`\lambda _1`$ is not easy in general. Therefore the following criteria may be helpful.

(iv) A general criterion of stability.

(a) A useful identity

For definiteness we keep in mind the microcanonical ensemble and the $`w(=\sigma )`$ introduced in the previous section but keep our new notations that refer to more general systems. Consider the parameter888Any other parameter might do. $`E`$ and introduce its conjugate with respect to $`w`$ rather than with respect to $`w_e`$. This defines a sort of inverse “temperature out of equilibrium” say $`\stackrel{~}{\beta }(E;x^a)`$. The derivative of $`w`$ given by equation (3.8) with respect to $`E`$ keeping all $`x^a`$’s and $`\zeta _q`$’s fixed:

$$\stackrel{~}{\beta }(E;x^a)=\frac{w}{E}=\beta (E)+\underset{a}{}\lambda _a\frac{dX^a}{dE}[x^aX^a(E)]+𝒪_2,\beta (E)=\frac{S}{E}.$$
(3.9)
Clearly the equilibrium value of $`\stackrel{~}{\beta }`$ is $`\beta `$ and the second order terms is to be interpreted as fluctuations of $`\beta `$ near equilibrium. When $`w=\sigma `$, $`\beta `$ is the inverse temperature; here $`\beta `$ is some “generalized” inverse temperature.

Let us calculate first order derivatives of $`\stackrel{~}{\beta }`$ with respect to $`x^a`$ and $`E`$ at the point of extremum. At that point derivatives of $`𝒪_2`$ in (3.9) are at least of order 1 and will thus add nothing to the derivatives at the point of extremum. The second terms of (3.9) contribute only through derivatives of $`[x^aX^a(E)]`$ because this quantity is zero at the extremum. One thus obtain:

$$\left(\frac{\stackrel{~}{\beta }}{x^a}\right)_e=\lambda _a\left(\frac{X^a}{E}\right)\mathrm{no}\mathrm{summation}\mathrm{on}a,\left(\frac{\stackrel{~}{\beta }}{E}\right)_e=\frac{\beta }{E}\underset{a}{}\lambda _a\left(\frac{X^a}{E}\right)^2.$$
(3.10)
Extracting $`X^a/E`$ from the first equality and inserting into the second gives the following expression for the slope of the “linear series” of conjugate parameters $`\beta (E)`$ with respect to $`w_e`$:

$$\frac{\beta }{E}=\left(\frac{\stackrel{~}{\beta }}{E}\right)_e+\underset{a}{}\frac{(\stackrel{~}{\beta }/x^a)_e^2}{\lambda _a}.$$
(3.11)
Let us emphasize again that the calculation is valid for any pair of conjugate parameters with respect to any $`w`$, $`e^w`$ being the statistical weight of any thermodynamic system, including one with short range forces between the particles. It also applies to any pair of conjugate parameters with respect to (minus) the potential energy of a mechanical system in equilibrium. It is however extremely useful to keep in mind the concrete example of the $`N`$ gravitating particles of mass $`m`$ with energy $`E`$ in a volume $`V`$.

(b) Poincaré’s criteria of stability.

There exists a classical result of Poincaré  which is discussed in various treatises on stability, see for instance Jeans  Ledoux  or Lyttleton , and was applied to self-gravitating systems in thermal equilibrium for the first time by Lynden-Bell and Wood . It says the following. A change of stability may only occur where two or more linear series of equilibria like $`\beta (E)`$ for instance999Any linear series, not only of conjugate thermodynamic functions like $`\beta (E)`$. For instance, one of the $`X^a`$’s as a function of one of the $`\zeta _q`$’s or a $`P(V)`$ diagram. have one equilibrium configuration in common (“bifurcation points”) or where two or more series merge into each other (“turning points”). When this happens stable equilibria may turn into unstable ones; reciprocally, unstable ones may become either stable or more unstable. Thus a change of sign of $`\lambda _1`$ and perhaps of more eigenvalues or stability coefficients appear only at bifurcations and turning points.

The Poincaré method of linear series does not say whether a change of stability actually occurs or not. It has been pointed out  that under very general conditions the stability or the number of unstable modes (the number of negative $`\lambda _a`$’s) can be deduced from the topological properties of series of equilibria, i.e. from purely thermodynamic considerations as we shall now review.

(c) Simple eigenvalues.

Suppose the spectrum of eigenvalues is simple that is non-degenerate: $`\lambda _1<\lambda _2<\lambda _3\mathrm{}.`$. This assumption is far more than we shall ever need but it makes further explanations simpler. In practice only some of the smaller $`\lambda _a`$’s need to be distinct. Consider a line $`\beta (E)`$ with points that correspond to stable configurations ($`\lambda _1>0`$) and follow the line towards a limit of stability. As we approach that limit $`\lambda _10`$ the $`1/\lambda _1`$ term in (3.11) begins to dominate and

$$\frac{\beta }{E}\frac{(\stackrel{~}{\beta }/x^1)_e^2}{\lambda _1}+\mathrm{}\mathrm{provided}(\stackrel{~}{\beta }/x^1)_e0.$$
(3.12)
We shall come back later to what happens if $`(\stackrel{~}{\beta }/x^1)_e0`$. If equation (3.12) holds, the slope of $`\beta (E)`$ has the same sign as $`\lambda _1`$. As a result:

(1) Two linear series of conjugate thermodynamic functions merge into each other at a point where there is a vertical tangent. A change of stability occurs thus necessarily there, at a turning point.

(2) The stable branch is the one with a positive slope near the vertical tangent. Another way to state this property is as follows. Stability is lost where the branch turns counterclockwise. A stable branch never turns clockwise.

(3) Say $`\lambda _1<0`$ and $`\lambda _20`$ along the $`\beta (E)`$ line. There is thus already one mode of instability. If a further mode of instability shows up it must be at a new turning point where $`\lambda _20`$. Near the new turning point the dominant contribution to (3.11) comes from the $`1/\lambda _2`$ term and equation (3.12) holds with an indice 2 instead of an indice 1.

(4) Say $`\lambda _1<0<\lambda _2`$ and $`\beta (E)`$ turns clockwise. This means that $`\lambda _10`$ again and that unstable configurations turn into stable ones on the side with positive slopes.

(5) If a stable linear series spirals inwards counterclockwise, stability is lost at the first turning point. We give an example of this case in section 4.

(6) If a change of stability occurs at a point where $`(\stackrel{~}{\beta }/x^1)_e0`$ and $`\lambda _10`$ the previous considerations become obviously invalid. Notice however that a small perturbation of the potential energy easily removes the coincidence. Indeed, let us change $`w`$ to $`w+ϵEx^1`$ say with $`|ϵ|<<1`$. The equilibrium equations (3.7) become

$$\frac{w}{x^1}+ϵE=0(|ϵ|<<1)\mathrm{and}\frac{w}{x^a}=0(a>1).$$
(3.13)
while

$$\left(\frac{\stackrel{~}{\beta }}{x^1}\right)_e=\left(\frac{^2[w+ϵEx^1]}{Ex^1}\right)_e\stackrel{\lambda _10}{}ϵ0\frac{\beta }{E}\frac{ϵ^2}{\lambda _1}.$$
(3.14)
This shows that such coincidences will show up in the mathematics but not in experimental observations  nor are they likely to appear in numerical calculations and many interesting models in astronomy and astrophysics are solved numerically.

(7) Finally if we know that one configuration is stable or we know its degree of instability we shall know about the stability or degree of instability of every other configuration of the linear series.

(d) Bifurcations and degenerate spectra.

Bifurcations like degeneracies of the spectrum of stability coefficients come often from an excess of symmetry. A characteristic property of bifurcations in problems of physical interest is that they are “unstable”. Small perturbations of the potential energy as just indicated remove unstable bifurcations and transform them into turning points ; such bifurcations do not show in experiments or in numerical solutions. We shall give an example of an unstable bifurcation in section 4. The theory of bifurcations101010The problem of distinguishing stable from unstable bifurcation offers little difficulty with linear series depending on one parameter . The unfloding of bifurcations with more than one parameter and several variables is much more complicated. Thom  classified all structurally stable inequivalent unfoldings or as it is called “elementary catastrophes” with up to two variables and four parameters. Great progress was done in this subject since 1975. Catastrophe theory has not been very useful in physics. It did not contribute a thing that was not already known otherwise. It provided however a sound mathematical basis to the theory of bifurcations and to the classification of inequivalent ones.  is of some interest in mathematical physics .

In the rest of this section we deal only with simple spectra and turning points at vertical tangents.

(v) $`C_V<0`$ and the nonequivalence of ensembles.

(a) About $`C_V<0`$.

We can now see that in a thermally stable isolated system but near instability $`(\lambda _10)`$ the heat capacity $`C_V=E/T`$ is negative in stable configurations and turns positive in unstable configurations by going through zero. Indeed, according to (3.12),

$$C_V=\frac{E}{T}=\beta ^2\frac{E}{\beta }\lambda _1\beta ^2\left(\frac{\stackrel{~}{\beta }}{x^1}\right)_e^2<0.$$
(3.15)
The same will of course happen in an unstable system on the verge of acquiring a second mode of instability $`(\lambda _1<0<\lambda _20)`$. $`C_V`$ may thus also be negative in unstable systems. The sign of $`C_V`$ is therefore not a criterion of stability like in classical thermodynamics.

The fundamental reason for this difference is the fact that the system is no subsystem of an ensemble. This is made obvious by the following standard explanation of why small subsystems are stable if their heat capacity is positive.

Consider an isolated system of energy $`E`$ made up of two systems of energy $`E_1`$ and $`E_2=EE_1`$ of comparable size. Let $`S_1(E_1)`$ and $`S_2(E_2)`$ be their entropies. The sum of the entropies is the entropy of the system out of equilibrium but for simplicity we write

$$S(E)=S_1(E_1)+S_2(E_2).$$
(3.16)
Equilibrium is reached when

$$\delta S=\left[\left(\frac{dS_1}{dE_1}\right)_e\left(\frac{dS_2}{dE_2}\right)_e\right]\delta E_1=(\beta _1\beta _2)\delta E_1=0\beta _1=\beta _2\beta .$$
(3.17)
To check stability we calculate second order variations of $`S`$ which at equilibrium must be negative:

$$(\delta ^2S)_e=\frac{1}{2}\left[\left(\frac{d\beta _1}{dE_1}\right)_e+\left(\frac{d\beta _2}{dE_2}\right)_e\right](\delta E_1)^2=\frac{1}{2}\beta ^2\left(\frac{1}{C_{V1}}+\frac{1}{C_{V2}}\right)(\delta E_1)^2<0.$$
(3.18)
Since $`\delta E_1=C_{V1}\delta T_1=C_{V1}\delta \beta _1/\beta ^2`$, we can rewrite the last expression as follows:

$$\left(\delta ^2S\right)_e=\frac{1}{2}\left(1+\frac{C_{V1}}{C_{V2}}\right)C_{V1}\left(\frac{d\beta _1}{\beta }\right)^2<0.$$
(3.19)
Now if system 1 is a small subsystem, $`E_1E_2`$, also $`C_{V1}C_{V2}`$ and

$$(\delta ^2S)_e\frac{1}{2}C_{V1}\left(\frac{\delta \beta _1}{\beta }\right)^2<0.$$
(3.20)
The stability conditions (3.20) implies that a stable small sub-system must have a positive heat capacity as is well known.

(b) About the non-equivalence of ensembles.

If instead of keeping the system isolated we keep it at constant temperature, the role of $`S`$ is then played by the thermodynamic potential $`w_{Ce}`$ or $`\beta F`$, see equation (2.42), and $`(\beta ,E)`$ is a pair of conjugate parameters with respect to $`\beta F`$. But $`E(\beta )`$ is the same line as $`\beta (E)`$ drawn in coordinates rotated $`90^{}`$ clockwise. Vertical tangents become horizontal ones and reciprocally. Let $`\lambda _{aC}`$ be the Poincaré coefficients of stability for the canonical ensemble. Near instability $`\lambda _{1C}0`$. The analog of (3.12) is now

$$\frac{(E)}{\beta }\frac{\left[(\stackrel{~}{E})/x_C^1\right]_e^2}{\lambda _{1C}}>0C_V=\frac{E}{T}=\beta ^2\frac{(E)}{\beta }>0.$$
(3.21)
This shows that near instability a stable canonical ensemble has $`C_V>0`$. Nevertheless unstable configurations may also have $`C_V>0`$. For instance if $`\lambda _{1C}<0<\lambda _{2C}`$ and $`\lambda _{2C}0`$. Thus $`C_V>0`$ is not a criterion of stability for a canonical ensemble either. It is surely a sufficient condition of stability near and indeed far away from instability.

The canonical ensemble and the isolated system or microcanonical ensemble cannot become unstable for the same equilibrium configuration unless vertical and horizontal tangents appear at the same point. This difference between ensembles is referred to as the non-equivalence of ensembles and is a property of non-extensive systems like self-gravitating ones. In extensive systems, different ensembles are equivalent . Why does the canonical ensemble get unstable and the other does not is nicely explained in the Lynden-Bell and Wood paper .

(vi) Summary and conclusions.

Under general conditions, mainly that the spectrum of Poincaré coefficients of stability be simple, stability limits and the number of unstable modes can be found for all calculated configurations using only linear series of thermodynamic functions of equilibrium configurations provided we know whether one configuration is stable or what is its degree of instability. Conjugate thermodynamic functions are particularly valuable in this respect as changes of stability show up at vertical tangents and the slopes in their vicinity have simple interpretations.

The thermodynamic criterion has a number of limitations. One rarely calculate all the sequences of equilibrium and therefore some bifurcations may not show up because the branch points are missing. Thus equilibria might become unstable and the system might choose to be in a more stable state which has not been calculated. A complete answer needs a detailed analysis of the second order terms of the Gibbs density of states function like the one started in subsection $`2(i)`$. The much simpler version of Padmanabhan  is also useful.

A less severe limitation is that the method is mathematically not a full proof because we must assume that the spectrum of stability coefficients is simple and at the same time that conditions like (3.12) hold. The method may however be applied to numerical solutions with a reasonable degree of confidence because this type of mathematical singularity will not show up in numerical calculations nor indeed in experiments.

Another general limitation of any thermodynamic criterion of stability is that we learn little about the nature of instabilities, triggering mechanisms, and what becomes of stable states which evolve through a series of quasi-equilibria along the linear series up to and beyond the limit of instability.

We have seen that the sign of heat capacities is not a criterion of stability and that in microcanonical and canonical ensembles near instability heat capacities have opposite signs. These ensembles are not equivalent as in classical thermodynamics. The proofs of these properties have a great degree of generality. The next section deals with fluctuations.

## 4 Fluctuations

(i) Fluctuations in self-gravitating systems.

The thermodynamics of fluctuations in stable equilibrium configurations becomes interesting near instability. Fluctuations can put the system out of equilibrium. Calculations of fluctuations are particularly interesting when dealing with metastable states of systems like the system of point particles in an enclosure considered in section 1. The entropy of such a system may be a “local” maximum smaller than the absolute maximum. In that case there are two maxima and between them there must a saddle point. Such a saddle point must show up close to a limit of stability because there the linear series turns counterclockwise. Figure 1 illustrates the situation with a linear series $`\beta (E)`$ for isolated systems. Let $`\beta `$ be the inverse temperature of a stable system and $`\beta ^{}<\beta `$ the inverse temperature of the unstable system with the same energy $`E`$. If the mean quadratic fluctuation of temperature or rather of its inverse $`<(\delta \beta )^2>`$ is equal to or greater than the square of the gap $`(\beta \beta ^{})^2`$, the statistical weight of a real fluctuation bigger than $`(\beta \beta ^{})`$ becomes important; this may put the system in a state where $`w`$ can increase to a bigger local maximum or to the global maximum. That is the entropy will have a chance to increase which it will surely do and equilibrium will be lost. It is thus interesting to evaluate the statistical weight of the fluctuations. Fluctuations can only displace the limit of stability towards higher values of $`\beta `$ and the question is how much higher. This is the subject of this section. The theory which has the same degree of generality as the theory of stability of section 2, comes from works by Okamoto, Parentani and myself .

(ii) The probability of a fluctuation near instability.

Figure 1 is a good starting point. The point of marginal stability has coordinates $`(E_m,\beta _m)`$ and a vertical tangent. Consider near that point a stable configuration with coordinates $`(E,\beta )`$. The dominant contribution to $`\mathrm{\Omega }`$ near instability, see (3.6), must come from $`\delta x^1`$ in $`w`$ since $`\mathrm{\Omega }\mathrm{}`$ when $`\lambda _10`$. Let us therefore integrate over all variables except $`x^1`$:

$$\mathrm{\Omega }(E)\frac{1}{𝒞}_{\mathrm{}}^+\mathrm{}e^{S\frac{1}{2}{\scriptscriptstyle \lambda _a(\delta x^a)^2}}\underset{a=1}{\overset{a=\mathrm{}}{}}dx^a\frac{e^S}{𝒞}\underset{a=2}{\overset{a=\mathrm{}}{}}\left(\frac{2\pi }{\lambda _a}\right)^{\frac{1}{2}}_{\mathrm{}}^+\mathrm{}e^{\frac{1}{2}\lambda _1(\delta x^1)^2}d(\delta x^1).$$
(4.1)
As long as the exponent is sufficiently steep (thus not too close to where $`\lambda _1=0`$) we may integrate over $`\delta x^1`$ between $`\pm \mathrm{}`$ as indicated.

Notice that $`\delta x^1=x^1X^1`$ where $`X^1`$ is the equilibrium value of $`x^1`$ at point $`(E,\beta )`$ and $`x^1`$ itself is some value out of equilibrium with the same energy $`E`$. We have a relation between $`\delta x^1`$ and the fluctuation of temperature $`\delta \beta =\stackrel{~}{\beta }(x^1;E)\beta `$ which to first order is simply

$$\delta \beta =\stackrel{~}{\beta }\beta \left(\frac{\stackrel{~}{\beta }}{x^1}\right)_e\delta x^1.$$
(4.2)
Replacing $`(\stackrel{~}{\beta }/x^1)_e`$ in terms of $`E/\beta `$ using the second of equations (3.12),

$$\delta x^1\frac{\delta \beta }{(\stackrel{~}{\beta }/x^1)_e}\pm \left(\frac{1}{\lambda _1}\frac{E}{\beta }\right)^{1/2}(\stackrel{~}{\beta }\beta ).$$
(4.3)
We use both these equalities to replace $`\delta x^1`$ in the integral of (4.1), the first equality in $`d(\delta x^1)`$ the second in the exponent, and obtain

$$\mathrm{\Omega }(E)\frac{e^S}{𝒞(\stackrel{~}{\beta }/x^1)_e}\underset{a=2}{\overset{a=\mathrm{}}{}}\left(\frac{2\pi }{\lambda _a}\right)^{\frac{1}{2}}_{\mathrm{}}^+\mathrm{}e^{\frac{1}{2}\frac{E}{\beta }(\stackrel{~}{\beta }\beta )^2}d(\stackrel{~}{\beta }\beta ).$$
(4.4)
We conclude that the statistical weight of a fluctuation of $`\beta `$ to between $`\stackrel{~}{\beta }`$ and $`\stackrel{~}{\beta }+d\stackrel{~}{\beta }`$ is proportional to the integrand of (4.4). Thus the probability $`d𝒫`$ of the fluctuation is given by the normalized expression:

$$d𝒫=\frac{1}{\sqrt{\pi }}e^{t^2}dt\mathrm{with}t^2=\frac{1}{2}\frac{E}{\beta }(\stackrel{~}{\beta }\beta )^2.$$
(4.5)
The mean quadratic fluctuation of temperature is thus given by

$$_{\mathrm{}}^+\mathrm{}\left(\frac{\delta \beta }{\beta }\right)^2𝑑𝒫=<\left(\frac{\delta \beta }{\beta }\right)^2>=<\left(\frac{\delta T}{T}\right)^2>=\frac{1}{(C_V)}>0.$$
(4.6)
Here we used (3.21) to get the expression on the right hand side. The formula is like that of Landau and Lifshitz  for mean quadratic fluctuations of temperature in a small subsystem; it has by necessity the opposite sign for $`C_V`$.

Fluctuations of other mean values near instability in different ensembles can be derived in a similar way.

The formula (4.5) may be made more specific because we are near the turning point where $`(E/\beta )_m=0`$. $`(E/\beta )`$ is thus small and to the lowest order in $`(\beta \beta _m)`$, assuming $`(^2E/\beta ^2)_m0`$,

$$\frac{E}{\beta }\left(\frac{^2E}{\beta ^2}\right)_m(\beta \beta _m)=\left(\beta ^3\frac{^2E}{\beta ^2}\right)_m\left(\frac{\beta \beta _m}{\beta _m}\right)\beta _m^2=\left(T\frac{C_V}{T}\right)_m\left(\frac{\beta \beta _m}{\beta _m}\right)\beta _m^2.$$
(4.7)
$`t^2`$ defined in (4.5) can also be written

$$t^2=\frac{1}{2}\left(T\frac{C_V}{T}\right)_m\left(\frac{\beta \beta _m}{\beta _m}\right)\left(\frac{\stackrel{~}{\beta }\beta }{\beta _m}\right)^2,$$
(4.8)
and the mean quadratic fluctuation satisfy the following equation:

$$\left(T\frac{C_V}{T}\right)_m\left(\frac{\beta \beta _m}{\beta _m}\right)<\left(\frac{\stackrel{~}{\beta }\beta }{\beta _m}\right)^2>=1.$$
(4.9)
If $`(^2E/\beta ^2)_m=0`$ a higher even derivative must be different from zero because $`E`$ is a minimum and the final results will not be very different; there will be higher powers of $`\beta \beta _m`$.

(iii) Stability limits induced by fluctuations.

We now define a new limit of stability, the point $`(\beta _l,E_l)`$ on the linear series in figure 1 where the mean quadratic fluctuation equals the square of the difference between $`\beta `$ and $`\beta ^{}`$ for reasons explained at the beginning of this section. Near the point of marginal stability $`(\beta _m,E_m)`$ we have approximately $`(\beta _l\beta _l^{})2(\beta _l\beta _m)`$; thus $`\beta _l`$ is defined by the condition that

$$<(\stackrel{~}{\beta }\beta _l)^2>4(\beta _l\beta _m)^2.$$
(4.10)
We shall obtain $`\beta _l`$ by inserting this value of the mean square fluctuations into (4.9); this gives

$$\left(4T\frac{C_V}{T}\right)_m\left(\frac{\beta _l\beta _m}{\beta _m}\right)^3=1\beta _l=\beta _m\left[1+\left(4T\frac{C_V}{T}\right)_m^{1/3}\right].$$
(4.11)
Notice that $`C_V`$ is of order $`N`$ and therefore the change from $`\beta _m`$ to $`\beta _l`$ is of order $`N^{1/3}`$ which is not necessarily very small. The corresponding value of $`E_l`$ is readily found by expanding $`E(\beta )`$ in a Taylor series to order two near $`\beta _m`$:

$$E_lE_m+\frac{1}{2}\left(\frac{^2E}{\beta ^2}\right)_m(\beta _l\beta _m)^2E_l=E_m+\frac{1}{8}\left(4T\frac{C_V}{T}\right)_m^{1/3}\beta _m^1.$$
(4.12)
To have some idea on how sharply defined the point $`(\beta _l,E_l)`$ is, consider an equilibrium configuration with $`\beta >\beta _l`$ and evaluate the weight of the probability for a fluctuation to induce instability, i.e. calculate $`e^{t^2}`$ for

$$\stackrel{~}{\beta }\beta 2(\beta \beta _m).$$
(4.13)
If we insert this value for $`\stackrel{~}{\beta }\beta `$ into (4.8), taking account of (4.11) we find that

$$t^2=\frac{1}{2}\left(\frac{\beta \beta _m}{\beta _l\beta _m}\right)^3.$$
(4.14)
For $`\beta \beta _m=\beta _l\beta _m`$, $`e^{t^2}0.6`$ but if $`\beta \beta _m`$ is twice as big as $`\beta _l\beta _m`$, $`e^{t^2}0.02`$ and if equal to three times that difference, $`e^{t^2}10^6`$. The exponential fall is relatively sharp and the new limit of stability $`(E_l,\beta _l)`$ induced by fluctuations is not too badly localized.

(iv) Summary and comments.

Fluctuations in non-extensive systems like self-gravitating ones have a peculiar behavior which has no analogue in classical thermodynamics where smooth turning points do not exist and ensembles are equivalent. Fluctuations can induce instability at lower temperatures (higher $`\beta `$’s) than the theoretical limit $`\beta _m`$. The relative change is of order $`N^{1/3}`$. $`N=10^3`$ is often considered in numerical calculations to study globular cluster models. Instability limits found in these models may be quite different from more realistic ones with $`N=10^5`$ or higher. Moreover, since $`N^{1/3}`$ is not necessarily a big number the stability limit may be sensitive to the values of $`(4TC_V/T)_m`$.

Fluctuation theory is useful when slow evolution occurs towards instability as we shall see in section 5.

## 5 Examples

In this section we give two examples in which thermodynamics of self-gravitating systems is used. The first one is that of isothermal spheres. Equilibrium configurations were obtained numerically. In this case linear series of conjugate thermodynamic functions correctly give stability limits. The calculation of fluctuations give useful additional information. The second example is that of liquid ellipsoids. Here equilibrium configurations are defined by quadratures and the multiple linear series are replete with bifurcations. The bifurcations are however unstable and may be lifted with small perturbations of the potential energy. Other examples are mentioned in less detail. In some of these examples, the stability theory of section 3 has been used and in others it has not but might have saved complicated calculations to obtain the same results had it been used.

(i) Isothermal spheres.

Isothermal spheres belong to standard literature in astronomy . This well posed mathematical problem can be solved numerically and provides a not too crude model for cores of globular clusters and played an important role in their understanding.

(a) Equilibrium configurations.

The problem at hand is the one treated in section $`2(\mathrm{𝑖𝑖𝑖})`$ with a spherical enclosure of radius $`r_B`$ say. Antonov  has given an elegant proof that in a sphere, spherical configurations maximize the entropy. Spherical equilibria have a gravitational potential $`U`$ which depends on the distance to the center $`r`$ only and Poisson’s equation (2.30) with a Boltzmann distribution (2.31) reduces in this case to

$$\frac{1}{r^2}\frac{d}{dr}\left[r^2\frac{dU(r)}{dr}\right]=4\pi G\rho (r)\rho (r)=m_{\mathrm{}}^+\mathrm{}fd^3p=m\left(\frac{2\pi m}{\beta }\right)^{3/2}e^{\alpha m\beta U}.$$
(5.1)
The boundary conditions for regular solutions are:

$$\left(\frac{dU}{dr}\right)_{r=0}=0\mathrm{and}U(r_B)U_B=\frac{GNm}{r_B}=\frac{GM}{r_B}.$$
(5.2)
Equation (5.1) has a singular solution

$$\rho _{sing}=\frac{1}{2\pi Gm\beta r^2}.$$
(5.3)
Regular solutions oscillate around this singular solution to which they become closer and closer as $`r`$ grows.

Excellent tables of numerical solutions have been given by Emden , more detailed ones are found in Chandrasekhar  with some additions in Chandrasekhar and Wares . The tables are still useful today even if Mathematica provides $`10^2`$ exact digits in a fraction of a second for solutions with density contrasts up to $`10^8`$ in no time.

The density contrast $`\rho (r=0)/\rho (r=r_B)`$ is a useful parameter which characterizes the relative strength of gravity. It is obtained by direct integration of equation (5.1) because $`U\mathrm{ln}\rho `$. Equilibrium configurations have been studied in detail in Chandrasekhar  but also in Lynden-Bell and Wood  and in Padmanabhan . A short review on isothermal spheres relevant to astronomy is given in Binney and Tremaine .

(b) Stable configurations.

($`\alpha `$) Isolated spheres and isothermal spheres in a heat bath.

Antonov  found analytically that stable isolated systems can exist for density contrasts111111In general we give numerical results with 3 digits. $`<709`$. A nice and simple analytical proof has been given by Padmanabhan . Lynden-Bell and Wood  found the same result using a combination of Poincaré’s turning point method and direct calculation of the entropy. They also found that equilibrium configurations maintained at constant temperature are stable for $`<32.1`$. It is worth noting that Padmanabhan’s method allows to calculate the density perturbations $`\delta \rho `$ that trigger instability. This gives interesting additional information that is beyond thermodynamic’s capability.

The method of linear series with conjugate thermodynamic functions give these results in one stroke and adds new interesting information. Figure 2 represents $`\beta m|U_B|`$ as a function of $`E/M|U_B|`$; $`N`$, $`r_B`$ and thus $`U_B`$ are constant.

For small $`\beta `$’s, i.e. high temperatures, gravity plays a minor role, the particles behave like a perfect gas and perfect gases are stable for arbitrary perturbations . The curve of figure 2 spirals inwards counterclockwise towards a point with coordinates $`E/M|U_B|=\frac{1}{4},\beta m|U_B|=2`$. If we follow the stable branch from low density contrasts, say from $`E=0`$ where $`6.85`$, to high density contrasts and negative energies, a vertical tangent appears at $`709`$. Thus isolated isothermal spheres become unstable for $`>709`$. This confirms once again Antonov’s result but figure 2 also shows that all configurations with a density contrast $`>709`$ are unstable with respect to spherically symmetric perturbations and that the number of unstable modes increases with $``$ because the linear series spirals inwards counterclockwise.

Notice as expected that stable configurations near the vertical tangent at $`709`$ have negative heat capacities. Notice however that unstable configurations with $`5.2210^3<<4.5010^4`$ have also $`C_V<0`$ though they are unstable and have one mode of instability: $`\lambda _1<0<\lambda _2`$.

Notice finally that all configurations are stable for non-spherically symmetric perturbations  but this has not and can not be proven by thermodynamic arguments since we have linear series of spherically symmetric equilibria only.

If we turn figure 2 clockwise $`90^{}`$ we are looking at $`E(\beta )`$. This is the appropriate linear series for canonical ensembles as we have seen in section 2. We can follow again the line of stable configurations from $`E=0`$ and $`=6.85`$ towards higher density contrasts with negative energies, we shall meet a vertical tangent at $`=32.1`$. Thus isothermal spheres in a heat bath become unstable for $`>32.1`$. This confirms Lynden-Bell and Wood’s result and shows in addition that all configurations with higher density contrasts are not only unstable but more and more so as $`E(\beta )`$ spirals inwards counterclockwise.

Notice as expected that stable configurations near the vertical tangent have $`C_V>0`$ but unstable configurations with two modes of instability $`\lambda _1<\lambda _2<0<\lambda _3`$ and density contrasts $`4.5010^4<<5.4510^5`$ have also $`C_V>0`$.

Chavanis used a method similar to that of Padmanabhan to calculate the density perturbations $`\delta \rho `$ that trigger instability in isothermal spheres at constant temperature  as well as in the ensembles studied in the next subsection .

($`\beta `$) Grand canonical ensembles and grand microcanonical ensembles.

The method of linear series with conjugate thermodynamic functions has also been applied to systems in which neither energy nor the number of particles are fixed. In grand canonical ensembles $`\alpha ,\beta ,r_B`$ are fixed. An appropriate linear series for this ensemble is, see section 3, the curve $`E(\beta )`$ at fixed $`\alpha `$ and $`r_B`$. This linear series, see also , is shown in figure 3A and 3B and must be looked at rotated $`90^{}`$ clockwise. What is represented is actually $`\beta ^{}=\beta /\beta _0`$ as a function of $`E^{}=(\beta _0^2Gm^2/r_B)E`$ with $`\beta _0=2\pi ^3m^7(8\pi Gr_B^2e^\alpha )^2`$. For $`E>0`$ , the non-dimensional quantities are of comparable magnitude. Notice that near $`E=0`$ there is a cross over but not a bifurcation with two different solutions whose stability or instability is known. For $`E<0`$ the curve winds in counterclockwise staying very near the vertical axis. The limit point has coordinates $`(E^{}=5.3710^3,\beta ^{}=13.6)`$. This is why the diagram has a different scale for $`E<0`$, figure 3a, and $`E>0`$, figure 3b.

At low density contrasts and high temperature the energy tends to zero; kinetic and potential energies are both high and of comparable magnitude. It has been shown  by calculating fluctuations of $`w_{GC}`$ near the extremum, see (2.40), that grand canonical ensembles are stable with respect to arbitrary perturbations121212Paper  contains correct calculations but the wrong interpretation. The authors believed they were dealing with isolated systems when they dealt in fact with a grand canonical ensemble. for $`<1.58`$. We can see on figure 3 that grand canonical ensembles have indeed a turning point at $`=1.58`$, $`E(\beta )`$ has a vertical tangent there and turns counterclockwise. At the next vertical tangent $`106`$ and the curve turns again counterclockwise and at higher density contrasts the curve is spiralling inwards. There is thus no stable equilibrium configuration for $`>1.58`$.

The situation is more eventful in the grand microcanonical ensemble in which $`E,\alpha ,r_B`$ are fixed. Here figure 3a and 3b are appropriate linear series for checking stability. At low density contrasts equilibrium configurations are certainly stable: the ensembles are more constrained than the grand canonical one which is stable. We see that stable configurations exist for density contrasts $`<1.66`$. Configurations with $`1.66<<11.6`$ have one mode of instability. There exists a second series of equilibrium configurations with $`11.6<<92.6`$. At higher density contrasts the curve spirals inwards counterclockwise and equilibrium configurations become more and more unstable.

(c) Fluctuations.

Fluctuations displace the limits of stability towards lower density contrasts. For instance in isolated isothermal spheres it was found  that quadratic fluctuations would induce instability at density contrasts $`_l709e^{3.30N^{1/3}}`$. If $`N=10^3`$, a number close to numerical experiments, $`_l510`$. The formula is still marginally valid for $`N=10`$ for which $`_l154`$. These are not small effects.

In this connection, it is worthwhile recalling Monaghan’s  application of the theory of hydrodynamic fluctuations to a self-gravitating gas. He showed that density fluctuations become large before the point of ordinary stability is reached. Thermodynamics supports Monaghan’s finding.

(d) Relevance to astronomy.

The simple model of isothermal spheres played an important role in understanding the structure of globular clusters. The model taken seriously by observers  is that of Michie  which was put to extensive use by King . The Michie-King model is a truncated Boltzmann distribution which in our notation is this, see (2.31):

$$f=A(e^{\beta E_0}e^{\beta E_c})\mathrm{for}E_0E_c\mathrm{and}f=0\mathrm{for}E_0>E_c.$$
(5.4)
$`E_c`$ is an energy cutoff that simulates the absence of high energy escapers. Stability limits for this model and variants thereof  do not change the general trend of equilibrium configurations . Linear series are always counterclockwise inwinding spirals and stable configurations exist up to some maximum density contrast that varies not very much.

The distribution (5.4) is a simple model for classifying observable parameters of globular clusters.

It is not clear what ensemble represents best thermal equilibrium in cores of globular clusters. However, that a stability limit exists is now commonly accepted. Globular clusters in quasi-thermal equilibrium evolve slowly towards greater entropy and density contrast due to stellar evaporation . They then reach a limit beyond which they cannot stay isothermal. The great contribution of Antonov  was to point out this instability.

What happens when a core becomes unstable was explained by Lynden-Bell . The explanation is now part of standard texts . Stars diffuse towards the center which becomes denser and denser while isothermal equilibrium is lost. This phenomena known as the gravothermal catastrophe leads eventually to core collapse. This has been confirmed by numerous studies mentioned in Meylan and Heggie .

Observational evidence by various authors, see in particular Trager, Dorjovsky and King’s paper , came in the late 1980’s when CCD observations allowed a systematic investigation of the inner surface brightness profile of globular clusters. These authors classified globular clusters into two different classes. About $`80\%`$ with a projected density profile that fitted Michie-King models and the rest with a density profile corresponding to a singular density in $`1/r^2`$. These cores are considered to have collapsed.

The work of Antonov and Lynden-Bell is the single most important contribution to observational astronomy based on thermodynamics of self-gravitating systems.

(ii) Maclaurin and Jacobi ellipsoids.

Liquid ellipsoids like isothermal spheres belong to standard literature. Chandrasekhar  devoted a whole book to the subject. Liquid ellipsoids have few independent variables and plenty of bifurcations. We give here one example in which a bifurcation is lifted with a small perturbation of the potential energy as discussed in Section 3. Solutions of Poisson’s equation of equilibrium which involve elliptic integrals are taken from Chandrasekhar’s book. The stability analysis presented here is taken from .

(a) Equilibrium configurations.

Consider a self-gravitating ellipsoid with uniform density $`\rho `$ and semi-axis $`abc`$. It rotates around the $`c`$-axis with uniform and constant angular velocity $`\mathrm{{\rm Y}}`$. The total mass $`M`$ and angular momentum $`L`$ are constant:

$$M=\frac{4\pi }{3}\rho abc,L=\frac{1}{5}M(a^2+b^2)\mathrm{{\rm Y}}I\mathrm{{\rm Y}}.$$
(5.5)
The motion described in comoving coordinates along $`a,b,c`$ has an effective potential energy $`V`$ which is, see for instance  :

$$V=\frac{1}{2}\rho Ud^3x+\frac{L^2}{2I}.$$
(5.6)
$`U`$ is the gravitational potential. There are three variables in $`V`$, $`a,b,c`$, and one constraint $`M`$. There are thus two independent variables, say,

$$x^1(\frac{a}{c})^2x^2(\frac{b}{c})^21.$$
(5.7)
The explicit form of $`V`$ is given in Chandrasekhar; in our notations,

$`V`$ $`=`$ $`{\displaystyle \frac{3GM^2}{10}}\left({\displaystyle \frac{4\pi \rho }{3M}}\right)^{1/3}F(x^1,x^2;s)\mathrm{with}s={\displaystyle \frac{25L^2}{3GM^3}}\left({\displaystyle \frac{4\pi \rho }{3M}}\right)^{1/3}\mathrm{and}`$
$`F`$ $`=`$ $`(x^1x^2)^{1/6}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\nu }{\sqrt{(1+\nu )(x^1+\nu )(x^2+\nu )}}}s{\displaystyle \frac{(x^1x^2)^{\frac{1}{3}}}{x^1+x^2}}.`$ (5.8)
Equilibrium configurations $`X^1(s)`$ and $`X^2(s)`$ are solutions of

$$\frac{V}{x^1}=\frac{V}{x^2}=0.$$
(5.9)
There exists one class of solutions to these equations for $`s<0.769`$: the Maclaurin spheroids with $`a=b`$ and an eccentricity $`e=(1\frac{c^2}{a^2})^{\frac{1}{2}}<0.813`$. For $`s0`$, $`\frac{a}{c}=\frac{b}{c}1`$; spheroids turn into spheres of infinite radius and zero mass. When $`s>0.769`$ there are two classes of solutions: Maclaurin spheroids with $`e>0.813`$ and Jacobi ellipsoids with non equal axis. At $`s=0.769`$ there is thus a bifurcation.

The solutions are shown and discussed in details in Chandrasekhar’s book together with their stability. Here we use the topology of linear series with conjugate variables with respect to minus the potential energy to find stability conditions.

(b) Stability conditions.

Let $`F_e(s)`$ represent the extremal values of $`F`$. $`F_e(s)`$ is equal to minus the potential energy times a constant. Equilibrium configurations are thus stable if $`F_e(s)`$ is a maximum of $`F`$. The conjugate parameter of $`s`$ with respect to $`F_e(s)`$ is

$$K(s)\frac{dF_e(s)}{ds}=\frac{(X^1X^2)^{\frac{1}{3}}}{X^1+X^2}.$$
(5.10)
The linear series reproduced from  is shown in figure 4. We see the bifurcation appearing at point $`B`$ where $`s=0.769`$. Since homogeneous spheres are stable  Maclaurin spheroids are stable for $`0<s<0.769`$ and $`0<e<0.813`$. A change of stability can only occur at point $`B`$ as we know but we cannot tell what happens at a bifurcation point just by looking at the linear series.

Let us then look as suggested in (3.13) at solutions of modified equations like, for instance,

$$\frac{V}{x^1}=ϵ\mathrm{and}\frac{V}{x^2}=0,$$
(5.11)
with $`ϵ`$ small. The effect of $`ϵ`$ is to break the symmetry and lift the bifurcation. Figure 4 represents also the linear series $`K(s)`$ for the perturbed solutions with $`ϵ=0.01`$. Any other small value of $`ϵ`$ would have a similar effect on the topology but curves with $`ϵ>0`$ cut those with $`ϵ=0`$ and the drawing is not so nice. We now see a continuous line (sequence 1) which connects smoothly for $`ϵ0`$ to the linear series of stable Maclaurin configurations for $`s<0.769`$ and to Jacobi ellipsoids which exists only for $`s>0.769`$. Sequence 1 has no vertical tangent. Thus by continuity Jacobi ellipsoids must be stable. Sequence 2 has a vertical tangent at point $`C`$. The $`CD`$ branch identifies with the stable sequence 1 when $`ϵ0`$, and the $`CE`$ branch represents necessarily unstable configurations with one unstable mode131313The two Poincaré coefficients of stability are different at point $`B`$, see .. Taking $`ϵ0`$, we conclude that Maclaurin spheroids with eccentricities $`e>0.813`$ are unstable.

The stability limits have been known since Poincaré’s time . The instability of Maclaurin spheroids for $`e>0.813`$ is secular; it shows up with a small viscous dissipation in a dynamical perturbation calculation .

(iii) Final comments and more examples.

(a) Final Comments

Thermodynamics of self-gravitational systems is helpful for calculating statistical equilibrium configurations, thermodynamic functions, stability limits and the effect of fluctuations near instability.

Globular clusters and other astronomical objects are only close to statistical equilibrium and evolve more often slowly. When evolution is slow enough these objects pass through a series of quasi-equilibrium configurations with ever increasing entropy. In this respect, linear series also mimic evolutionary tracks. Turning points indicate the limits beyond which evolution proceeds in a non isothermal way, generally towards higher central densities and core collapse. The fluctuations move somewhat the stability limits to lower temperatures.

What happens next, once thermal equilibrium no more exists, is of immense interest but does not belong to this paper. Thermodynamics is a small chapter in theoretical astronomy though an interesting one.

(b) More applications

Thermodynamic methods as described in the present work have been used in various studies which we shall now mention briefly. One dimensional isothermal parallel sheets have been shown to be stable to one dimensional perturbations . Isothermal axially-symmetric equilibria of gravitating rods are generally stable with respect to cylindrically symmetric perturbations . There is no gravothermal catastrophe as in isothermal spheres but no equilibria exist below some finite temperature. In contact with a heat bath slightly colder than that, the system collapses slowly, giving up an unlimited amount of energy. The effect of a mass spectrum on the stability limits of isothermal spheres has also been analyzed .

Isothermal spheres with a Fermi-Dirac distribution behave in some ways like systems of hard spheres particles: they have phase transitions from gaseous to core-halo structure. Stability limits with a Fermi-Dirac distribution were analyzed by Chavanis and Sommeria , see also . A cutoff in potential has a similar effect . Another example in which there may be phase transitions from white dwarfs to neutron stars but where the turning point method has not been used, though it would have proven extremely useful, is the study by Harrison, Thorne Wakano and Wheeler  of cold catalyzed matter. They minimize the mass-energy $`M`$. The only independent parameter is the number of baryons $`N`$. A plot of the conjugate parameter $`(M)/N`$ versus $`N`$ using their own table on page 152-153 would have shown the stability conditions and in particular the number of unstable modes in one stroke  as is clear from figure 5.

Chavanis  also studied isothermal spheres in Newtonian mechanics but with a relativistic equation of state. He found that instability sets in at smaller density contrasts but stronger binding energies.

Stability conditions in relativistic spheres with energy cutoffs  have also not been treated by the turning point method but it is obvious from figures in  that the their results might have been obtained in a more simple way.

Since Bekenstein  attributed an entropy to black holes and Hawking  found that they emit black body radiation of quantum origin, studies of their stability in various surroundings and with various charges have flourishes. After the seminal papers by Hawking  and by Gibbons and Perry on the thermodynamics of black holes, there appeared a series of papers based on linear series devoted to Kerr black holes  and Kerr-Newman black holes . A more detailed analysis of black holes in a cavity which takes account of the effect of fluctuations started in  was analyzed in greater details in . All those papers neglect the mass-energy of the cavity, usually a thin shell, that is supposed to play the role of a massless reflective wall. It has been pointed out  that thin shells which do not allow infinite tensions141414Thin shells that satisfy the dominant energy condition $`\rho <P<\rho `$ with $`\rho >0`$. must contain at least $`30\%`$ of the total mass-energy. This may affect considerably the stability limits and throws some doubts on considerations, with sometimes far reaching consequences, that were derived from black holes in equilibrium with radiation in weightless cavities .

Acknowledgements

This work grew out of a series of lectures at the Laboratoire de Physique Théorique at the University of Orsay in November 2001. The lectures entitled “Thermodynamique des gaz autogravitants” covered a great deal more material than the present much centered work. I wish to thank Jean-Philippe Uzan who gave me the opportunity to rethink this interesting subject. Useful conversations ensued with many people. I thank in particular Nathalie Deruelle and Bernard Jancovici whose clarifying remarks were particularly helpful.

Donald Lynden-Bell read a first draft of this paper and made highly valuable comments which were incorporated in the text. Pierre-Henri Chavanis who read the second draft contributed additional useful remarks which I also incorporated.

I am also grateful to Evgeni Sorkin, Ehud Nakar and Jonathan Oren for their help with the cosmetics of the figures.

I thank the Laboratoire de Physique Théorique at Orsay for a warm welcome.

Appendices

## Appendix A Proof that (2.23) can be transformed into (2.28)

We start from (2.23) with $`H=H_K+H_P`$ defined in (2.21). The trick consists in introducing first a “continuous density” $`\rho (𝐫)`$:

$$\rho (𝐫)=\underset{i}{}m\delta (𝐫𝐫_i).$$
(A.1)
The potential energy $`H_P`$ can then be written

$$H_P=\frac{1}{2}_V\frac{G\rho (𝐫)\rho (𝐫^{})}{|𝐫𝐫^{}|}d^3rd^3r^{}.$$
(A.2)
Next we discretize the whole space, cutting it into small cubes with an indice $`a`$ and volume $`\kappa ^3`$. Replacing $`\rho (𝐫)d^3r`$ in (A.2) in each cube with indices $`a`$ by $`\rho ^a\kappa ^3`$ the potential energy takes the form

$$H_P=\frac{1}{2}\underset{a,b}{}u_{ab}\rho ^a\rho ^b.$$
(A.3)
The matrix $`u_{ab}`$ is the discrete inverse of the Laplacian and more precisely

$$\underset{a,b\mathrm{continuum}}{lim}u_{ab}^1=\frac{\mathrm{\Delta }}{4\pi G}\delta (𝐫𝐫^{}).$$
(A.4)
Let $`W_a`$ represent the mean field in discretized form. Define a set of variables $`X_a`$ to be integrated out soon:

$$X_a=W_a\beta ^{1/2}u_{ab}\rho ^b.$$
(A.5)
The following identity151515This old identity due to Stratonovich  has proved useful in many-body physics. is readily constructed with (A.3) and (A.5) by elementary algebraic manipulations:

$$\beta H_P+\frac{1}{2}u^{1ab}X_aX_b=\beta ^{1/2}W_a\rho ^a+\frac{1}{2}u^{1ab}W_aW_b.$$
(A.6)
Then taking the exponent of both sides of this expression and integrating over the whole domain of variation of the $`W_a`$’s and $`X_a`$’s leads to:

$$e^{\beta H_P}=\frac{1}{[\mathrm{det}(2\pi u_{ab})]^{\frac{1}{2}}}_{\mathrm{}}^+\mathrm{}e^{\left(\beta ^{\frac{1}{2}}W_a\rho ^a+\frac{1}{2}u^{1ab}W_aW_b\right)}\underset{a}{}dW_a.$$
(A.7)
We now return to continuity, taking into account formula (A.2) and (A.4); this lets us write $`e^{\beta H_P}`$ in the following form:

$$e^{\beta H_P}=\frac{1}{}_{\mathrm{}}^+\mathrm{}e^{\left[m\underset{i}{}\beta ^{\frac{1}{2}}W(𝐫_i)+\frac{1}{8\pi G}_{\mathrm{}}^+\mathrm{}W(𝐫)\mathrm{\Delta }W(𝐫)d^3r\right]}𝒟W,$$
(A.8)
where

$$=\underset{a,b\mathrm{continuum}}{lim}\sqrt{det(2\pi u_{ab})}\text{and}𝒟W=\underset{a\mathrm{continuum}}{lim}\underset{a}{}dW_a.$$
(A.9)
$``$ is certainly divergent; $`u_{ab}`$ tends to $`1/|𝐫𝐫^{}|`$. We assume that a short distance cutoff will make $``$ convergent. The nature of the cutoff is unimportant in our approximation - see section 1(iii)(c). The expression for $`e^{\beta H_P}`$ obtained in (A.8) we put back into $`\mathrm{\Omega }(E)`$ in (2.23). It takes now little work to see that $`\mathrm{\Omega }`$ can be written as follows

$$\mathrm{\Omega }=\frac{1}{(2\pi i)}_{bi\mathrm{}}^{b+i\mathrm{}}_{\mathrm{}}^+\mathrm{}e^{\left[\beta E+\frac{1}{8\pi G}_{\mathrm{}}^+\mathrm{}W(𝐫)\mathrm{\Delta }W(𝐫)d^3r\right]}\mathrm{\Psi }𝑑\beta 𝒟W,$$
(A.10)
where

$$\mathrm{\Psi }=\frac{1}{N!}\left(e^{\beta [\frac{1}{2m}p^2+m\beta ^{\frac{1}{2}}W(𝐫)]}𝑑\omega \right)^N\mathrm{is}\mathrm{of}\mathrm{the}\mathrm{form}\mathrm{\Psi }=\frac{1}{N!}\psi ^N.$$
(A.11)
This is the of the same form as $`\stackrel{~}{\mathrm{\Psi }}`$ in (2.6). We may thus replace $`\mathrm{\Psi }`$ with the same approximate expression as the one obtained for $`\stackrel{~}{\mathrm{\Psi }}`$ in (2.6) to (2.9) :

$$\mathrm{\Psi }=\frac{1}{N!}\psi ^N\frac{1}{2\pi i}_{ai\mathrm{}}^{a+i\mathrm{}}\left[\underset{N^{}=0}{\overset{\mathrm{}}{}}\frac{1}{N^{}!}\psi ^N{}_{}{}^{}e_{}^{\alpha (N^{}N)}\right]d\alpha =\frac{1}{2\pi i}_{ai\mathrm{}}^{a+i\mathrm{}}[\underset{N^{}=0}{\overset{\mathrm{}}{}}\frac{1}{N^{}!}f^N{}_{}{}^{}]e^{\alpha N}d\alpha ,$$
(A.12)
in which

$$f=e^{\alpha \beta \left[\frac{1}{2m}p^2+m\beta ^{\frac{1}{2}}W(𝐫)\right]}.$$
(A.13)
We can also write $`\mathrm{\Psi }`$ like this

$$\mathrm{\Psi }=\frac{1}{2\pi i}_{ai\mathrm{}}^{a+i\mathrm{}}e^{\left(\alpha N+{\scriptscriptstyle f𝑑\omega }\right)}𝑑\alpha ;$$
(A.14)
replacing $`\mathrm{\Psi }`$ by this expression in (A.10) gives $`\mathrm{\Omega }(E)`$ as written in equation (2.28) with $`\sigma `$ shown in equation (2.29).

We followed closely the calculations161616Paper  contains correct calculations but the wrong interpretation. The authors believed they were dealing with isolated systems when they dealt in fact with a grand canonical ensemble. in  and  . We left aside all sorts of intricacies about existence and convergence problems in functional integrations. These were given more carefully considerations in  with appropriate references.

## Appendix B Calculation of $`\mathrm{\Omega }(E)`$ to order two in $`\delta W`$ \- Formula (2.35)

We start from $`\sigma (\alpha ,\beta ,W)`$ given in equation (2.29) which we expand in a Taylor series to order two in $`\delta \alpha =\alpha \alpha _e,\delta \beta =\beta \beta _e,\delta W=WW_e`$ around their extremal values $`\alpha _e,\beta _e,W_e`$ defined by (2.30)-(2.32). There is no reason à priori to limit $`\delta W`$ to continuous perturbations except for the fact  that non-continuous functions form a subset of measure zero. Here are the first order derivatives:

$`{\displaystyle \frac{\sigma }{\alpha }}`$ $`=`$ $`N+{\displaystyle f𝑑\omega },{\displaystyle \frac{\sigma }{\beta }}=E{\displaystyle \left[\frac{1}{2m}p^2+\frac{1}{2}m\beta ^{\frac{1}{2}}W(𝐫)\right]f𝑑\omega },`$
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{\sigma }{W}}\delta Wd^3r`$ $`=`$ $`m\beta ^{\frac{1}{2}}{\displaystyle \delta Wf𝑑\omega }+{\displaystyle \frac{1}{4\pi G}}{\displaystyle _{\mathrm{}}^+\mathrm{}}\delta W\mathrm{\Delta }Wd^3r.`$ (B.1)
It will be of great help to define with an over-bar the mean values in the phase space of one particle. For instance:

$$\overline{X}=\frac{1}{N}Xf_e𝑑\omega $$
(B.2)
Thus, the mean value of the total energy, equation (2.32), can be written as follows:

$$E=E_K+E_P=N\overline{\left(\frac{p^2}{2m}\right)}+\frac{1}{2}mN\overline{U}=N\left(\frac{3}{2\beta _e}+\frac{1}{2}m\overline{U}\right)\frac{\beta _eE}{N}=\frac{3}{2}+\frac{1}{2}m\beta _e\overline{U}.$$
(B.3)
This expression will soon be useful. Another useful relation is the virial equality:

$$E+E_K=3P_bV\stackrel{(\text{3.2})}{}\frac{\beta _eE}{N}+\frac{3}{2}=3\frac{\rho _b}{(mN/V)}=3\frac{\rho _b}{\rho _{mean}}.$$
(B.4)
Here now are the second derivatives of $`\sigma `$ with respect to $`\alpha ,\beta ,W`$ calculated at the point of extremum:

$`\left({\displaystyle \frac{^2\sigma }{^2\alpha }}\right)_e=N`$ , $`\left({\displaystyle \frac{^2\sigma }{\alpha \beta }}\right)_e=E,`$
$`\left({\displaystyle \frac{^2\sigma }{^2\beta }}\right)_e`$ $`=`$ $`{\displaystyle \frac{N}{4}}\left[{\displaystyle \frac{15}{\beta _e^2}}+{\displaystyle \frac{7m}{\beta _e}}\overline{U}+m^2(\overline{U})^2+m^2\overline{\left(U\overline{U}\right)^2}\right],`$
$`{\displaystyle \left(\frac{^2\sigma }{\alpha W}\right)_e\delta Wd^3r}`$ $`=`$ $`Nm\beta _e^{\frac{1}{2}}\overline{\delta W},`$
$`{\displaystyle \left(\frac{^2\sigma }{\beta W}\right)_e\delta Wd^3r}`$ $`=`$ $`N\left[m\beta _e^{\frac{1}{2}}\overline{\delta W}+{\displaystyle \frac{1}{2}}m^2\beta _e^{\frac{1}{2}}\overline{(U\delta W)}\right],`$
$`{\displaystyle \delta W\left(\frac{^2\sigma }{WW^{}}\right)_e\delta W^{}d^3rd^3r^{}}`$ $`=`$ $`Nm^2\beta _e\overline{(\delta W)^2}+{\displaystyle \frac{1}{4\pi G}}{\displaystyle \delta W\mathrm{\Delta }\delta Wd^3r}.`$ (B.5)
With these derivatives we may write the expression for $`\delta ^2\sigma `$ at the extremum:

$`2\delta ^2\sigma `$ $`=`$ $`\underset{1}{\underset{}{N(\delta \alpha )^2}}+\underset{2}{\underset{}{(2E\delta \alpha \delta \beta )}}`$ (B.6)
$`+`$ $`\underset{3}{\underset{}{(N/4)\left[15\beta _e^2+7m\beta _e^1\overline{U}+m^2(\overline{U})^2+m^2\overline{(U\overline{U})^2}\right](\delta \beta )^2}}`$
$`+`$ $`\underset{4}{\underset{}{(2Nm\beta _e^{\frac{1}{2}}\delta \alpha \overline{\delta W})}}+\underset{5}{\underset{}{N\left[2m\beta _e^{\frac{1}{2}}\overline{\delta W}+m^2\beta _e^{\frac{1}{2}}\overline{(U\delta W)}\right]\delta \beta }}+\underset{6}{\underset{}{Nm^2\beta _e\overline{(\delta W)^2}}}`$
$`+`$ $`\underset{7}{\underset{}{(4\pi G)^1{\displaystyle _{\mathrm{}}^+\mathrm{}}\delta W\mathrm{\Delta }\delta Wd^3r}}.`$
We now use (B.3) and (B.4) to simplify some of the terms. Our aim is to reduce $`\delta ^2\sigma `$ to a sum of squares that can be separately integrated in $`\alpha ,\beta `$ and $`W`$ or equivalently in $`\delta \alpha ,\delta \beta `$ and $`\delta W`$. We assume as usual that the exponent is steep and that order two is a good approximation. $`\delta \alpha `$ and $`\delta \beta `$ vary along a line parallel to the imaginary axis between $`\pm \mathrm{}`$.

First we write

$$\underset{1+2+3}{\underset{}{\mathrm{Terms}}}=\underset{1^{}}{\underset{}{N\left(\delta \alpha \frac{\beta E}{N}\frac{\delta \beta }{\beta _e}\right)^2}}+\underset{2^{}}{\underset{}{Nb^2\left(\frac{\delta \beta }{\beta _e}\right)^2}}b^2=\frac{1}{4}\left[\frac{6\rho _b}{\rho _{mean}}+m^2\beta _e^2\overline{(U\overline{U})^2}\right]>0.$$
(B.7)
Second,

$`\underset{4+5}{\underset{}{\mathrm{Terms}}}`$ $`=`$ $`\underset{3^{}}{\underset{}{2Nm\beta _e^{\frac{1}{2}}\left(\delta \alpha {\displaystyle \frac{\beta _eE}{N}}{\displaystyle \frac{\delta \beta }{\beta _e}}\right)\overline{\delta W}}}+\underset{4^{}}{\underset{}{2N\overline{(g\delta W)}{\displaystyle \frac{\delta \beta }{\beta _e}}}},`$
$`\mathrm{where}\overline{(g\delta W)}`$ $`=`$ $`m\beta _e^{\frac{1}{2}}\left(1{\displaystyle \frac{\beta _eE}{N}}\right)\overline{\delta W}+{\displaystyle \frac{1}{2}}m^2\beta _e^{\frac{3}{2}}\overline{(U\delta W)}.`$ (B.8)
With these two results we can now see that the quadratic sum $`\delta ^2\sigma `$ reduces to this:

$`2\delta ^2\sigma `$ $`=`$ $`\underset{1^{}}{\underset{}{N\left(\delta \alpha {\displaystyle \frac{\beta _eE}{N}}{\displaystyle \frac{\delta \beta }{\beta _e}}\right)^2}}+\underset{2^{}}{\underset{}{Nb^2\left({\displaystyle \frac{\delta \beta }{\beta _e}}\right)^2}}+\underset{3^{}}{\underset{}{\left[2Nm\beta _e^{\frac{1}{2}}\left(\delta \alpha {\displaystyle \frac{\beta _eE}{N}}{\displaystyle \frac{\delta \beta }{\beta _e}}\right)\overline{\delta W}\right]}}`$ (B.9)
$`+`$ $`\underset{4^{}}{\underset{}{2N\overline{(g\delta W)}{\displaystyle \frac{\delta \beta }{\beta _e}}}}+\underset{6}{\underset{}{Nm^2\beta _e\overline{(\delta W)^2}}}+\underset{7}{\underset{}{(4\pi G)^1{\displaystyle \delta W\mathrm{\Delta }\delta Wd^3r}}}.`$
The sum of $`1^{}`$ and $`3^{}`$, it is a difference of positive quantities:

$`\underset{1^{}+3^{}}{\underset{}{\mathrm{Terms}}}`$ $`=`$ $`N\left(\delta \alpha {\displaystyle \frac{\beta _eE}{N}}{\displaystyle \frac{\delta \beta }{\beta _e}}\right)^22Nm\beta ^{\frac{1}{2}}\left(\delta \alpha {\displaystyle \frac{\beta _eE}{N}}{\displaystyle \frac{\delta \beta }{\beta _e}}\right)\overline{\delta W}`$ (B.10)
$`=`$ $`N\left(\delta \alpha {\displaystyle \frac{\beta _eE}{N}}{\displaystyle \frac{\delta \beta }{\beta _e}}m\beta _e^{\frac{1}{2}}\overline{\delta W}\right)^2Nm^2\beta _e(\overline{\delta W})^2,`$
next we sum $`2^{}`$ and $`4^{}`$ which is of the same type as $`1^{}+3^{}`$:

$$\underset{2^{}+4^{}}{\underset{}{\mathrm{Terms}}}=Nb^2\left(\frac{\delta \beta }{\beta _e}\right)^2+2N\overline{(g\delta W)}\frac{\delta \beta }{\beta _e}=N\left[b\frac{\delta \beta }{\beta _e}+\frac{\overline{(g\delta W)}}{b}\right]^2N\left[\frac{\overline{(g\delta W)}}{b}\right]^2.$$
(B.11)
We now observe that $`2\delta ^2\sigma =\mathrm{eq}.(\text{B.10})+\mathrm{eq}.(\text{B.11})+\mathrm{Terms}(6+7)`$. Let us introduce a couple of real variables which vary between $`\pm \mathrm{}`$, similar to those introduced in equations (2.17):

$$i\alpha ^{}=\delta \alpha \frac{\beta _e}{N}\frac{\delta \beta }{\beta _e}m\beta _e^{\frac{1}{2}}\overline{\delta W}\mathrm{and}i\beta ^{}=b\frac{\delta \beta }{\beta _e}+\frac{\overline{(g\delta W)}}{b}.$$
(B.12)
In terms of $`\alpha ^{}`$, $`\beta ^{}`$ and $`\delta W`$, $`\delta ^2\sigma `$ reduces to

$`\delta ^2\sigma `$ $`=`$ $`{\displaystyle \frac{1}{2}}N(\alpha ^2+\beta ^2)+\delta ^2\sigma ^{}\mathrm{in}\mathrm{which}`$
$`\delta ^2\sigma ^{}`$ $`=`$ $`{\displaystyle \frac{1}{8\pi G}}{\displaystyle \delta W\mathrm{\Delta }\delta Wd^3r}+{\displaystyle \frac{N}{2}}\left[m^2\beta _e\overline{(\delta W\overline{\delta W})^2}b^2\{\overline{(g\delta W)}\}^2\right].`$ (B.13)
With this result we can calculate $`w`$ to order two starting from (2.28), (2.29); inserting $`\sigma =S+\delta ^2\sigma `$ into (2.28) we obtain

$$\mathrm{\Omega }(E)\frac{1}{(2\pi i)^2}_{\mathrm{}}^+\mathrm{}𝒟W_{\mathrm{}}^+\mathrm{}_{\mathrm{}}^+\mathrm{}e^{S+\delta ^2\sigma }𝑑\alpha 𝑑\beta =\frac{\beta _ee^S}{2\pi bN}_{\mathrm{}}^+\mathrm{}𝒟We^{\delta ^2\sigma ^{}}.$$
(B.14)
This is equivalent to formula (2.35).

It is worth noticing that $`\delta ^2\sigma `$, see (B.9), is equal to terms $`(6+7)`$ minus two squares, see (B.10) and (B.11). This implies that:

$$\delta ^2\sigma ^{}\mathrm{Terms}(6+7)=\frac{1}{2}_{\mathrm{}}^+\mathrm{}\delta W\left(\frac{\mathrm{\Delta }}{4\pi G}+m\beta _e\rho \right)\delta Wd^3r=\delta ^2\sigma _{GC}$$
(B.15)
$`\delta ^2\sigma _{GC}`$ is the second order term of the exponent for the grand canonical ensemble $`\sigma _{GC}`$ introduced in equation (2.39). The spectrum of eigenvalues of this quadratic form has been studied in detail in  for isothermal spheres where it was shown that in a stable system $`\delta ^2\sigma _{GC}\mathrm{}`$. Thus $`\delta ^2\sigma \mathrm{}`$ as well for such configurations. See section 1(iii)(c) about this divergence.

Figure captions

Figure 1: This figure illustrates the point made in section 4 subsection (i) where full details are given.

Figure 2: The figure represents the pair of conjugate thermodynamic functions $`\beta (E)`$ in appropriate units for isothermal spheres described in section 4(i)$`(\alpha )`$. The linear series gives stability limits in either isolated spheres or spheres in a heat bath. The parameter along the line, at vertical and at horizontal tangents is the density contrast $`=\rho (center)/\rho (boundary)`$.

Figure 3a and 3b: These figures represent the same pair of conjugate thermodynamic functions $`\beta (E)`$ as in figure 2 with different units (defined in the text) and is appropriate to detect stability limits of isothermal spheres in a grand canonical or a grand micro-canonical ensemble as described in section 4(i)$`(\beta )`$. Figure 3a is for $`0.008<E<0`$. Figure 3b is for $`0<E<0.6`$. For density contrasts $`1`$, $`E0`$ and $`\beta \mathrm{}`$. Thus the line with density contrasts $`>1.66`$ will, before $`E`$ turns negative, cuts the line that comes from infinity .

Figure 4: This figure represents the pair of conjugate variables for Maclaurin and Jacobi ellipsoids described in 4(ii). The figure is reproduced from  in which our small $`s`$ was represented by $`S`$. $`S`$ here represents the entropy.

Figure 5: This is a topologically correct plot of $`(E)/N`$ as a function of the number of baryons in units of $`N_{}=10^{57}`$. The parameter along the line is the central density in $`grcm^3`$. The linear series shows six consecutive turning points calculated by Harrison. There are two stable branches corresponding respectively to cold white dwarfs and to neutron stars. More details are given in .