A&A 432, 117-138 (2005)
DOI: 10.1051/0004-6361:20041114

On the lifetime of metastable states in self-gravitating systems

P. H. Chavanis

Laboratoire de Physique Théorique, Université Paul Sabatier, 118 route de Narbonne, 31062 Toulouse, France

Received 17 April 2004 / Accepted 3 November 2004

Abstract
We discuss the physical basis of the statistical mechanics of self-gravitating systems. We show the correspondance between statistical mechanics methods based on the evaluation of the density of states and partition function and thermodynamical methods based on the optimization of a thermodynamical potential (entropy or free energy). We address the question of the thermodynamic limit of self-gravitating systems, the justification of the mean-field approximation, the validity of the saddle point approximation near the transition point, the lifetime of metastable states and the fluctuations in isothermal spheres. In particular, we emphasize the tremendously long lifetime of metastable states of self-gravitating systems which increases exponentially with the number of particles N except in the vicinity of the critical point. More specifically, using an adaptation of the Kramers formula justified by a kinetic theory, we show that the lifetime of a metastable state scales as ${\rm e}^{N\Delta s}$ in microcanonical ensemble and  ${\rm e}^{N\Delta j}$ in canonical ensemble, where $\Delta s$ and $\Delta j$ are the barriers of entropy and free energy $j=s-\beta \epsilon$ per particle respectively. The physical caloric curve must take these metastable states (local entropy maxima) into account. As a result, it becomes multi-valued and leads to microcanonical phase transitions and "dinosaur's necks'' (Chavanis 2002b, [arXiv:astroph/0205426]; Chavanis & Rieutord 2003, A&A, 412, 1). The consideration of metastable states answers the critics raised by D.H.E. Gross [cond-mat/0307535/0403582].

Key words: gravitation - stellar dynamics - globular clusters: general - methods: miscellaneous

   
1 Introduction

The statistical mechanics of self-gravitating systems has a long history starting with the seminal papers of Antonov (1962) and Lynden-Bell & Wood (1968). A statistical mechanics approach is particularly relevant to describe the late stages of "small'' groups of stars ( $N\sim 10^{6}$), such as globular clusters, which evolve under the influence of stellar encounters ("collisional'' relaxation). Apart from astrophysical applications, the statistical mechanics of stellar systems is of great interest in physics because it differs in many respects from that of more familiar systems with short-range interactions (Padmanabhan 1990; Chavanis 2002c). In particular, for systems with long-range interactions, the thermodynamical ensembles are not equivalent, negative specific heats are allowed in the microcanonical ensemble (but not in the canonical ensemble) and metastable equilibrium states can have tremendously long lifetimes making them of considerable interest.

Two types of approaches have been developed to determine the statistical equilibrium state of a self-gravitating system. In the thermodynamical approach, one determines the most probable distribution of particles by maximizing the Boltzmann entropy at fixed mass and energy in the microcanonical ensemble or by minimizing the free energy F=E-TS at fixed mass and temperature in the canonical ensemble (Lynden-Bell & Wood 1968; Katz 1978; Chavanis 2002a). This approach is the simplest and the most illuminating. In addition, it is directly related to kinetic theories (based on the Landau or on the Fokker-Planck equation) for which the Boltzmann entropy (or the Boltzmann free energy) plays the role of a Lyapunov functional and satisfies a H-theorem. Alternatively, in the statistical mechanics approach, one starts from the density of states or partition function, transforms it into a functional integral and uses a saddle point approximation valid in a properly defined thermodynamic limit (Horwitz & Katz 1978; de Vega & Sanchez 2002; Katz 2003).

In the first part of this paper, we discuss the connexion between these two procedures. We remain at a heuristic level, stressing more the physical ideas than the mathematical formalism. In Sect. 2, we introduce the entropy by a combinatorial analysis. In order to regularize the problem at short distances, we consider either the case of self-gravitating fermions or the case of self-gravitating particles with a soften potential. We also discuss the thermodynamic limit of the classical and quantum self-gravitating gas. In Sect. 3, we show the relation between the density of states g(E) and the entropy functional S[f] and between the partition function $Z(\beta)$ and the free energy functional $J[f]=S[f]-\beta E[f]$. In the thermodynamic limit, the saddle point approximation amounts to maximizing the entropy at fixed mass and energy (microcanonical ensemble) or to minimizing the free energy at fixed mass and temperature (canonical ensemble). In Sect. 4, we discuss the notion of canonical and microcanonical phase transitions in self-gravitating systems. We perform the (standard) horizontal and (less standard) vertical Maxwell constructions and discuss the validity of the saddle point approximation near the transition point for finite N systems. These results (e.g., microcanonical first order phase transitions) are relatively new in statistical mechanics and still subject to controversy (Gross 2003, 2004). Therefore, we provide a relatively detailed discussion of these issues.

In the second part of the paper, we emphasize the importance of metastable states in astrophysics and show how they can be taken into account in the statistical approach. In Sect. 5, we use the Kramers formula to estimate the lifetime of a metastable state. We show that the lifetime of a metastable state scales as ${\rm e}^{N\Delta s}$ in microcanonical ensemble and  ${\rm e}^{N\Delta j}$ in canonical ensemble, where $\Delta s$ and $\Delta j$ are the barriers of entropy and free energy $j=s-\beta \epsilon$ per particle respectively. Therefore, the typical lifetime of a metastable state scales as  ${\rm e}^{N}$ except in the vicinity of the critical point $E_{\rm c}$ (Antonov energy) or $T_{\rm c}$ (Emden-Jeans temperature). We explicitly compute the barriers of entropy and free energy close to the critical point for classical self-gravitating particles (stars). The very long lifetime of metastable states, scaling as ${\rm e}^{N}$, was pointed out by Chavanis & Rieutord (2003) and the difficulty of a stellar system to overcome the entropic barrier and collapse was qualitatively discussed in Chavanis & Sommeria (1998). We here improve these arguments by developing a theory of fluctuations in isothermal spheres, following the approach of Katz & Okamoto (2000). We also determine how finite N effects affect the collapse temperature and the collapse energy. Finally, in Sect. 6, we derive a Fokker-Planck equation for the evolution of the distribution of energies P(E,t) in the canonical ensemble and make contact with the standard Kramers problem. We determine the typical lifetime of a metastable state by calculating the escape time accross a barrier of free energy.

   
2 The most probable distribution

   
2.1 The Fermi-Dirac distribution

We consider a system of N particles confined within a spherical box of radius R and interacting via Newtonian gravity. Let $f({\vec r},{\vec v},t)$ denote the distribution function of the system, i.e. $f({\vec r},{\vec v},t){\rm d}^{3}{\vec r} {\rm d}^{3}{\vec v}$ gives the mass of particles whose position and velocity are in the cell $({\vec r},{\vec v};{\vec r}+{\rm d}^{3}{\vec r},{\vec v}+{\rm d}^{3}{\vec v})$ at time t. The integral of f over the velocity determines the spatial density

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

and the total mass of the configuration is given by

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

where the integral extends over the entire domain. On the other hand, in the meanfield approximation, the total energy of the system can be expressed as

$$% E={1\over 2}\int fv^{2}{\rm d}^{3}{\vec r}{\rm d}^{3}{\vec v}+{1\over 2}\int\rho\Phi {\rm d}^{3}{\vec r}=K+W,$$ (3)

where K is the kinetic energy and W the potential energy. The meanfield expression of the potential energy is obtained from the exact expression

 

$$\biggl \langle -{1\over 2}\sum_{i\neq j}{Gm^{2}\over \vert{\vec r}_{i}-{\vec r}_{j}\vert}\biggr \rangle$$                                     W =

$$-{1\over 2}GN(N-1)m^{2}\int {P_{2}({\vec r}_{1},{\vec r}_{2})\ove... ...t{\vec r}_{1}-{\vec r}_{2}\vert}{\rm d}^{3}{\vec r}_{1}{\rm d}^{3}{\vec r}_{2},$$ = (4)

by approximating the two-body distribution function $P_{2}({\vec r}_{1},{\vec r}_{2})$ by the product of two one-body distribution functions $P_{1}({\vec r}_{1})$ $\times$ $P_{1}({\vec r}_{2})$ and using $\rho({\vec r})=Nm P_{1}$. For self-gravitating systems, this mean-field approximation is exact in a proper thermodynamic limit $N\rightarrow +\infty$ with $\eta=\beta GMm/R$ and $\Lambda=-ER/GM^{2}$ fixed (see Appendix A). The gravitational potential $\Phi=-G\int \rho({\vec r}')/\vert{\vec r}-{\vec r}'\vert{\rm d}^{3}{\vec r}'$ is solution of the Newton-Poisson equation

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

In order to regularize the problem at short distances, we shall invoke quantum mechanics and use the Pauli exclusion principle. The Pauli exclusion principle is a fundamental concept in physics and it has also applications in astrophysics, e.g. in white dwarf and neutron stars. Therefore, it can be considered as a physically relevant small-scale regularization for compact objects (e.g., Chavanis 2002d). We wish to determine the most probable distribution of self-gravitating fermions at statistical equilibrium. To that purpose, we divide the individual phase space $\lbrace {\vec r},{\vec v}\rbrace$into a very large number of microcells with size (h/m)³ where his the Planck constant (the mass m of the particles arises because we use ${\vec v}$ instead of ${\vec p}$ as a phase space coordinate). A microcell is occupied either by 0 or 1 fermion (or g=2s+1fermions if we account for the spin). We shall now group these microcells into macrocells each of which contains many microcells but remains nevertheless small compared to the phase-space extension of the whole system. We call $\nu$ the number of microcells in a macrocell. Consider the configuration $\lbrace n_i \rbrace$ where there are n₁ fermions in the 1st macrocell, n₂ in the 2nd macrocell etc., each occupying one of the $\nu$ microcells with no cohabitation. The number of ways of assigning a microcell to the first element of a macrocell is $\nu$, to the second $\nu -1$ etc. Since the particles are indistinguishable, the number of ways of assigning microcells to all n_i particles in a macrocell is thus

$$% {1\over n_i!}{\times} {\nu!\over (\nu-n_i)!}\cdot$$ (6)

To obtain the number of microstates corresponding to the macrostate  $\lbrace n_i \rbrace$ defined by the number of fermions n_i in each macrocell (irrespective of their precise position in the cell), we need to take the product of terms such as Eq. (6) over all macrocells. Thus, the number of microstates corresponding to the macrostate  $\lbrace n_i \rbrace$, i.e. the probability of the state  $\lbrace n_i \rbrace$, is

$$% W(\lbrace n_i \rbrace)=\prod_i {\nu!\over n_i!(\nu-n_i)!}\cdot$$ (7)

This is the Fermi-Dirac statistics. As is customary, we define the entropy of the state  $\lbrace n_i \rbrace$ by

$$% S(\lbrace n_i \rbrace)=\ln W(\lbrace n_i \rbrace).$$ (8)

It is convenient here to return to a representation in terms of the distribution function giving the phase-space density in the ith macrocell

$$% f_i=f({\vec r}_i,{\vec v}_i)={n_i \ m\over \nu \ ({h\over m})^3}={n_i\eta_0\over \nu},$$ (9)

where we have defined $\eta_0=m^{4}/h^3$, which represents the maximum value of f due to Pauli's exclusion principle. Now, using the Stirling formula, we have

 

$$% \ln W(\lbrace n_i \rbrace)\simeq \sum_i \nu (\ln\nu-1)-\nu\bigg... ...iggl (1-{f_i\over \eta_0}\biggr )\biggr\rbrace-1\biggr\rbrack\biggr\rbrace\cdot$$ (10)

Passing to the continuum limit $\nu\rightarrow 0$, we obtain the usual expression of the Fermi-Dirac entropy

$$% S=-k_B\int \biggl\lbrace {f\over\eta_{0}}\ln {f\over\eta_{0... ...rm d}^{3}{\vec r}{\rm d}^{3}{\vec v}\over ({h\over m})^3}\cdot$$ (11)

If we take into account the spin of the particles, the above expression remains valid but the maximum value of the distribution function is now $\eta_{0}=g m^{4}/h^{3}$, where g=2s+1 is the spin multiplicity of the quantum states (the phase space element has also to be multiplied by g). In the non-degenerate (or classical) limit $f\ll\eta_0$, the Fermi-Dirac entropy (11) reduces to the Boltzmann entropy

$$% S=-k_{\rm B}\int {f\over m}\biggl\lbrack \ln\biggl ({f h^3\... ...iggr )-1\biggr \rbrack {\rm d}^{3}{\vec r}{\rm d}^{3}{\vec v}.$$ (12)

Now that the entropy has been precisely justified, the statistical equilibrium state (most probable state) of self-gravitating fermions is obtained by maximizing the Fermi-Dirac entropy (11) at fixed mass (2) and energy (3):

$$% {\rm Max}\quad S[f]\quad \vert \quad E[f]=E,\ M[f]=M.$$ (13)

Introducing Lagrange multipliers 1/T (inverse temperature) and $\mu$ (chemical potential) to satisfy these constraints, and writing the variational principle in the form

$$% \delta S-{1\over T}\ \delta E+{\mu\over T} \delta N=0,$$ (14)

we find that the critical points of entropy correspond to the Fermi-Dirac distribution

$$% f={\eta_{0}\over 1+\lambda {\rm e}~^{\beta m ({v^{2}\over 2}+\Phi)}},$$ (15)

where $\lambda={\rm e}^{-\beta \mu}$ is a strictly positive constant (inverse fugacity) and $\beta={1\over k_{\rm B} T}$ is the inverse temperature. Clearly, the distribution function satisfies $f\le \eta_{0}$, which is a consequence of Pauli's exclusion principle.

So far, we have assumed that the system is isolated so that the energy is conserved. If now the system is in contact with a thermal bath (e.g., a radiation background) fixing the temperature, the statistical equilibrium state minimizes the free energy F=E-TS, or maximizes the Massieu function $J=S-\beta E$, at fixed mass and temperature:

$$% {\rm Max}\quad J[f]\quad \vert\quad M[f]=M.$$ (16)

Introducing Lagrange multipliers and writing the variational principle in the form

$$% \delta J+{\mu\over T} \delta N=0,$$ (17)

we find that the critical points of free energy are again given by the Fermi-Dirac distribution (15). Therefore, the critical points (first variations) of the variational problems (13) and (16) are the same. However, the stability of the system (regarding the second variations) can be different in microcanonical and canonical ensembles (see, e.g., Chavanis 2002b). When this happens, we speak of a situation of ensemble inequivalence. The stability of the system can be determined by a graphical construction, by simply plotting the series of equilibria $\beta(E)$ and using the turning point method of Katz (1978, 2003). Inequivalence of statistical ensembles occurs when the series of equilibria presents turning points or bifurcations.

   
2.2 Classical particles with soften gravitational potential

We now consider a system of classical self-gravitating particles, like stars in globular clusters. In order to make the problem of statistical mechanics well-posed mathematically (see below), we introduce a soften potential of the form

$$% u({\vec r}-{\vec r}')={-Gm^{2}\over \sqrt{({\vec r}-{\vec r}')^{2}+r_{0}^{2}}},$$ (18)

where r₀ is the soften radius. As we shall see, the soften radius r₀ plays a role similar to the inverse of $\eta_{0}$, the maximum phase space density, in the case of self-gravitating fermions. As said previously, this soften radius is introduced in order to pose the problem correctly. However, we shall argue in the sequel that this small-scale cut-off is irrelevant for the structure of stellar systems.

We wish to determine the most probable distribution of stars at statistical equilibrium (Ogorodnikov 1965). To that purpose, we divide the individual phase space $\lbrace {\vec r},{\vec v}\rbrace$ into a very large number of microcells with size (h/m)³ where h is a constant with dimension of angular momentum. Of course, quantum mechanics is not relevant for stellar systems so that h should not be confused with the Planck constant in the present context. For classical systems, a microcell can be occupied by an arbitrary number of particles. Adapting the counting analysis of Sect. 2.1 to the present context, the number of microstates corresponding to the macrostate  $\lbrace n_i \rbrace$, i.e. its probability, is

$$% W(\lbrace n_i \rbrace)=N!\prod_i {\nu^{n_{i}}\over n_i!}\cdot$$ (19)

This is the Maxwell-Boltzmann statistics. If we define the entropy of the state  $\lbrace n_i \rbrace$ by

$$% S(\lbrace n_i \rbrace)=\ln W(\lbrace n_i \rbrace),$$ (20)

and take the continuum limit, we obtain the usual expression of the Boltzmann entropy

$$% S=-k_{B} \int {f\over m}\ln\biggl ({f h^3\over N m^{4}}\biggr ) {\rm d}^{3}{\vec r}{\rm d}^{3}{\vec v}.$$ (21)

Note that it differs from the expression (12) obtained from the Fermi-Dirac entropy. This is of course related to the Gibbs paradox in standard thermodynamics (Huang 1963). In the absence of self-gravity, Eq. (21) reduces to the awkward expression

$$% S=Nk_{\rm B}\ln\biggl \lbrack V\biggl ({4\pi m\over 3h^{2}}{E\over N}\biggr )^{3/2}\biggr\rbrack+{3\over 2}N,$$ (22)

which is clearly non-extensive. By constrast, Eq. (12) leads to the Sackur-Tetrode formula

$$% S=Nk_{\rm B} \ln\biggl \lbrack {V\over N}\biggl ({4\pi m\over 3h^{2}}{E\over N}\biggr )^{3/2}\biggr\rbrack+{5\over 2}N,$$ (23)

which is extensive. As is well-known, the origin of this discrepency is due to the indiscernability of the particles and to the presence of the factor N! in the Maxwell statistics (19). For a molecular gas, the Gibbs paradox is usualy solved by invoking quantum mechanics. For a system of stars, one cannot use this argument. We shall consider that the stars are discernable and use the expression (21) for the entropy. However, this choice does not affect the structure of the equilibrium state as we shall see in the sequel.

The most probable distribution of stars at statistical equilibrium is now obtained by maximizing the Boltzmann entropy (21) at fixed mass and energy. This yields the Maxwell-Boltzmann distribution

$$% f=A {\rm e}~^{-\beta m ({v^{2}\over 2}+\Phi)},$$ (24)

where $\Phi$ is related to the density $\rho$ by

$$% \Phi=-G\int {\rho({\vec r}')\over \sqrt{({\vec r}-{\vec r}')^{2}+\epsilon_{0}^{2}}}{\rm d}^{3}{\vec r}'.$$ (25)

The microcanonical ensemble is the correct description of stellar systems which form an isolated Hamiltonian system in a first approximation. We can also consider the case of self-gravitating systems in contact with a thermal bath of non-gravitational origin which imposes its temperature T. For such systems, the correct description is the canonical ensemble and the statistical equilibrium state is obtained by minimizing the Boltzmann free energy F=E-TS at fixed mass. The canonical ensemble is also the correct description of a gas of self-gravitating Brownian particles (Chavanis et al. 2002). In this model, the friction and the stochastic fluctuations can mimick the influence of an external medium (thermostat) to which the system of origin is coupled.

   
2.3 Thermodynamic limit of self-gravitating systems

We introduce dimensionless variables such that ${\vec r}=R{\vec r}'$, ${\vec v}=U{\vec v}'$ and f=(M/R³U³)f' where R is the box radius, M is the mass of the system and $U\equiv (GM/R)^{1/2}$ is a typical velocity obtained by a Virial type argument (or dimensional analysis). For self-gravitating fermions, the entropy (11) can be expressed as

$$% S=-Nk_{\rm B}\mu\int {\rm d}^{3}{\vec r}'{\rm d}^{3}{\vec v... ...r \mu}\biggr )\ln\biggl (1-{f'\over \mu}\biggr )\biggr\rbrace,$$ (26)

where $\mu=\eta_0\sqrt{G^3 M R^3}$ is the degeneracy parameter (Chavanis & Sommeria 1998). Writing $\mu=(R/R_{*})^{3/2}$ with R_*=h²/GM^1/3m^8/3, we note that the degeneracy parameter is the ratio, to the power 3/2, of the system's radius divided by the radius R_* of a "white dwarf star'' (i.e. a completely degenerate ball of fermions) with mass M. The classical limit corresponds to $R\gg R_{*}$, i.e. $\mu\rightarrow \infty$. The conservation of mass is equivalent to

$$% \int {f'} {\rm d}^{3}{\vec r}'{\rm d}^{3}{\vec v}'=1,$$ (27)

and the conservation of energy is equivalent to

$$% {ER\over GM^{2}}=\int f' {v^{'2}\over 2}{\rm d}^{3}{\vec r}... ...r'}_{2}\vert}{\rm d}^{3}{\vec r}'_{1}{\rm d}^{3}{\vec r}'_{2}.$$ (28)

Finally, the Massieu function can be written

$$% J=N(s[f']+\eta \Lambda[f']),$$ (29)

where s=S/N, $\eta=\beta GMm/R$ and $\Lambda=-ER/GM^{2}$. We define the thermodynamic limit as $N\rightarrow +\infty$ such that $\mu=\eta_0\sqrt{G^3 M R^3}$, $\Lambda=-ER/GM^{2}$ and $\eta=\beta GMm/R$ are fixed. Coming back to physical quantities, it makes sense to fix h, m and G. Then, we have the scalings $R\sim N^{-1/3}$, $T\sim N^{4/3}$, $E\sim N^{7/3}$, $S\sim N$ and $J\sim N$ as $N\rightarrow +\infty$ (the free energy F scales as N^7/3). This is the quantum thermodynamic limit (QTL) for the self-gravitating gas (Chavanis 2002b; Chavanis & Rieutord 2003). This thermodynamic limit is relevant for compact objects with small radii $R\sim N^{-1/3}\ll 1$such as white dwarfs, neutron stars, fermion balls etc. The usual thermodynamic limit $N,R\rightarrow +\infty$ with N/R³ constant is clearly not relevant for inhomogeneous systems whose energy is non-additive (Padmanabhan 1990).

For classical particles with soften potential, the entropy (21) can be expressed as

$$% S=-N\int {f'}\biggl\lbrack \ln\biggl ({f'\over N\nu}\biggr )-1\biggr\rbrack {\rm d}^{3}{\vec r}'{\rm d}^{3}{\vec v}',$$ (30)

where $\nu\equiv {m^{4}\over h^{3}}\sqrt{G^{3}MR^{3}}$ is the counterpart of the degeneracy parameter. For classical particles, we see that it does not play any fundamental role in determining the structure of the system since it just appears as an additional constant term (independent on f) in the entropy. If we only consider the part of entropy that depends on the distribution function, we get

$$% S_{R}=-N\int {f'}\ln f' {\rm d}^{3}{\vec r}'{\rm d}^{3}{\vec v}'.$$ (31)

This is the relevant part of the entropy functional considered by Antonov (1962) and Lynden-Bell & Wood (1968). We can therefore write S[f]=S_R[f]+S_I where S_I is the constant part (irrelevant). The conservation of mass is equivalent to Eq. (27) and the conservation of energy is equivalent to

$$% {ER\over GM^2}=\int f' {v^{'2}\over 2}{\rm d}^{3}{\vec r}'{... ...epsilon^{2}}}{\rm d}^{3}{\vec r}'_{1}{\rm d}^{3}{\vec r}'_{2},$$ (32)

where $\epsilon=r_{0}/R$. As before, the Massieu function is given by Eq. (29). We define the thermodynamic limit as $N\rightarrow +\infty$ such that $\Lambda=-ER/GM^{2}$, $\eta=\beta GMm/R$ and $\epsilon=r_{0}/R$ are fixed. Coming back to physical quantities, it makes sense to fix r₀, m and G. Then, we have the scalings (Chavanis & Rieutord 2003) $R\sim 1$, $E\sim N^{2}$, $T\sim N$, $S_{R}\sim N$ and $J_{R}\sim N$ as $N\rightarrow +\infty$ (the free energy F_R scales as N²). These scalings imply that $\nu\sim N^{1/2}$ (if we fix h). Therefore, the (irrelevant) constant part of the entropy per particle diverges logarithmically as $S_{I}/N\sim \ln\nu\sim {1\over 2}\ln N\rightarrow +\infty$. This does not seem to be a crucial problem since this diverging term does not depend on f and therefore does not affects the structure of the equilibrium state. However, in a strict sense, there is no thermodynamic limit for classical self-gravitating particles with soften potential. This contrasts with the case of self-gravitating fermions that possess a rigorous thermodynamic limit (QTL).

Let us finally consider the case of classical self-gravitating particles without small-scale cut-off. The entropy is given by the Boltzmann formula (21). When r₀=0, we know that the Boltzmann entropy has no global maximum at fixed mass and energy (Antonov 1962). However, for sufficiently high energies, it has local entropy maxima that describe metastable gaseous states. The thermodynamic limit in that context corresponds to $N\rightarrow +\infty$ such that $\Lambda=-ER/GM^{2}$ and $\eta=\beta GMm/R$ are of order unity. If we fix m, G and T, we have the scalings $R\sim N$, $E\sim N$, $S\sim N$, $J\sim N$ and $F\sim N$ as $N\rightarrow +\infty$. This is the classical thermodynamic limit (CTL), or dilute limit, for the self-gravitating gas (de Vega & Sanchez 2002). Physically, it describes metastable gaseous states that are not affected by the small-scale cut-off (Chavanis & Rieutord 2003). As we shall see, these metastable states have considerably long lifetimes so that this thermodynamic limit is relevant for classical objects with large radii $R\sim N\gg 1$ such as globular clusters.

   
3 Connexion with statistical mechanics

   
3.1 Series of equilibria and metastable states

Figure 1: Series of equilibria for self-gravitating fermions with small cut-off/large $\mu$/large system size R. It has a Z-shape structure (dinosaur's neck). There can be several values of inverse temperature $\beta$ for a given energy E. They correspond to local maxima (LEM), global maxima (GEM) or saddle points (SP) of entropy S[f]. Similar remarks apply in the canonical ensemble where the role of (E, S) and ($\beta$, J) is reversed.

Open with DEXTER

Figure 2: Series of equilibria for self-gravitating fermions with large cut-off/small $\mu$/small system size R. It has a N-shape structure. There is only one value of inverse temperature $\beta$ for a given energy E in the microcanonical ensemble. It corresponds to a global maximum of entropy (GEM). By contrast, there are several values of energy E for a given $\beta$ in the canonical ensemble. They correspond to local maxima (LFEM), global maxima (GFEM) or saddle points (SP) of free energy J[f].

Open with DEXTER

	Figure 3: Series of equilibria for self-gravitating classical particles without small-scale cut-off. There exists only local entropy maxima (metastable state) or unstable saddle points of entropy in MCE. Similar remarks hold in CE for the free energy.
Open with DEXTER

The critical points of entropy S[f] at fixed E and M (i.e., the distribution functions  $f({\vec r},{\vec v})$ which cancel the first order variations of S at fixed E, M) form a series of equilibria parameterized, for example, by the density contrast ${\cal R}=\rho(0)/\rho(R)$ between the center and the edge of the system (see Chavanis 2002b). At each point in the series of equilibria corresponds a temperature $\beta$ and an energy E. In this approach, $\beta$ is the Lagrange multiplier associated with the conservation of energy in the variational problem (14). It has also the interpretation of a kinetic temperature in the Fermi-Dirac distribution (15). We can thus plot $\beta(E)$ along the series of equilibria. The form of this "caloric curve'' depends on the value of the degeneracy parameter $\mu$ in the case of fermions (Chavanis 2002b) and on the soften radius $\epsilon$ for regularized classical systems (Chavanis & Ispolatov 2002). It also depends on the dimension of space D (Sire & Chavanis 2002; Chavanis 2004a). In D=3, the caloric curve has a Z-shape (see Fig. 1) for small cut-off and a N-shape (see Fig. 2) for large cut-off (for no cut-off, we recover the well-known spiral of Fig. 3). There can be several values of temperature $\beta$ for the same energy E because the variational problem (13) can have several solutions: a local entropy maximum (metastable state), a global entropy maximum, and one or several saddle points. We must represent all these solutions on the caloric curve because local entropy maxima (metastable states) are in general more physical than global entropy maxima for the timescales achieved in astrophysics. Indeed, the system can remain frozen in a metastable gaseous phase for a very long time. This is the case, in particular, for globular clusters and for the gaseous phase of fermionic matter (at high energy and high temperature). The time required for a system placed in a metastable gaseous state to collapse is in general tremendously long and increases exponentially with the number N of particles. Thus, $t_{{\rm life}}\rightarrow +\infty$ in the thermodynamic limit $N\rightarrow +\infty$. The robustness of metastable states is due to the long-range nature of the gravitational potential. Therefore, at high temperatures and high energies, the global entropy maximum is not physically relevant. Condensed objects (e.g., planets, stars, white dwarfs, fermion balls,...) only form below a critical energy $E_{\rm c}$ (Antonov energy) or below a critical temperature $T_{\rm c}$ (Jeans temperature), when the gaseous metastable phase ceases to exist. The point where the metastable phase disappears is called a spinodal point.

   
3.2 Microcanonical ensemble

Let us explain things differently so as to make a closer contact with statistical mechanics. In statistical mechanics, one usually starts with the density of states

 

$$% g(E)=\int \delta \lbrack E-H({\vec r}_1,...,{\vec r}_N,{\vec v}... ...od_{i=1}^N {{\rm d}^{3}{\vec r}_i {\rm d}^3 {\vec v}_i\over ({h\over m})^{3N}},$$ (33)

where H is the Hamiltonian. For our system

$$% H={1\over 2}\sum_{i=1}^{N}mv_i^2-\sum_{i<j}{Gm^2\over \vert{\vec r}_i-{\vec r}_j\vert}\cdot$$ (34)

The density of states is the normalization factor of the N-body microcanonical distribution

$$% P_{N}({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N})={1\over g(E)}\delta (E-H),$$ (35)

stating that all accessible microstates are equiprobable.

Introducing the probability $W(\lbrace n_i\rbrace)$ of the state $\lbrace n_i \rbrace$, we can rewrite the density of states in the form

$$% g(E)=\sum_{E(\lbrace n_{i}\rbrace )=E}W(\lbrace n_i\rbrace),$$ (36)

where the sum runs over all macrostates with energy E. Introducing the entropy $S=\ln W$ of the state  $\lbrace n_i \rbrace$ and taking the continuum limit, the density of states can be expressed formally as

$$% g(E)=\int {\cal D}f {\rm e}^{S[f]}\delta (E-E[f])\delta (M-M[f]),$$ (37)

where the sum runs over all distribution functions and S[f] is the Fermi-Dirac entropy (11) if the particles are fermions and the Boltzmann entropy (21) for classical particles (in that case, the gravitational potential must be regularized otherwise the density of states (33) diverges). The probability of the distribution function  $f({\vec r},{\vec v})$ is proportional to  ${\rm e}^{S[f]}$. We now define the microcanonical entropy by $S_{{\rm micro}}(E)=\ln g(E)$ and the microcanonical temperature by $\beta_{{\rm micro}}=dS_{{\rm micro}}(E)/dE$. By definition, the caloric curve $\beta _{{\rm micro}}(E)$ is uni-valued (Gross 2003). In the thermodynamic limit defined in Sect. 2.3, the entropy S[f] scales as $\sim$N, that is S[f]=Ns[f] where $s\sim 1$ is the entropy per particle. Therefore, for $N\rightarrow +\infty$, the integral in Eq. (37) is dominated by the state  $f_{{\rm global}}({\vec r},{\vec v})$ which is the global maximum of S[f] at fixed M and E (rigorously speaking, we should work with the dimensionless quantities defined in Sect. 2.3 to get rid of the $N\rightarrow +\infty$ dependence). Then, $g(E)\simeq {\rm e}^{S[f_{{\rm global}}]}$, $S_{{\rm micro}}\simeq S[f_{{\rm global}}]$ and $\beta_{{\rm micro}}=\delta S/\delta E=\beta$. However, this approach fails to take into account metastable states (local maxima of S[f] at fixed M and E), which are of considerable interest in astrophysics. Indeed, equilibrium statistical mechanics tells nothing about timescales; a kinetic theory is required in that case. As explained above, these metastable states can persist for very long times. They correspond to the observed "diluted'' structures in the universe (e.g., globular clusters). Therefore, the caloric curve  $\beta _{{\rm micro}}(E)$ does not describe the system adequately. The series of equilibria $\beta(E)$ contain more information as they show local and global entropy maxima (as well as unstable saddle points). The curve  $\beta _{{\rm micro}}(E)$ can be deduced from $\beta(E)$ by keeping only global entropy maxima (see Fig. 4). If we adopt this procedure, we find that the system exhibits a first order microcanonical phase transition (provided that the system size $\mu$is sufficiently large) at a transition energy  $E_{\rm t}(\mu )$ where the gaseous phase and the condensed phase have the same entropy (see Chavanis 2002b). In the strict thermodynamic limit $N\rightarrow +\infty$, this phase transition is marked by a discontinuity of temperature. In fact, for finite N systems, the saddle point approximation $g(E)\simeq {\rm e}^{S[f_{{\rm global}}]}$ breaks down near the transition energy (Chavanis & Ispolatov 2002). This is because the contribution of the local entropy maximum in the functional integral (37) becomes more and more important as we approach $E_{\rm t}$. For the saddle point approximation to be valid, the number of particles must scale as $N\sim \vert\Lambda-\Lambda_{\rm t}\vert^{-1}$ for $\Lambda\rightarrow \Lambda_{\rm t}$ (see Sect. 4.2). Thus, for large but finite N, the temperature variation is sharp but remains continuous at the transition. We again emphasize that, due to the existence of metastable states, this phase transition may not be physically relevant. The true phase transition (gravothermal catastrophe) will rather take place at, or near, the spinodal point $E_{\rm c}$ (Antonov energy) where the metastable branch disappears. Estimating the lifetime of a metastable state by the Kramers formula  $t_{{\rm life}}\sim {\rm e}^{\Delta S}$, where $\Delta S$ is the height of the entropic barrier (difference of entropy between the local maximum and the saddle point), we find that  $t_{{\rm life}} \sim {\rm exp}\lbrack{2\lambda'N(\Lambda_{\rm c}-\Lambda)^{3/2}}\rbrack$ with $\lambda'\simeq 0.863159...$ (see Sect. 5). Except in the vicinity of the critical point  $\Lambda_{\rm c}$, the lifetime of a gaseous metastable state is enormous as it increases exponentially with the number of particles. Thus, metastable states have a considerable interest in astrophysics.

Figure 4: Strict microcanonical caloric curve for $\mu =10^{5}$ (the control parameter is the energy E). This figure is obtained from Fig. 1 by keeping only global entropy maxima. It corresponds therefore to the true caloric curve  $\beta _{{\rm micro}}(E)$ which is univalued (Gross 2003). For $N\rightarrow +\infty$, there is a discontinuity of temperature at the transition energy  $E_{\rm t}(\mu )$. For finite N systems, this discontinuity is smoothed out. Although this caloric curve is correct in a strict sense, it is not physical because it ignores metastable states that have an infinite lifetime in the thermodynamic limit. The physical caloric curve is obtained from Fig. 1 by discarding the unstable saddle points of entropy that form the intermediate branch (see Fig. 5).

Open with DEXTER

Figure 5: Physical microcanonical caloric curve for $\mu =10^{5}$ (the control parameter is the energy E). This figure displays stable and metastable equilibrium states. The metastable states have infinite lifetime in the thermodynamic limit so they are physically relevant. The first order microcanonical phase transition of Fig. 4 does not take place in practice. Due to the existence of metastable states, the system displays a microcanonical hysteretic cycle marked by a "collapse'' and an "explosion'' at the spinodal points where the branch of metastable states disapears (Chavanis & Rieutord 2003).

Open with DEXTER

If we now consider the case of classical particles ( $\hbar \rightarrow 0$ or $\mu\rightarrow +\infty$), the transition energy  $E_{\rm t}(\mu )$ is rejected to $+\infty$ so that the whole branch of gaseous states is metastable. This "no cut-off'' limit is relevant to classical objects such as globular clusters or to the interstellar medium because, for these systems, the size of the particles (stars and atoms) clearly does not matter. In that case, the series of equilibria $\beta(E)$forms a spiral (see Fig. 3) indicating the existence of one local entropy maximum and one (or several) saddle points of entropy for a given energy (Lynden-Bell & Wood 1968). This spiral is the limiting form, for $\mu\rightarrow +\infty$, of the fermionic caloric curve (see Fig. 11 in Chavanis 2002b). In this limit, the branch of "collapsed'' states (condensed phase) coincides with the x-axis where $\beta=0$. It corresponds to configurations made of two particles in contact ($\sim$binary star) surrounded by a hot halo with $T\rightarrow +\infty$. This "binary+halo'' configuration has an infinite entropy so, in a sense, it is the most probable configuration in the microcanonical ensemble (see Appendix A of Sire & Chavanis 2002). However, for sufficiently large energies (above the Antonov point), these configurations must be discarded "by hands'' because they are reached for inaccessibly large times. Therefore, for classical particles, the physical caloric curve  $\beta_{{\rm physical}}(E)$is obtained by taking $g(E)\simeq {\rm e}^{S[f_{{\rm local}}]}$ and $S_{{\rm physical}}\simeq S[f_{{\rm local}}]$ where  $f_{{\rm local}}$ is the local entropy maximum at fixed mass and energy (Chavanis 2003).

   
3.3 Canonical ensemble

In the canonical ensemble, the object of fundamental importance is the partition function

$$% Z(\beta)=\int {\rm e}^{-\beta H({\vec r}_1,...,{\vec r}_N,{... ...}^{D}{\vec r}_i {\rm d}^D {\vec v}_i\over ({h\over m})^{DN}},$$ (38)

which is the normalization of the N-body canonical distribution

$$% P_{N}({\vec r}_1,{\vec v}_1,...,{\vec r}_N,{\vec v}_N)={1\over Z(\beta)} {\rm e}^{-\beta H}.$$ (39)

The partition function can be rewritten

 

$$% Z(\beta) \sum_{\lbrace n_{i}\rbrace } {\rm e}^{-\beta E(\lbrace n_{i}\rbrace )} W(\lbrace n_{i}\rbrace )$$ =

$$\sum_{\lbrace n_{i}\rbrace } {\rm e}^{S(\lbrace n_{i}\rbrace )-\b... ...{i}\rbrace )} =\sum_{\lbrace n_{i}\rbrace } {\rm e}^{J(\lbrace n_{i}\rbrace )},$$ = (40)

where $J(\lbrace n_{i}\rbrace )=S(\lbrace n_{i}\rbrace )-\beta E(\lbrace n_{i}\rbrace )$ is the "free energy'' of the macrostate $\lbrace n_i \rbrace$ and the sum runs over all macrostates (in the present context, $J[f]=S[f]-\beta E[f]$ is a more natural functional than the usual free energy F[f]=E[f]-TS[f]). Taking the continuum limit, the partition function can be expressed formally as

$$% Z(\beta)=\int {\cal D}f {\rm e}^{J[f]}\delta (M-M[f]),$$ (41)

where the sum runs over all distribution functions and J[f] is the Fermi-Dirac free energy if the particles are fermions and the Boltzmann free energy for classical particles (in that case, the gravitational potential must be regularized otherwise the partition function (38) diverges). Note that Eq. (41) can also be obtained from the exact formula

$$% Z(\beta)=\int {\rm e}^{-\beta E}g(E){\rm d}E,$$ (42)

by substituting Eq. (37) for g(E) and carrying out the integration over E. We now define the canonical free energy by $F_{{\rm cano}}=-(1/\beta)\ln Z$. The average energy of the system at temperature T can be written $\langle E\rangle_{{\rm cano}}=-\partial\ln Z/\partial\beta$. By definition, the caloric curve $\langle E\rangle_{{\rm cano}}(\beta)$ is uni-valued. In the thermodynamic limit defined in Sect. 2.3, the free energy J[f] scales as $\sim$N. Therefore, for $N\rightarrow +\infty$, the integral in Eq. (41) is dominated by the state $f_{{\rm global}}({\vec r},{\vec v})$ which is the global maximum of J[f] at fixed M. Then, $Z(\beta)\simeq {\rm e}^{J[f_{{\rm global}}]}$, $F_{{\rm cano}}\simeq -(1/\beta)J[f_{{\rm global}}]=E[f_{{\rm global}}]-TS[f_{{\rm global}}]$ and $\langle E\rangle_{{\rm cano}}=-\delta J/\delta\beta=E$. Metastable states (local maxima of J[f] at fixed M) can be taken into account by plotting the full curve $E(\beta)$. It is obtained from $\beta(E)$ defined in Sect. 3.1 by simply reversing the graph since the critical points of the variational problems (13) and (16) are the same (see Sect. 2.1). The curve $\langle E\rangle_{{\rm cano}}(\beta)$ can be deduced from $E(\beta)$ by keeping only global maxima of free energy J (see Fig. 6). If we adopt this procedure, we find that the system exhibits a first order canonical phase transition at a transition temperature  $T_{\rm t}(\mu )$where the gaseous phase and the condensed phase have the same free energy (see Chavanis 2002b). This phase transition is marked by a discontinuity of energy (latent heat). In fact, the saddle point approximation $Z(\beta)\simeq {\rm e}^{J[f_{{\rm global}}]}$ breaks down near the transition temperature. For large but finite N, the energy variation is sharp but remains continuous at the transition. For the saddle point approximation to be valid, the number of particles must scale as $N\sim \vert\eta-\eta_{\rm t}\vert^{-1}$ for $\eta\rightarrow \eta_{\rm t}$(see Sect. 4.2). We again emphasize that, due to the existence of metastable states, this phase transition may not be relevant and that the physical phase transition (isothermal collapse) takes place at, or near, the spinodal point $T_{\rm c}$ (Jeans temperature) where the metastable branch disappears (see Fig. 8). Estimating the lifetime of a metastable state by the Kramers formula  $t_{{\rm life}}\sim {\rm e}^{\Delta F/k_{\rm B}T}$, where $\Delta F$ is the height of the potential barrier (difference of free energy between the local maximum and the saddle point), we find that $t_{{\rm life}} \sim {\rm exp}\lbrace 2\lambda N (\eta_{\rm c}-\eta)^{3/2}\rbrace$ with $\lambda\simeq 0.16979815...$ (see Sect. 5). Except in the vicinity of the critical point  $\eta _{\rm c}$, the lifetime of a gaseous metastable state is enormous as it increases exponentially with the number of particles. Metastable states are therefore highly robust. For classical objects ( $\hbar \rightarrow 0$), the transition temperature  $T_{\rm t}(\mu )$ is rejected to $+\infty$ so that the whole branch of gaseous states is metastable. In that case, the series of equilibria $E(\beta)$ forms a spiral (see Fig. 3) indicating the existence of one local minimum of free energy F and one (or several) saddle points of free energy for a given temperature. In the classical limit, the branch of "collapsed'' states (condensed phase) is rejected to $E\rightarrow -\infty$. It corresponds to configurations where all the particles have collapsed at r=0. This "Dirac peak'' configuration has an infinite free energy $F=-\infty$ (due to the infinite binding energy) so, in a sense, it is the most probable configuration in the canonical ensemble (see Appendix B of Sire & Chavanis 2002). This differs from the global entropy maximum made of a binary star surrounded by a hot halo in the microcanonical ensemble. For sufficiently large temperatures (above the Emden-Jeans point), these configurations must be discarded "by hands'' because they are reached for inaccessibly large times. Therefore, for classical particles, the physical caloric curve  $E_{{\rm physical}}(\beta)$is obtained by taking $Z(\beta)\simeq {\rm e}^{J[f_{{\rm local}}]}$ and $J_{{\rm physical}}\simeq J[f_{{\rm local}}]$ where  $f_{{\rm local}}$ is the local maximum of free energy J at fixed mass and temperature (Chavanis 2003).

Figure 6: Strict canonical caloric curve for $\mu =10^{3}$ (the control parameter is the inverse temperature $\beta$). This figure is obtained from Fig. 2 by keeping only global maxima of free energy J. It corresponds therefore to the true canonical caloric curve which is univalued and does not display negative specific heats contrary to the corresponding microcanonical caloric curve (see Fig. 7). For $N\rightarrow +\infty$, there is a discontinuity of energy (latent heat) at the transition temperature  $T_{\rm t}(\mu )$. For finite N systems, this discontinuity is smoothed-out. Although this canonical caloric curve is correct in a strict sense, it is not physical because it ignores metastable states that have an infinite lifetime in the thermodynamic limit $N\rightarrow +\infty$. The physical canonical caloric curve is obtained from Fig. 2 by discarding the unstable saddle points of free energy that form the intermediate branch with negative specific heats (see Fig. 8).

Open with DEXTER

Figure 7: Microcanonical caloric curve for the same cut-off/degeneracy parameter $\mu$/system's radius R as Fig. 6 (the control parameter is the energy E). It has a N-shape structure and displays a stable region with negative specific heats. All the equilibria are global maxima of entropy. In the canonical ensemble the region of negative specific heats is replaced by a phase transition (see Fig. 6). The curves of Figs. 6 and 7 are similar to those obtained by Padmanabhan (1990) with his toy model consisting in N=2 stars in gravitational interaction.

Open with DEXTER

   
3.4 Grand canonical ensemble

In the grand canonical ensemble, the partition function is

 

$$% Z_{{\rm GC}}(\beta,\mu)=\sum_{N=0}^{+\infty} {\rm e}^{N\mu\over k_{\rm B}T} Z_N(\beta).$$ (43)

Using Eq. (41), we get

 

$$% Z_{{\rm GC}} \sum_{N=0}^{+\infty}\int {\cal D}f {\rm e}^{J[f]}{\rm e}^{N\mu\over k_{\rm B}T}\delta (N-N[f])$$ =

$$\int {\cal D}f {\rm e}^{J[f]+\beta \mu N[f]}=\int {\cal D}f {\rm e}^{G[f]},$$ = (44)

where $G[f]=J[f]+\beta \mu N[f]$ is the grand potential. Of course, the expression (41) for Z_N is correct only for $N\gg 1$. However, the contribution of small N terms in the grand partition function (43) is expected to be weak so that Eq. (44) provides a good approximation of the series. Now, the grand potential G[f] scales as $\sim$ $N_{0}\equiv {R\over \beta Gm^{2}}$. Therefore, in the thermodynamic limit $R\rightarrow +\infty$with fixed $\beta$, G (gravitational constant) and m, the partition function  $Z_{{\rm GC}}$ is dominated by the distribution f which maximizes the grand potential G[f] at fixed $\beta$ and $\mu$. This problem has been considered for classical particles in D=3 (Chavanis 2003). We shall reserve for a future work the study of self-gravitating fermions in the grand canonical ensemble.

Figure 8: Physical canonical caloric curve for $\mu =10^{3}$ (the control parameter is the inverse temperature $\beta$). This figure displays stable and metastable equilibrium states. The metastable states have infinite lifetime in the thermodynamic limit so they are physically relevant. The first order canonical phase transition of Fig. 6 does not take place in practice. Due to the existence of metastable states the system displays a canonical hysteretic cycle marked by a "collapse'' and an "explosion'' at the spinodal points where the branch of metastable states disapears (Chavanis & Rieutord 2003).

Open with DEXTER

   
4 First order microcanonical and canonical phase transitions

   
4.1 Maxwell constructions and critical points

The deformation of the caloric curve when we vary the degeneracy parameter $\mu$/ system size R is represented in Figs. 9 and 10. Similar curves are obtained for a hard sphere gas or a soften potential (Chavanis & Ispolatov 2002) instead of fermions. In that case, $1/\mu$ plays the role of the cut-off radius a or soften radius r₀.

For $\mu<\mu_{{\rm CTP}}\simeq 82.5$, the curve $\beta(E)$ defining the series of equilibria is monotonic, so there is no phase transition. For $\mu>\mu_{{\rm CTP}}$, the curve $E(\beta)$ is multivalued so that a canonical first order phase transition is expected. The temperature of transition in the canonical ensemble can be obtained by a Maxwell construction as for the familiar Van der Waals gas. The equal area Maxwell condition A₁=A₂ (see Fig. 9) can be expressed as

$$% \int_{E_{A}}^{E_{\rm C}}(\beta-\beta_{t}){\rm d}E=0,$$ (45)

where E_A is the energy of the gaseous phase and $E_{\rm C}$ the energy of the condensed phase at the transition temperature $T_{\rm t}$. Since ${\rm d}S=\beta {\rm d}E$, one has

$$% S_{\rm C}-S_{A}-\beta_{\rm t}(E_{\rm C}-E_{A})=0.$$ (46)

Introducing the free energy $J=S-\beta E$, we verify that the Maxwell construction is equivalent to the equality of the free energy of the two phases at the transition:

$$% J_{A}=J_{\rm C}.$$ (47)

If we keep only global maxima of free energy (in the canonical ensemble) as in Fig. 6, the winding branch has to be replaced by a horizontal plateau. We see in Fig. 9 that the extent of the plateau decreases as $\mu$ decreases (dashed line). At the canonical critical point  $\mu_{{\rm CTP}}$, the plateau disappears and the curve presents an inflexion point.

For $\mu>\mu_{{\rm MTP}}\simeq 2600$, the curve $\beta(E)$ is multivalued so that a microcanonical first order phase transition is expected to occur (in addition to the canonical first order phase transition described previously). The energy of transition can be obtained by a vertical Maxwell construction. The equal area Maxwell condition A₁=A₂ (see Fig. 10) can be expressed as

$$% \int_{\beta_{A}}^{\beta_{\rm C}}(E-E_{\rm t}){\rm d}\beta=0,$$ (48)

where T_A and $T_{\rm C}$ are the temperatures of the two phases at the transition energy $E_{\rm t}$. Since ${\rm d}J=-E {\rm d}\beta$, one has

$$% J_{\rm C}-J_{A}+E_{\rm t}(\beta_{\rm C}-\beta_{A})=0.$$ (49)

Thus, the Maxwell construction is equivalent to the equality of the entropy of the two phases at the transition:

$$% S_{A}=S_{\rm C}.$$ (50)

If we keep only global maxima of entropy (in the microcanonical ensemble) as in Fig. 4, the winding branch has to be replaced by a vertical plateau. We see that the extent of the plateau decreases as $\mu$ decreases (dashed line). At the microcanonical critical point  $\mu_{{\rm MTP}}$, the plateau disappears and the curve presents an inflexion point.

	Figure 9: Horizontal Maxwell plateau associated with a canonical first order phase transition.
Open with DEXTER

	Figure 10: Vertical Maxwell plateau associated with a microcanonical first order phase transition.
Open with DEXTER

Therefore, for $\mu>\mu_{{\rm MTP}}$, we expect a microcanonical and a canonical first order phase transition, for $\mu_{{\rm CTP}}<\mu<\mu_{{\rm MTP}}$only a canonical first order phase transition and for $\mu<\mu_{{\rm CTP}}$no phase transition at all. We emphasize, however, that due to the presence of long-lived metastable states, the first order phase transitions and the plateaux are not relevant for the timescales of interest (see Sect. 3). Only the zeroth order phase transitions (gravothermal catastrophe and isothermal collapse) marked by a discontinuity of entropy or free energy are physically relevant.

Figure 11: Caloric curve of the self-gravitating Fermi gas with $\mu =1000$. In the canonical ensemble, the temperature of transition is determined by a Maxwell construction. For $\eta <\eta _{\rm t}$ ( $T>T_{\rm t}$), the gaseous states are global maxima of free energy J[f] and the condensed states local maxima (see the dashed curve). The situation is reversed for $\eta >\eta _{\rm t}$( $T<T_{\rm t}$).

Open with DEXTER

   
4.2 Validity of the saddle point approximation near the transition point

In this section, we discuss the validity of the saddle point approximation near the transition point. In the canonical ensemble, the partition function can be written

$$% Z(\beta)=\int {\rm e}^{-\beta E}g(E){\rm d}E,$$ (51)

where g(E) is the density of states with energy E. Introducing the entropy $S(E)=\ln g(E)$, we can rewrite the partition function as

 

$$% Z(\beta)=\int {\rm e}^{S(E)-\beta E}{\rm d}E=\int {\rm e}^{J(E)}{\rm d}E=\int {\rm e}^{Nj(E)}{\rm d}E,$$ (52)

where $J(E)=S(E)-\beta E$ is the "free energy''. As explained in Sect. 3, the equilibrium states correspond to maxima of J. The condition J'(E)=0, leading to $S'(E)=\beta$, determines a series of equilibria $E(\beta)$. We shall consider the case where the series of equilibria $E(\beta)$ has the N-shape structure of Fig. 11. The transition temperature  $T_{\rm t}(\mu )$ is determined by a Maxwell construction (see Sect. 4.1). From Eq. (52), the distribution of energies at temperature T is given by

$$% P(E)={1\over Z(\beta)}{\rm e}^{J(E)}.$$ (53)

It has a bimodal structure as shown in Fig. 12. The energies where P(E) is maximum (denoted E₁ and E₂) correspond to the stable (GFEM) and metastable (LFEM) states in Fig. 11. The energy where P(E) is minimum (denoted E_*) corresponds to the unstable states (SP) in Fig. 11. For  $T>T_{\rm t}$, gaseous configurations (high energies) are more probable than condensed configurations (low energies). This is the opposite for $T<T_{\rm t}$. Note that the notion of "more probable'' is delicate since the system can remain blocked in a metastable ("less probable'') state for very long time, making that state the physically most relevant state.

Figure 12: Distribution of energies $P(E)\propto {\rm exp}(S(E)-\beta E)$ in the canonical ensemble for different values of temperature (here $\mu =1000$). The entropy S(E) has been calculated in the mean-field approximation, using Eq. (66). The distribution of energies has a bimodal structure characteristic of a first order canonical phase transition. For $\eta =\eta _{\rm t}=1.052$, the two phases have the same probability. For $\eta =1<\eta _{\rm t}$, the gaseous phase (high energies) is the "most probable'' and for $\eta =1.1>\eta _{\rm t}$, the condensed phase (low energies) is the "most probable'' in a strict sense.

Open with DEXTER

For $N\rightarrow +\infty$, the partition function can be approximated by

$$% Z(\beta)={\rm e}^{Nj(E_{1})}+{\rm e}^{Nj(E_{2})},$$ (54)

where E₁ and E₂ are the energies at which J(E) is maximum. We now wish to obtain the strict caloric curve $\langle E\rangle_{{\rm cano}}(\beta)$ defined in Sect. 3. When T is not too close from the transition temperature $T_{\rm t}$ and $N\rightarrow +\infty$, we need just keep the contribution of the global maximum of free energy as explained previously. To investigate the situation close to the transition temperature, we rewrite the partition function (54) as

$$% Z(\beta)={\rm e}^{Nj(E_{1})}\biggl\lbrack 1+{\rm e}^{N(j(E_{2})-j(E_{1}))}\biggr\rbrack.$$ (55)

Now, close to the transition point, we have $j(E_{1})=j(E_{2})+\lambda^{2}(\eta-\eta_{t})$, where $\lambda$ is a constant of order unity depending on $\mu$ (see Fig. 11, dashed line). Therefore,

$$% Z(\beta)={\rm e}^{Nj(E_{1})}\biggl\lbrack 1+{\rm e}^{-N\lambda^{2}\Delta\eta}\biggr\rbrack,$$ (56)

where $\Delta\eta=\eta-\eta_{\rm t}$. Thus, the saddle point approximation is valid for $N\vert\eta-\eta_{\rm t}\vert\gg 1$. This requires increasing large values of N as we approach the transition temperature $T_{\rm t}$.

	Figure 13: Caloric curve of the self-gravitating Fermi gas with $\mu =10^{5}$. In the microcanonical ensemble, the energy of transition is determined by a Maxwell construction. For $\epsilon >\epsilon _{\rm t}$, the gaseous states are global maxima of entropy S[f] and the condensed states are local maxima (see the dashed curve). The situation is reversed for $\epsilon <\epsilon _{\rm t}$.
Open with DEXTER

In the microcanonical ensemble, the density of states can be written

$$% g(E)=\int {\cal D}f {\rm e}^{S[f]}\delta (E-E[f])\delta (M-M[f]).$$ (57)

We shall consider the case where the series of equilibria $\beta(E)$has the S-shape structure of Fig. 13. The transition energy  $E_{\rm t}(\mu )$ is determined by a Maxwell construction (see Sect. 4.1). For $N\rightarrow +\infty$, the density of states can be approximated by

$$% g(E)={\rm e}^{Ns[f_{1}]}+{\rm e}^{Ns[f_{2}]},$$ (58)

where f₁ and f₂ are the distribution functions corresponding to the global and local maxima of entropy (we shall denote by f_* the distribution function corresponding to the saddle point of entropy). The corresponding temperatures $\beta_{1}$, $\beta_{2}$ and $\beta_{*}$ form the series of equilibria $\beta(E)$ in Fig. 13. We now wish to obtain the strict caloric curve  $\beta _{{\rm micro}}(E)$ defined in Sect. 3. When E is not too close from the transition energy $E_{\rm t}$ and $N\rightarrow +\infty$, we need just keep the contribution of the global maximum of entropy. To investigate the situation close to the transition energy, we rewrite the density of states (58) as

$$% g(E)={\rm e}^{Ns[f_{1}]}\biggl\lbrack 1+{\rm e}^{N(s[f_{2}]-s[f_{1}])}\biggr\rbrack.$$ (59)

Now, close to the transition point, we have $s[f_1]=s[f_2]+\lambda'^{2}(\Lambda-\Lambda_{\rm t})$, where $\lambda'$ is a constant of order unity depending on $\mu$ (see Fig. 13). Therefore,

$$% g(E)={\rm e}^{Ns[f_{1}]}\biggl\lbrack 1+{\rm e}^{-N\lambda^2\Delta\Lambda}\biggr\rbrack,$$ (60)

where $\Delta\Lambda=\Lambda-\Lambda_{\rm t}$. Thus, the saddle point approximation is valid for $N\vert\Lambda-\Lambda_{\rm t}\vert\gg 1$. This requires increasing large values of N as we approach the transition point $E_{\rm t}$.

   
4.3 Comparison with other works

To conclude this section, we shall briefly compare our results with those of previous works. A complementary discussion is given in Chavanis (2002b).

The statistical mechanics of the classical self-gravitating gas ($\hbar=0$, no small scale cut-off) has been first considered by Antonov (1962), Lynden-Bell & Wood (1968) and Katz (1978) by using thermodynamical potentials. Lynden-Bell & Wood (1968) realized the importance of negative specific heats as a driving mechanism for the gravothermal catastrophe. The paradox posed by the negative heat capacity in statistical mechanics was solved by Thirring (1970) in terms of ensemble inequivalence and further discussed by Lynden-Bell & Lynden-Bell (1977) with the aid of a toy model. The occurence of negative heat capacities in the reverse gradient isotherms of all first order phase transitions was demonstrated by Lynden-Bell (1999). These results have been reviewed by Padmanabhan (1990) and complements have been given by de Vega & Sanchez (2002) and Chavanis (2002a, 2003).

The hard sphere problem was explored by Aronson & Hansen (1972) who gave an account of the phase transitions which had been suggested but not demonstrated for that system. A similar discussion has been given by Hertel & Thirring (1971) for self-gravitating fermions. These authors consider a large small-scale cut-off (or a small system size) and discuss first order phase transitions and Maxwell constructions in the canonical ensemble. This corresponds to Fig. 9. The distinction between the Emden-Jeans point $T_{\rm c}$ and the transition point  $T_{\rm t}(\mu )$ was emphasized by Kiessling (1989) and Stahl et al. (1994). They argued that the instability should take place at the transition point (where the free energies of the two phases coincide), i.e. at a higher temperature than the Emden-Jeans temperature: $T_{\rm t}>T_{\rm c}$. This is however not the case in practice because of the long lifetime of metastable states that we discuss here. The Emden-Jeans point $T_{\rm c}$, which really corresponds to Jeans instability threshold for box confined systems (see Chavanis 2002a) is the most physically relevant collapse point.

Stahl et al. (1994) considered the case of a low small-scale cut-off and discussed first order phase transitions in the microcanonical ensemble (see also Youngkins & Miller 2000). However, they did not give the caloric curves $\beta(E)$ and the vertical Maxwell constructions corresponding to Fig. 10. These curves are important to settle the thermodynamical stability of the system using the standard turning point argument (Katz 1978, 2003). Chavanis & Sommeria (1998) gave the series of equilibria of Fermi-Dirac spheres but did not discuss the nature of phase transitions in detail. A more complete discussion was given by Chavanis (2002b) for self-gravitating fermions and by Chavanis & Ispolatov (2002) for a soften potential. The phase diagrams in the $(\mu,T)$ and $(\mu,E)$planes were given by Chavanis & Rieutord (2003). These studies show how the series of equilibria depend on the "effective'' cut-off parameter $\mu$, leading to Z-shape and N-shape curves. Therefore, they make the bridge between the classical spiral of Lynden-Bell & Wood (no cut-off) and the N-shape caloric curve studied by Aronson & Hansen (large cut-off).

The effect of the dimension of space D on the gravitational phase transitions has been studied by Sire & Chavanis (2002) and Chavanis (2004a). Finally, gravitational phase transitions in rotating systems have been discussed by Votyakov et al. (2002) for a model on a lattice with a large spacing and by Chavanis & Rieutord (2003) for fermions. Dynamical models of gravitational phase transitions have been considered by Ispolatov & Karttunen (2003) in the microcanonical ensemble (by using molecular dynamics simulations) and by Chavanis et al. (2002) in the canonical ensemble (by introducing a model of self-gravitating Brownian particles). More recent simulations have been performed in the case of self-gravitating Brownian fermions (Chavanis et al. 2004).

   
5 The persistence of metastable states

   
5.1 Typical lifetime of a metastable state

In this section, we estimate the lifetime of a metastable state by using an adaptation of the Kramers formula (Risken 1989). We start first by the canonical ensemble. Close to the critical temperature $T_{\rm c}$, see Fig. 11, the free energy F[f] has one global minimum  $F_{{\rm global}}$ (stable GFEM), one local minimum  $F_{{\rm local}}$ (metastable LFEM) and one saddle point  $F_{{\rm saddle}}$(unstable SP). We call $E_{{\rm local}}$ the energy of the metastable equilibrium state and  $E_{{\rm saddle}}$ the energy of the unstable saddle point. For a system initially prepared at  $E_{{\rm local}}$, the probability of the energy E is $P(E)\sim {\rm e}^{-(F(E)-F_{{\rm local}})/k_{\rm B}T}$. Now, if the energy fluctuations drive the system past  $E_{{\rm saddle}}$, it will collapse. Therefore, the lifetime of the metastable state can be estimated by  $t_{{\rm life}}\sim 1/P(E_{{\rm saddle}})$ or

$$% t_{{\rm life}}\sim {\rm e}^{\Delta F/k_{\rm B}T},$$ (61)

where $\Delta F=\vert F_{{\rm saddle}}-F_{{\rm local}}\vert$ is the barrier of potential appropriate to our problem. Noting that

$$% t_{{\rm life}}\sim {\rm e}^{N \Delta j},$$ (62)

we conclude that, except in the vicinity of the critical point $T_{\rm c}$, the lifetime of a metastable state scales as ${\rm exp}(N)$. Therefore, metastable states (LFEM) are extremely robust in astrophysics and cannot be neglected, even if there exists states with lower free energy (GFEM). To investigate the situation close to the critical point $T_{\rm c}$, we shall calculate the barrier of potential $\Delta j$ for the classical self-gravitating gas ( $\hbar \rightarrow 0$). To that purpose, we use the results derived in a preceding paper (Chavanis 2002a).

We recall that the series of equilibria is parameterized by $\alpha=(4\pi G\beta\rho_{0})^{1/2}R$, where $\rho_{0}$ is the central density. Introducing the Milne variables

$$% u={\xi {\rm e}^{-\psi}\over \psi'}, \qquad v=\xi\psi',$$ (63)

and noting $u_{0}=u(\alpha)$ and $v_{0}=v(\alpha)$ their values at the edge of the confining box, the thermodynamical parameters of the self-gravitating gas are given by

$$% \eta=v_{0},$$ (64)

$$% \Lambda={1\over v_{0}}\biggl ({3\over 2}-u_{0}\biggr ),$$ (65)

$$% s\equiv {S-S_{0}\over Nk_{\rm B}}=-{1\over 2}\ln v_{0}-\ln (u_{0}v_{0})+v_{0}+2u_{0}-{3},$$ (66)

where

$$% {S_{0}\over Nk_{\rm B}}=\ln\mu+{1\over 2}\ln\pi-\ln 2-{1\over 2}\cdot$$ (67)

The free energy per particle $j=s+\eta\Lambda$ is

$$% j=-{1\over 2}\ln v_{0}-\ln(u_{0}v_{0})+v_{0}+u_{0}-{3\over 2}\cdot$$ (68)

We now expand the Milne variables close to the critical point  $\alpha_{\rm c}$. Recalling that $u(\alpha_{\rm c})=1$, $v(\alpha_{\rm c})=\eta_{\rm c}$ and $v'(\alpha_{\rm c})=0$, we get

 

$$% u_{0}=1+{1\over\alpha_{\rm c}}\left(2-\eta_{\rm c}\right)\epsil... ...ta_{\rm c}\right)\left(\eta_{\rm c}^{2}+2\eta_{\rm c}-4\right)\epsilon^{3}+...,$$ (69)

 

$$% v_{0}=\eta_{\rm c}+{\eta_{\rm c}\over 2\alpha_{\rm c}^{2}}(2-\e... ...m c}\over 6\alpha_{\rm c}^{3}}(2-\eta_{\rm c})(2+\eta_{\rm c})\epsilon^{3}+...,$$ (70)

where $\epsilon=\alpha-\alpha_{\rm c}$. We also recall that $\alpha_{\rm c}\simeq 8.993195...$ and $\eta_{\rm c}\simeq 2.517551...$ Substituting Eqs. (69) and (70) in Eq. (68), we find that

 

$$% j=j_{0}-{1\over 4\alpha_{\rm c}^{2}}(\eta_{\rm c}-2)\epsilon^{2}+{1\over 12\alpha_{\rm c}^{3}}(\eta_{\rm c}-2)(10-3\eta_{\rm c})\epsilon^{3}+...,$$ (71)

where $j_{0}=\eta_{\rm c}-(3/2)\ln\eta_{\rm c}-1/2$. On the other hand, recalling that $\eta=v(\alpha)$, we get

$$% \eta_{\rm c}-\eta={\eta_{\rm c}(\eta_{\rm c}-2)\over 2\alph... ...{\rm c}^{2}-4\right)\over 6\alpha_{\rm c}^{3}}\epsilon^{3}+...$$ (72)

Eliminating $\epsilon$ from the foregoing relations, we finally obtain

$$% j=j_{0}+(\eta-\eta_{\rm c})\biggl\lbrack K\pm\lambda(\eta_{\rm c}-\eta)^{1/2}\biggr\rbrack$$ (73)

with $K=1/2\eta_{\rm c}$ and $\lambda=(2/3)\sqrt{2(\eta_{\rm c}-2)/\eta_{\rm c}^{3}}$. We note that $K=\Lambda_{0}$ where $\Lambda_{0}$ is the normalized energy at the critical temperature (this comes from the fact that $\delta J=-E\delta\beta$). Equation (73) reproduces the cusp at $\eta=\eta_{\rm c}$ formed by the curve $j(\eta)$ in Fig. 11. Since ${\rm d}J=-E {\rm d}\beta$, $J(\alpha)$ and  $\eta(\alpha)$ are extrema at the same points in the series of equilibria, which is at the origin of the cusp. From Eq. (73), the barrier of free energy near the critical point  $\eta _{\rm c}$ is given by

$$% \Delta j=2\lambda(\eta_{\rm c}-\eta)^{3/2}.$$ (74)

Accordingly, the lifetime of the metastable state scales as

$$% t_{{\rm life}}\sim {\rm e}^{2\lambda N (\eta_{\rm c}-\eta)^{3/2}}$$ (75)

with $\lambda\simeq 0.16979815...$ Metastable states are therefore robust if $N(\eta_{\rm c}-\eta)^{3/2}\gg 1$. Note that if we had estimated the lifetime of the metastable state by $t_{{\rm life}}\sim {\rm e}^{\Delta E/k_{\rm B}T}$, we would have obtained $t_{{\rm life}}\sim {\rm e}^{2\lambda'' N (\eta_{\rm c}-\eta)^{1/2}}$ with $\lambda''=\lbrack 2(\eta_{\rm c}-2)/\eta_{\rm c}\rbrack^{1/2}\simeq 0.64121317...$ We see that entropic effects modify the power of  $\eta_{\rm c}-\eta$ in the expression of the metastable state lifetime. Returning to Eq. (75), the temperature of collapse $T_{\rm l}$ taking into account finite N effects can be estimated by $2\lambda N (\eta_{\rm c}-\eta_{\rm l})^{3/2}\sim 1$. This leads to

$$% \eta_{\rm l}=\eta_{\rm c}\biggl\lbrack 1-{1\over\eta_{\rm c}}\biggl ({1\over 2\lambda}\biggr )^{2/3}N^{-2/3}\biggr\rbrack.$$ (76)

A numerical application gives

$$% \eta_{\rm l}=2.517\ (1-0.816\ N^{-2/3}).$$ (77)

In the microcanonical ensemble, the lifetime of a metastable state can be estimated by

$$% t_{{\rm life}}\sim {\rm e}^{\Delta S}\sim {\rm e}^{N\Delta s},$$ (78)

where $\Delta S$ is the entropic barrier. By performing a study similar to the previous one close to the turning point of energy  $\Lambda_{\rm c}$ (see Fig. 13), we find that

$$% s=s_{0}+(\Lambda_{\rm c}-\Lambda)\biggl\lbrack \eta_{0}\pm \lambda'(\Lambda_{\rm c}-\Lambda)^{1/2}\biggr\rbrack,$$ (79)

with s₀=-0.192962..., $\eta_{0}=2.03085...$ and $\lambda'=0.863159...$ Therefore,

$$% \Delta s=2\lambda'(\Lambda_{\rm c}-\Lambda)^{3/2},$$ (80)

and

$$% t_{{\rm life}}\sim {\rm e}^{2\lambda'N(\Lambda_{\rm c}-\Lambda)^{3/2}}.$$ (81)

The energy of collapse $E_{\rm l}$ taking into account finite N effects can be estimated by $2\lambda'N(\Lambda_{\rm c}-\Lambda)^{3/2}\sim 1$. This leads to

$$% \Lambda_{\rm l}=\Lambda_{\rm c}\biggl\lbrack 1-{1\over\Lamb... ...\biggl ({1\over 2\lambda'}\biggr )^{2/3}N^{-2/3}\biggr\rbrack.$$ (82)

A numerical application gives

$$% \Lambda_{\rm l}=0.335\ (1-2.077\ N^{-2/3}).$$ (83)

This corresponds to a density contrast

$$% {\cal R}_{\rm l}=708.6\ (1-6.014\ N^{-1/3}).$$ (84)

These results are similar to those found by Katz & Okamoto (2000) by analyzing the microcanonical fluctuations of isothermal spheres. In particular, the scaling with N is the same. In the following section, we apply their theory of fluctuations to the canonical ensemble and show the relation with the preceding approach.

   
5.2 Canonical fluctuations in isothermal spheres

The canonical partition function can be written

$$% Z(\beta)=\int {\rm e}^{J(E)}{\rm d}E,$$ (85)

where $J(E)=S(E)-\beta E$. As before, we shall consider the situation where J(E) has two maxima (stable) and one minimum (unstable). This corresponds to the caloric curve of Fig. 11. We shall be interested here by the metastable state, i.e. the local maximum of J(E). Thus, we eliminate "by hands'' the contribution of the global maximum in the integral. Expanding J(E) around the local maximum up to second order, we obtain

$$% Z(\beta)={\rm e}^{J(E)}\int {\rm e}^{{1\over 2}J''(E)(\delta E)^{2}}{\rm d}(\delta E),$$ (86)

where E now refers to the energy of the metastable state and $\delta E$ is a fluctuation around equilibrium. According to the foregoing equation, the probability of the fluctuation $\delta E$ is given by

$$% P(\delta E)\sim {\rm e}^{{1\over 2}J''(E)(\delta E)^{2}}.$$ (87)

Noting that J''(E)=S''(E) and recalling that ${\partial S\over\partial E}=\beta(E)$ at equilibrium, we can rewrite the foregoing expression as

$$% P(\delta E)={1\over\sqrt{2\pi}}\biggl \vert{{\rm d}\beta\ov... ...\rm e}^{{1\over 2}{{\rm d}\beta\over {\rm d}E}(\delta E)^{2}}.$$ (88)

This formula shows that only equilibrium states with positive specific heats $C={\rm d}E/{\rm d}T=-\beta^{2}{{\rm d}E\over {\rm d}\beta}>0$ are stable in the canonical ensemble. The variance of the fluctuations of energy is given by the usual formula

$$% \left\langle (\delta E)^{2}\right\rangle=-{1\over {{\rm d}\beta\over {\rm d}E}}=T^{2}C\ge 0,$$ (89)

which can also be directly derived from the exact N-body distribution function (39) in a classical way. These results are valid on the whole gaseous branch of Fig. 11. If we now examine the situation close to the critical point ( $\beta_{\rm c}, E_{0}$) where $({\rm d}\beta/{\rm d}E)_{\rm c}=0$, we have to first order

$$% \biggl \vert{{\rm d}\beta\over {\rm d}E}\biggr \vert=\biggl... ... d}^{2}\beta\over {\rm d}E^{2}}\biggr \vert _{\rm c}(E-E_{0}).$$ (90)

We note E' the energy of the unstable saddle point of free energy J at temperature T. Close to the critical point ( $\beta_{\rm c}, E_{0}$) we can approximate the caloric curve $\beta(E)$ by a parabole. Thus, E' is related to E, the energy of the local maximum of free energy J at temperature T, by

 

E-E'=2(E-E₀). (91)

Using the criterion of Katz & Okamoto (2000), adapted to the canonical ensemble, we define the temperature of collapse as the temperature $T_{\rm l}$ at which the typical fluctuations of energy $\delta E=\sqrt{\langle (\delta E)^{2}\rangle}$ are of the same order as the difference E-E'. Indeed, as we approach the critical point ( $\beta_{\rm c}, E_{0}$) the fluctuations of energy become so important (since the specific heat diverges) that the system can overcome the barrier of potential played by the saddle point of free energy and collapse (eventually reaching the global maximum of J). Thus, for finite N systems, gravitational collapse can occur before the ending ( $\beta_{\rm c}, E_{0}$) of the metastable branch (spinodal point). According to the preceding criterion, the temperature of collapse is determined by the condition

$$% \langle (\delta E)^{2}\rangle_{\beta_{\rm l}}=4(E_{\rm l}-E_{0})^{2}.$$ (92)

Using Eqs. (89) and (90), we obtain

$$% E_{\rm l}=E_{0}+ \biggl \lbrace 4\biggl \vert{{\rm d}^{2}\beta\over {\rm d}E^{2}}\biggr \vert _{\rm c}\biggr \rbrace^{-1/3}.$$ (93)

It is more logical to write this criterion in terms of the temperature $T_{\rm l}$. Close to the critical point, we have

$$% \beta_{\rm c}-\beta={1\over 2}\biggl \vert{{\rm d}^{2}\beta\over {\rm d}E^{2}}\biggr \vert _{\rm c}(E-E_{0})^{2}.$$ (94)

Thus, using Eq. (93), we get

$$% \beta_{\rm l}=\beta_{\rm c}-{1\over 8} \biggl \lbrace 4\big... ...a\over {\rm d}E^{2}}\biggr \vert _{\rm c}\biggr \rbrace^{1/3}.$$ (95)

Figure 14: Caloric curve showing metastable states (local maxima of J) and unstable equilibria (saddle points of J) in the canonical ensemble for classical particles ( $\hbar \rightarrow 0$). As we approach the critical temperature  $\eta _{\rm c}$, the fluctuations of energy increase considerably allowing the system to collapse before the end of the metastable branch (spinodal point). The temperature of collapse $T_{\rm l}$ is the temperature at which the typical fluctuations of energy are of the same order as the energy difference  $\vert\epsilon '-\epsilon \vert$.

Open with DEXTER

We can obtain more explicit results in the case of classical isothermal spheres (see Fig. 14). Introducing the usual dimensionless parameters $\Lambda=-ER/GM^{2}$ and $\eta=\beta GMm/R$, the foregoing relation can be rewritten

$$% \eta_{\rm l}=\eta_{\rm c}-{1\over 8N^{2/3}} \biggl \lbrace ... ... {\rm d}\Lambda^{2}}\biggr \vert _{\rm c}\biggr \rbrace^{1/3},$$ (96)

Now, we have (see Chavanis 2002a)

$$% \biggl ({{\rm d}^{2}\eta\over {\rm d}\Lambda^{2}}\biggr )_{\rm c}=-{\eta_{\rm c}^{3}\over\eta_{\rm c}-2}\cdot$$ (97)

Therefore, the temperature of collapse for finite N systems is given by

$$% \eta_{\rm l}=\eta_{\rm c}\biggl\lbrack 1-{1\over 8N^{2/3}}\biggl ({4\over \eta_{\rm c}-2}\biggr )^{1/3}\biggr\rbrack \cdot$$ (98)

A numerical application gives

$$% \eta_{\rm l}=2.517\ (1-0.247\ N^{-2/3}).$$ (99)

This estimate can be compared with Eq. (77) which has the same scaling with N. Of course, we should not give too much credit on the numerical factor in front of N^-2/3 since the criterion for determining $T_{\rm l}$ is essentially phenomenological. The corresponding energy  $\Lambda_{\rm l}$ can be deduced from Eq. (94) and is given by

$$% \Lambda_{\rm l}=\Lambda_{0}\biggl\lbrack 1- 2\biggl ({\eta_{\rm c}-2\over 4}\biggr )^{1/3}N^{-1/3}\biggr \rbrack.$$ (100)

It may also be of interest to determine the corresponding density contrast  ${\cal R}_{\rm l}$. To that purpose, we start from the formula (see Chavanis 2002a)

$$% {\cal R}={\rm e}^{\psi(\alpha)}={\alpha^{2}\over u_{0}v_{0}},$$ (101)

and we use the expansion (69) and (70) of the Milne variables. This yields

$$% {\cal R}_{\rm l}={\cal R}_{\rm c}\biggl (1+{\eta_{\rm c}\over\alpha_{\rm c}}\epsilon_{l}\biggr ),$$ (102)

where ${\cal R}_{\rm c}=\alpha_{\rm c}^{2}/\eta_{\rm c}$. Now, $\epsilon_{\rm l}$ is determined by Eqs. (72) and (98) yielding

$$% \epsilon_{\rm l}={-\alpha_{\rm c}\over 4^{1/3}N^{1/3}(\eta_{\rm c}-2)^{2/3}}\cdot$$ (103)

Substituting this result in Eq. (102) we finally obtain

$$% {\cal R}_{\rm l}={\cal R}_{\rm c}\biggl \lbrace 1-{\eta_{\r... ... \lbrack 2(\eta_{\rm c}-2)\rbrack^{2/3}}N^{-1/3}\biggr\rbrace,$$ (104)

or, with numerical values,

$$% {\cal R}_{\rm l}=32.1\ (1-2.45\ N^{-1/3}).$$ (105)

The above theory thus predicts that, for finite N systems, the collapse should take place slightly before the canonical spinodal point  $\eta _{\rm c}$ due to the enhancement of energy fluctuations as we approach this critical point. Similar conclusions can be obtained by developing a theory of hydrodynamic fluctuations showing that the density fluctuations become large before the point of ordinary instability is reached (Monaghan 1978). The Monte Carlo simulations of de Vega & Sanchez (2002) in the canonical ensemble (with N=2000particles) reveal that the collapse indeed takes place before the critical point. However, the collapse occurs apparently at the point where the isothermal compressibility $\kappa_{T}=-(1/V)(\partial V/\partial p)_{T}$ diverges. This corresponds to an inverse normalized temperature $\eta_{T}=2.43450...$ It is sensibly smaller than the value obtained from Eq. (99) with N=2000. The same discrepency with the prediction of Katz & Okamoto (2000) is found in the microcanonical ensemble by de Vega & Sanchez (2002). These results seem to indicate that the higher collapse temperature and energy found in Monte Carlo simulations are not due to finite N effects. They seem to be independent on N and correspond to other critical points that do not coincide with the spinodal point. We intend to perform independent Monte Carlo simulations to check these results.

   
5.3 General expression of the potential barrier

We can use the preceding approach to obtain a simple approximate expression of the potential barrier close to the critical point  $\eta _{\rm c}$. Consider a system at fixed temperature T and denote by E its equilibrium energy such that J'(E)=0 (we consider here that E is the energy of the metastable state). For a fluctuation  $E+\delta E$, we have seen that the variation of free energy $\delta J={1\over 2}J''(E)(\delta E)^{2}$ can be expressed as

$$% \delta J={1\over 4}\beta''(E_{0})(E-E')(\delta E)^{2}$$ (106)

where we have assumed that E is close to the critical point E₀and, as before, E' is the energy of the minimum of J at temperature T. We can use this expression to estimate the potential barrier $\Delta J=J_{{\rm local}}-J_{{\rm saddle}}$. Thus, setting $\delta E=E-E'$, we get

$$% \Delta J={1\over 4}\vert\beta''(E_{0})\vert(E-E')^{3}.$$ (107)

Using Eq. (94) to express this relation in terms of the temperature, we finally obtain

$$% \Delta J=\sqrt{32\over \vert\beta''(E_{0})\vert}(\beta_{\rm c}-\beta)^{3/2} \qquad ({\rm approx.}).$$ (108)

This is an estimate, because the curve J(E) is not just a parabole between E and E' as it has been assumed to get Eq. (106) from Eq. (91). Using Eq. (97), this approximate expression can be written

$$% \Delta j=6\lambda (\eta_{\rm c}-\eta)^{3/2} \qquad ({\rm approx.}).$$ (109)

It differs from the exact expression (74) by a factor 3. Noting that $\Delta J={1\over 2}(\Delta E)^{2}/\langle (\delta E)^{2}\rangle$, according to Eqs. (88) and (89), the criterion (92) of Katz & Okamoto (2000) discussed previously is equivalent to $\Delta J\sim 1/2$. Alternatively, writing $t_{{\rm life}}\sim {\rm e}^{\Delta J}$, the criterion that we have used to get Eq. (76) corresponds to $\Delta J\sim 1$. Therefore, on a qualitative point of view, the approaches of Sects. 5.1 and 5.2 are equivalent. They slightly differ on the details (definition of collapse temperature, estimate of $\Delta J$) explaining why Eqs. (77) and (99) are not exactly the same.

We can also try to calculate $\Delta J$ by working directly on the series of equilibria $\beta(E)$. Taking E as a control parameter, we have $J(E)=S(E)-\beta(E)E$. Noting that $J'(E)=-\beta'(E)E$, $J''(E)=-\beta''(E)E-\beta'(E)$ and $J'''(E)=-\beta'''(E)E-2\beta''(E)$, and expanding J(E) close to the critical point where $\beta'(E_{0})=0$, we get

 

$$% J(E)=J(E_{0})-{1\over 2}\beta''(E_{0})E_{0}(E-E_{0})^{2} -{1\over 3!}\lbrack \beta'''(E_{0})E_{0}+2\beta''(E_{0})\rbrack (E-E_{0})^{3}+...$$ (110)

Using the relation

 

$$% \beta-\beta_{\rm c}= {1\over 2}\beta''(E_{0})(E-E_{0})^{2}+{1\over 3!}\beta'''(E_{0})(E-E_{0})^{3}+...$$ (111)

between the temperature and the energy close to the critical point, we find that

 

$$% J(E)=J(E_{0})-E_{0}(\beta-\beta_{\rm c}) \pm {1\over 3} \beta''... ...iggl\lbrack {2(\beta-\beta_{\rm c})\over \beta''(E_{0})}\biggr\rbrack^{3/2}+...$$ (112)

Therefore, close to the critical point, the barrier of free energy is exactly given by

$$% \Delta J={1\over 3}\sqrt{32\over \vert\beta''(E_{0})\vert}(\beta_{\rm c}-\beta)^{3/2} \qquad ({\rm exact}),$$ (113)

which returns Eq. (74).

We can use the same type of approach in the microcanonical ensemble to obtain a simple expression of the entropic barrier close to the critical point  $(\Lambda_{\rm c},\eta_0)$. Consider a system at fixed energy E and denote by $\beta$ its equilibrium temperature (we consider here that $\beta$ is the temperature of the metastable state). For a fluctuation  $\beta+\delta \beta$, the variation of entropy can be expressed as (see Katz & Okamoto 2000 for details)

$$% \delta S={1\over 2}E'(\beta)(\delta\beta)^{2}.$$ (114)

Assuming that $\beta$ is close to the critical point $\beta_{0}$, and using arguments similar to those developed previously, we get

$$% \delta S={1\over 4}E''(\beta_{0})(\beta-\beta')(\delta\beta)^{2}$$ (115)

where $\beta'$ is the inverse temperature of the saddle point of entropy. We can use this expression to estimate the entropic barrier $\Delta S=S_{{\rm local}}-S_{{\rm saddle}}$. Thus, setting $\delta \beta=\beta-\beta'$, we find that

$$% \Delta S={1\over 4}E''(\beta_{0})(\beta-\beta')^{3}.$$ (116)

Using

$$% E-E_{\rm c}\simeq {1\over 2}E''(\beta_{0})(\beta-\beta_{0})^{2},$$ (117)

to express the temperature as a function of the energy close to the critical point, we finally obtain

$$% \Delta S=\sqrt{32\over E''(\beta_{0})}(E-E_{\rm c})^{3/2} \qquad ({\rm approx.}).$$ (118)

We note that Eqs. (108) and (118) are symmetrical provided that we interchange E and $\beta$. Evaluating numerically ${\rm d}^{2}\Lambda/{\rm d}\eta^{2}$ at the critical point  $(\Lambda_{\rm c},\eta_0)$, this approximate expression can be written

$$% \Delta s=6\lambda' (\Lambda-\Lambda_{\rm c})^{3/2} \qquad ({\rm approx.}).$$ (119)

It differs from the exact expression (80) by a factor 3.

We can also try to calculate $\Delta S$ directly from the series of equilibria $E(\beta)$. Taking $\beta$ as a control parameter, we have $S(\beta)=J(\beta)+\beta E(\beta)$ and we recall that ${\rm d}S=\beta {\rm d}E$and ${\rm d}J=-E {\rm d}\beta$. Thus $S'(\beta)=\beta E'(\beta)$, $S''(\beta)=E'(\beta)+\beta E''(\beta)$ and $S'''(\beta)=2E''(\beta)+\beta E'''(\beta)$. Expanding $S(\beta)$close to the critical point where  $E'(\beta_{0})=0$, we get

 

$$% S(\beta)=S(\beta_{0})+{1\over 2}E''(\beta_{0})\beta_{0}(\beta-\... ...brack 2E''(\beta_{0})+\beta_{0}E'''(\beta_{0})\rbrack (\beta-\beta_{0})^{3}+...$$ (120)

Using the relation

 

$$% E-E_{\rm c}= {1\over 2}E''(\beta_{0})(\beta-\beta_{0})^{2}+{1\over 3!}E'''(\beta_{0})(\beta-\beta_{0})^{3}+...$$ (121)

between the energy and the temperature close to the critical point, we find that

 

$$% S(\beta)\!=\!S(\beta_{0})\!+\!\beta_{0}(E-E_{\rm c})\pm {1\over... ...t)\biggl\lbrack {2(E\!-\!E_{\rm c})\over E''(\beta_{0})}\biggr\rbrack^{3/2}+...$$ (122)

Therefore, close to the critical point, the barrier of entropy is exactly given by

$$% \Delta S={1\over 3}\sqrt{32\over E''(\beta_{0})}(E-E_{\rm c})^{3/2} \qquad ({\rm exact}),$$ (123)

which returns Eq. (80).

   
6 Relation to the Kramers problem

   
6.1 The Fokker-Planck equation

In the preceding section, we have used the Kramers formula to estimate the lifetime of metastable states in self-gravitating systems. We would like now to justify this formula from first principles. In order to determine the lifetime of a metastable state, we need to introduce a dynamical model. In the canonical ensemble, we can consider a system of self-gravitating Brownian particles (Chavanis et al. 2002) described by the stochastic equations

$$% {{\rm d}{\vec r}_{i}\over {\rm d}t}={\vec v}_{i},$$ (124)

$$% {{\rm d}{\vec v}_{i}\over {\rm d}t}=-\nabla_{i}U({\vec r}_{1},...,{\vec r}_{N})-\xi{\vec v}_{i}+\sqrt{2D}{\vec R}_{i}(t),$$ (125)

where $-\xi {\vec v}_{i}$ is a friction force and ${\vec R}_{i}(t)$ is a white noise satisfying $\langle {\vec R}_{i}(t)\rangle={\bf0}$ and $\langle {R}_{a,i}(t){R}_{b,j}(t')\rangle=\delta_{ij}\delta_{ab}\delta(t-t')$, where a,b=1,2,3 refer to the coordinates of space and i,j=1,...,Nto the particles. The particles interact via the gravitational potential $U({\vec r}_{1},...,{\vec r}_{N})=\sum_{i<j}u({\vec r}_{i}-{\vec r}_{j})$ where $u({\vec r}_{i}-{\vec r}_{j})=-G/\vert{\vec r}_{i}-{\vec r}_{j}\vert$. The inverse temperature $\beta=1/T$ is related to the diffusion coefficient through the Einstein relation $\xi=D\beta$.

Using standard stochastic processes, we can derive the N-body Fokker-Planck equation (see Chavanis 2004b)

 

$$% {\partial P_{N}\over\partial t}+\sum_{i=1}^{N}\biggl ({\vec v}_... ... {\partial P_{N}\over\partial {\vec v}_{i}}+\xi P_{N}{\vec v}_{i}\biggr\rbrack,$$ (126)

where ${\vec F}_{i}=-\nabla_{i}U({\vec r}_{1},...,{\vec r}_{N})$ is the gravitational force acting on particle i and  $P_{N}({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N},t)$ is the N-body distribution function. Its stationary states correspond to the canonical distribution

 

$$% P_{N}({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N})=... ...e \sum_{i=1}^{N}{v_{i}^{2}\over 2}+U({\bf r}_{1},...,{\bf r}_{N})\bigr\rbrace}.$$ (127)

If we implement a mean-field approximation (see Chavanis 2004b), we can show that the distribution function $f({\vec r},{\vec v},t)=NP_{1}$ is solution of the Kramers-Poisson system. However, this is not the approach that we shall consider here.

We wish to obtain the time evolution of the distribution of energies P(E,t). To that purpose, we shall follow a method similar to that developed by Kramers (1940) in his investigation of the escape of Brownian particles over a potential barrier. The difference is that we work here in a 6N dimensional phase space. Assuming that $P_{N}({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N},t)$ depends only on energy $E=\sum_{i=1}^{N}{v_{i}^{2}\over 2}+U({\vec r}_{1},...,{\vec r}_{N})$ and time t, and averaging the Fokker-Planck Eq. (126) on the hypersurface of energy E, we show in Appendix B that

$$% g(E){\partial P_{N}\over\partial t}(E,t)=3M{\partial\over\p... ...rtial P_{N}\over\partial E}+\beta P_{N}\biggr )\biggr \rbrack,$$ (128)

where g(E) is the density of states and I(E) is the phase space hypervolume with energy less than E (thus $g(E)={\rm d}I/{\rm d}E$). Now, the distribution of energies is given by

 

P(E,t)=P_N(E,t)g(E). (129)

At equilibrium, using Eq. (127), we have

$$% P(E)={1\over Z(\beta)}g(E){\rm e}^{-\beta E},$$ (130)

which returns Eq. (53). Out of equilibrium, substituting Eqs. (129) into (128) and simplifying the resulting expressions, we finally obtain

$$% {\partial P\over\partial t}={\partial\over\partial E}\biggl... ...\partial P\over\partial E}+\beta PF'(E)\biggr )\biggr \rbrack,$$ (131)

where D(E)=3M I(E)/g(E) and $F(E)=E-TS(E)=E-T\ln g(E)$ is the free energy. This is similar to the Fokker-Planck equation describing the stochastic motion of a particle in a potential where the energy E plays the role of the position x and where the free energy F(E) plays the role of the potential U(x). In the following, we shall assume that the free energy F(E) has a local minimum at E_A (metastable), a local maximum at E_B (unstable) and a global minimum at $E_{\rm C}$. A typical situation is illustrated in Fig. 15. We shall prepare a large number ${\cal N}\gg 1$of systems close to the energy E_A with the canonical distribution (130). Thus, ${\cal N}{\times} P(E,t) {\rm d}E$ gives the number of systems with energy between E and  $E+{\rm d}E$ at time t. As time goes on, a fraction of these systems reaches the energy E_B and undergoes gravitational collapse towards $E_{\rm C}$. Therefore, we adopt the boundary condition

 

P(E_B,t)=0. (132)

Our aim, now, is to estimate the current of diffusion past E_B and the typical lifetime of metastable states.

	Figure 15: Illustration of the barrier of free energy corresponding to Fig. 12. Point A corresponds to the gaseous metastable state and point B is the unstable solution creating a barrier for reaching the condensed state C.
Open with DEXTER

   
6.2 The stationary solutions

The stationary solutions of Eq. (131) are of the form

$$% {\partial P\over\partial E}+\beta PF'(E)=-{J\over D(E)},$$ (133)

where J<0 is the current of diffusion in energy space. Using the boundary condition (132), the solution of Eq. (133) reads

$$% P(E)=J {\rm e}^{-\beta F(E)}\int_{E}^{E_{B}}{{\rm e}~^{\beta F(x)}\over D(x)}{\rm d}x.$$ (134)

The current of diffusion can therefore be expressed as

$$% J={P_{A}{\rm e}~^{\beta F(E_{A})}\over\int_{E_{A}}^{E_{B}}{{\rm e}~^{\beta F(E)}\over D(E)}{\rm d}E}\cdot$$ (135)

To estimate the probability P_A, we shall approximate the curve F(E) close to A by a parabole. We thus make the expansion

$$% F(E)=F(E_{A})+{1\over 2TC_{A}}(E-E_{A})^{2}+...$$ (136)

where we have used $F''(E)=-T\beta'(E)=1/TC$ where $C={\rm d}E/{\rm d}T$ is the specific heat. Therefore,

 

$${1\over Z}{\rm e}^{-\beta F(E_{A})}$$                              P_A =

$$\simeq {1\over\int_{-\infty}^{+\infty}{\rm e}^{-{\beta^{2}\over 2C_{A}}(E-E_{A})^{2}}{\rm d}E}={\beta\over\sqrt{2\pi C_{A}}}\cdot$$ (137)

On the other hand, the integral in Eq. (135) is dominated by the value of the integrand close to B. Making the same quadratic expansion as in Eq. (136), we get

 

$$% \int_{E_{A}}^{E_{B}}{{\rm e}~^{\beta F(E)}\over D(E)}{\rm d}E \simeq -{1\over 2}{\times} {{\rm e}~^{\beta F(E_{B})}\over D(E_{B})}\int_{-\infty}^{+\infty}{\rm e}^{{\beta^{2}\over 2C_{B}}(E-E_{B})^{2}}{\rm d}E$$

$$-{{\rm e}~^{\beta F(E_{B})}\over 2 D(E_{B})}{\sqrt{2\pi \vert C_{B}\vert}\over\beta},$$ = (138)

where we recall that C_B<0 for the unstable solution. We thus obtain the expression of the current

$$% J=-{\beta^{2}D(E_{B})\over \pi \sqrt{C_{A}\vert C_{B}\vert}} {\rm e}^{-\beta (F_{B}-F_{A})}.$$ (139)

This expression involving the barrier of free energy $\Delta F$, is similar to the one obtained by Kramers (1940) in his classical study. In our case, the parameters have a thermodynamical interpretation while Kramers considers a dynamical system in a potential U(x).

   
6.3 The escape time

The preceding approach assumes that the population of systems that are introduced at A is continuously renewed so as to counterbalance the population of systems that are lost at B and maintain a stationary regime. We shall now relax this simplifying assumption and look for decaying solutions of Eq. (131) of the form

$$% P(E,t)={\rm e}^{-\lambda t}h(E),$$ (140)

where h(E) satisfies the differential equation

$$% {{\rm d}\over {\rm d}E}\biggl\lbrack D\biggl ({{\rm d}h\over {\rm d}E}+\beta h F'(E)\biggr )\biggr\rbrack=-\lambda h.$$ (141)

Assuming that close to E_A the system is at equilibrium, the foregoing equation can be integrated into

$$% {{\rm d}h\over {\rm d}E}+\beta h F'(E)=-{\lambda \over D(E)}\int_{E_{A}}^{E} h(E'){\rm d}E'.$$ (142)

In usual situations, the eigenvalue $\lambda$ is expected to be small as it corresponds to the inverse lifetime of the metastable states. We thus consider the perturbative expansion of h in powers of $\lambda$ and write

$$% h=h_{0}+\lambda h_{1}+...$$ (143)

Substituting this expansion in Eq. (142) and identifying terms of equal order, we obtain the differential equations

$$% {{\rm d}h_{0}\over {\rm d}E}+\beta h_{0} F'(E)=0,$$ (144)

$$% {{\rm d}h_{1}\over {\rm d}E}+\beta h_{1} F'(E)=-{1 \over D(E)}\int_{E_{A}}^{E} h_{0}(E'){\rm d}E'.$$ (145)

The first equation integrates into

$$% h_{0}=K{\rm e}^{-\beta F(E)}.$$ (146)

Substituting this result in Eq. (145), we get

$$% {{\rm d}h_{1}\over {\rm d}E}+\beta h_{1} F'(E)=-{K \over D(E)}\int_{E_{A}}^{E} {\rm e}^{-\beta F(E')} {\rm d}E'.$$ (147)

The solution of this first order differential equation can be written

$$% h_{1}=-\chi(E)K{\rm e}^{-\beta F(E)},$$ (148)

where the function $\chi$ is defined by

$$% \chi'(E)={{\rm e}~^{\beta F(E)}\over D(E)}\int_{E_{A}}^{E}{\rm e}^{-\beta F(E')}{\rm d}E',$$ (149)

with $\chi(E_{A})=0$. Therefore, in the approximation $\lambda\ll 1$, the solution of Eq. (131) is

$$% P(E,t)=K{\rm e}^{-\lambda t} {\rm e}^{-\beta F(E)}\bigl\lbrack 1-\lambda \chi(E)\bigr\rbrack.$$ (150)

The eigenvalue $\lambda$ is determined by the boundary condition (132) yielding

$$% \lambda={1\over\chi(E_{B})}\cdot$$ (151)

Therefore, the lifetime of the metastable state is given by

$$% t_{{\rm life}}\sim\lambda^{-1}= \chi(E_{B}).$$ (152)

We can now try to simplify this expression. First, we approximate Eq. (149) by

 

$$% \chi'(E) -{{\rm e}~^{\beta F(E)}\over D(E)}\int_{-\infty}^{+\infty}{\rm e}^{-\beta F(E_{A})}{\rm e}^{{-\beta^{2}\over 2C_{A}}(E-E_{A})^{2}}{\rm d}E$$ =

$$-{\sqrt{2\pi C_{A}}\over \beta}{\rm e}^{-\beta F(E_{A})}{{\rm e}~^{\beta F(E)}\over D(E)}\cdot$$ = (153)

After integrating and setting E=E_B, we obtain

$$% \chi(E_{B})=-{\sqrt{2\pi C_{A}}\over \beta}{\rm e}^{-\beta ... ...t_{E_{A}}^{E_{B}}{{\rm e}~^{\beta F(E')}\over D(E')}{\rm d}E'.$$ (154)

With the additional approximation

 

$$% \chi(E_{B})={\sqrt{2\pi C_{A}}\over \beta}{\rm e}^{-\beta F(E_{... ...er D(E_{B})}{\rm e}^{-{\beta^{2}\over 2\vert C_{B}\vert}(E-E_{B})^{2}}{\rm d}E,$$ (155)

we finally get

 

$$% \lambda={\beta^{2}D(E_{B})\over \pi \sqrt{C_{A}\vert C_{B}\vert}} {\rm e}^{-\beta (F_{B}-F_{A})}.$$ (156)

Therefore, the lifetime of a metastable state behaves as

 

$$% t_{{\rm life}}\sim {\rm e}^{\Delta F\over kT}\sim {\rm e}^{\Delta J}.$$ (157)

We note that the expression of $\lambda$ is similar to the expression (139) obtained for the current J. The connexion is the following. In the non-stationary case, the current of diffusion at E_B is

 

$$-D(E_{B})\biggl ({\partial P\over\partial E}+\beta PF'(E)\biggr )_{E_{B}}$$                                J_B =

$$-\lambda {\rm e}^{-\lambda t} D(E_{B})\biggl ({\partial h_{1}\over\partial E}+\beta h_{1}F'(E)\biggr )_{E_{B}}$$ =

$$\lambda {\rm e}^{-\lambda t} \int_{E_{A}}^{E_{B}} h_{0}(E'){\rm d}E'\simeq \lambda \int_{E_{A}}^{E_{B}} P(E',t){\rm d}E'.$$ = (158)

Hence, normalizing the current by the exponential decay of the density probability, we get

 

$$% {J_{B}\over {\rm e}^{-\lambda t}}=-\lambda,$$ (159)

which is equivalent to Eq. (139).

In the preceding analysis, we have worked in the canonical ensemble because the Brownian model (124) and (125) is easier to study than the N-stars Hamiltonian model, while exhibiting qualitatively the same phenomena (phase transitions, metastable states etc.). We expect to have symmetric expressions in the microcanonical ensemble with the correspondance $E\leftrightarrow \beta$ and $S\leftrightarrow J$. This study is left for a future work.

   
7 Conclusion

In this paper, we have completed previous investigations concerning the statistical mechanics of self-gravitating systems in microcanonical and canonical ensembles. The microcanonical ensemble is the proper description of isolated Hamiltonian systems such as globular clusters (Binney & Tremaine 1987). The canonical ensemble is relevant for systems in contact with a heat bath of non-gravitational origin. It is also the proper description of stochastically forced systems such as self-gravitating Brownian particles (Chavanis et al. 2002). We have justified the mean-field approximation, in a proper thermodynamic limit $N\rightarrow +\infty$ with $\eta=\beta GMm/R$ and $\epsilon=ER/GM^{2}$ fixed, from the equilibrium BBGKY hierarchy (Appendix A). In this thermodynamic limit, the equilibrium state is determined by a maximization problem: the maximization of entropy at fixed mass and energy in the microcanonical ensemble and the minimization of free energy at fixed mass and temperature in the canonical ensemble. This determines the most probable macroscopic distribution of particles at equilibrium. This can also be seen as a saddle point approximation in the functional integral formulation of the density of states and partition function. The probability of the distribution function is given by a Cramer formula ${\cal P}[f]\sim {\rm e}^{Ns[f]}$ or ${\cal P}[f]\sim {\rm e}^{Nj[f]}$, which constitutes a result of Large Deviations (for a regularized gravitational potential). We have shown that the saddle point approximation is less and less accurate close to the transition point since the condition $N\vert E-E_{\rm t}\vert\gg 1$ (in microcanonical ensemble) or $N\vert T-T_{\rm t}\vert\gg 1$ (in canonical ensemble) must be satisfied.

We have also argued that the lifetime of metastable states (local entropy maxima) scales as ${\rm exp}(N)$ due to the long-range nature of the interaction. Therefore, the importance of these metastable states is considerable and they cannot be simply ignored. Metastable states are in fact stable and they correspond to observed structures in the universe such as globular clusters. The preceding estimate must, however, be revised close to the critical point. By solving a Fokker-Planck equation, we have shown that the lifetime of metastable states is given by the Kramers formula involving the barrier of entropy or free energy. These barriers have been calculated exactly close to the Antonov energy $E_{\rm c}$ (in microcanonical ensemble) and close to the Jeans-Emden temperature $T_{\rm c}$ (in canonical ensemble). We have obtained the estimates $t_{{\rm life}}\sim {\rm exp}\lbrace 1.726\ N (\Lambda_{\rm c}-\Lambda)^{3/2}\rbrace$ (in microcanonical ensemble) and $t_{{\rm life}}\sim {\rm exp}\lbrace 0.339\ N (\eta_{\rm c}-\eta)^{3/2}\rbrace$ (in canonical ensemble) so that the lifetime decreases as we approach $E_{\rm c}$ or $T_{\rm c}$. This implies that the collapse will take place slightly above $E_{\rm c}$ or $T_{\rm c}$ at an energy $\Lambda_{\rm l}=\Lambda_{\rm c}(1-2.077\ N^{-2/3})$ or temperature $\eta_{\rm l}=\eta_{\rm c}(1-0.816\ N^{-2/3})$. Similar conclusions have been reached by Katz & Okamoto (2000). Yet, these predictions do not seem to be consistent with the Monte Carlo simulations of de Vega & Sanchez (2002), although they find that the collapse indeed takes place slightly before the critical point. Independent simulations are under preparation to check that point. In a recent work (Chavanis 2004b; and Appendix A), we have shown that the mean-field approximation becomes inaccurate close to the critical point. This may be a reason for the discrepency between the theory of fluctuations exposed here and the numerical simulations of de Vega & Sanchez (2002).

Finally, a part of our discussion was devoted to answer the critics raised by Gross (2003, 2004) in recent comments. This author argues that the microcanonical entropy  $S_{{\rm micro}}(E)$ and the microcanonical temperature  $\beta _{{\rm micro}}(E)$ must be single valued. This is true in a strict sense, but the problem is richer than that because of the existence of long-lived metastable states. Therefore, the physical caloric curve/series of equilibria $\beta(E)$ is multi-valued and leads to "dinosaur's necks'' and special "microcanonical phase transitions'' (Chavanis 2002b). This is specific to systems with long-range interactions in view of the long lifetime of metastable states (local entropy maxima). These results have stimulated a general classification of phase transitions by Bouchet & Barré (2005). Microcanonical phase transitions (as in Fig. 1) have not been fully appreciated by Gross and his collaborators because their studies (e.g., Votyakov et al. 2002) consider a large small-scale cut-off for which the caloric curve looks like Fig. 2 and is univalued. If these authors reduce their small-scale cut-offs, they will see "dinosaurs'' appear!

On a mathematical point of view, the statistical mechanics of the classical self-gravitating gas ($\hbar=0$) poses problem because the density of states diverges even if the system is enclosed within a box, due to the formation of binaries. Our approach suggests the following procedure to rigorously carry out a statistical mechanics: first, introduce a soften potential with small-scale cut-off $\epsilon$. We can then perform the statistical mechanics properly. Keep only the metastable states defined as local maxima of entropy. For $N\rightarrow +\infty$, they are physically relevant since they have infinite lifetime. Finally, take the limit $\epsilon\rightarrow 0$ (the structure of the gaseous metastable states is independent on the cut-off for $\epsilon\ll 1$). We note that the order of the limits $N\rightarrow +\infty$ and $\epsilon\rightarrow 0$ is not interchangeable. It should be made clear that these metastable states exist only in "boxes'', or for truncated isothermal systems, at less than the Antonov critical radius (this concerns about $80\%$ of tidally truncated globular clusters). Now, real star clusters do not leave in boxes so it is always the evolution beyond the Antonov point (core collapse) that matters eventually especially as the mass segregation decreases the Antonov radius. This ultimate evolution leads to the formation of binaries.

Acknowledgements

I am grateful to J. Katz for stimulating discussions. I also thank the referee for additional comments. While this paper was in course of redaction, I became aware of a preprint by Antoni et al. [cond-mat/0401177] (Europhys. Lett., 66, 645 (2004)) where similar arguments about the lifetime of metastable states for a toy model with long-range interactions have been developed independently.

   
Appendix A: Justification of the mean-field approximation from the equilibrium BBGKY hierarchy

In this Appendix, we show that the mean-field approximation is exact for self-gravitating systems in a properly defined thermodynamic limit $N\rightarrow +\infty$ with $\eta=\beta GMm/R$ and $\Lambda=-ER/GM^{2}$ fixed. In the canonical ensemble, the equilibrium N-body distribution function is given by

$$% P_{N}={1\over Z(\beta)}{\rm e}^{-\beta U({\vec r}_{1},...,{\vec r}_{N})},$$ (A.1)

where we only consider the configurational part (the velocity part, which is just a product of Maxwellians, is trivial). Here, $U({\vec r}_{1},...,{\vec r}_{N})=\sum_{i<j}u_{ij}$ where $u_{ij}=-G/\vert{\vec r}_{i}-{\vec r}_{j}\vert$ is the gravitational potential (rigorously speaking, we must use a soften potential in order to regularize the partition function). Taking the derivative of Eq. (A.1) with respect to  ${\vec r}_1$, we get

$$% {\partial P_{N}\over\partial {\vec r}_{1}}=-\beta P_{N}{\partial U\over\partial {\vec r}_{1}}\cdot$$ (A.2)

From this relation we can obtain the full equilibrium BBGKY hierarchy for the reduced distribution functions (Chavanis 2004b). Restricting ourselves to the one and two-body distribution functions

$$% P_{j}({\vec r}_{1},...,{\vec r}_{j})=\int P_{N}({\vec r}_{1... ...ec r}_{N}){\rm d}^{3}{\vec r}_{j+1}...{\rm d}^{3}{\vec r}_{N},$$ (A.3)

with j=1,2, we find

$$% {\partial P_{1}\over\partial {\vec r}_{1}}({\vec r}_{1})=-\... ...tial u_{12}\over\partial {\vec r}_{1}}{\rm d}^{3}{\vec r}_{2},$$ (A.4)

 

$$% {\partial P_{2}\over\partial {\vec r}_{1}}({\vec r}_{1},{\vec r... ...\vec r}_{3}){\partial u_{13}\over\partial {\vec r}_{1}}{\rm d}^{3}{\vec r}_{3}.$$ (A.5)

As is well-known, each equation of the hierarchy involves the next order distribution function. We now decompose the two-body and three-body distribution functions in the suggestive forms

$$% P_{2}({\vec r}_{1},{\vec r}_{2}) = P_{1}({\vec r}_{1})P_{1}({\vec r}_{2})+P_{2}'({\vec r}_{1},{\vec r}_{2}),$$ (A.6)

 

$$P_{3}({\vec r}_{1},{\vec r}_{2},{\vec r}_{3})=P_{1}({\vec r}_{1})... ...\vec r}_{3})P_{1}({\vec r}_{1})+P_{3}'({\vec r}_{1},{\vec r}_{2},{\vec r}_{3}),$$ (A.7)

where P_n' are the cumulants. We shall consider the thermodynamic limit $N\rightarrow +\infty$ with fixed $\eta=\beta GNm^{2}/R$. In this limit, it can be shown that the non trivial correlations P'_nare of order  N^-(n-1). Here, we shall just establish this result for the two-body distribution function P'₂. Substituting the decompositions (A.6) and (A.7) in Eqs. (A.4) and (A.5) and assuming that P'₃, of order N^-2, is negligible (this corresponds to the Kirkwood approximation in plasma physics) the equation for the two-body distribution function becomes after simplification

 

$${\partial P_{2}'\over\partial {\vec r}_{1}}({\vec r}_{1},{\vec r}... ...{1}({\vec r}_{1})P_{1}({\vec r}_{2}){\partial u_{12}\over\partial {\vec r}_{1}}$$

$$\qquad\qquad -\beta P'_{2}({\vec r}_{1},{\vec r}_{2}){\partial u_... ...{\vec r}_{3}){\partial u_{13}\over\partial {\vec r}_{1}}{\rm d}^{3}{\vec r}_{3}$$

$$\qquad\qquad -\beta NP_{1}({\vec r}_{1})\int P'_{2}({\vec r}_{2},{\vec r}_{3}){\partial u_{13}\over\partial {\vec r}_{1}}{\rm d}^{3}{\vec r}_{3}.$$ (A.8)

We thus find that $P_{1}\sim 1$ and $P_2'\sim\beta u\sim \beta Gm^{2}/R= \eta/N=O(1/N)$. Therefore, in the limit $N\rightarrow +\infty$, the two-body distribution function is the product of two one-body distribution functions:

$$% P_2({\vec r}_1,{\vec r}_2)=P_1 ({\vec r}_1)P_1 ({\vec r}_2)+O(1/N).$$ (A.9)

This justifies the exactness of the mean-field approximation for self-gravitating systems. Note that $\eta/N$ can be identified as the gravitational version of the "plasma parameter''. This is similar to the remark of Lundgren & Pointin (1977) for the point vortex gas. Now, plugging this result in Eq. (A.4), we find that, for $N\rightarrow +\infty$,

$$% {\partial P_{1}\over\partial {\vec r}_{1}}({\vec r}_{1})=-\... ...tial u_{12}\over\partial {\vec r}_{1}}{\rm d}^{3}{\vec r}_{2}.$$ (A.10)

Integrating with respect to ${\vec r}_1$ and introducing the mean density $\rho({\vec r})=\langle \sum_{i}m\delta({\vec r}_{i}-{\vec r})\rangle=NmP_{1}({\vec r})$, we obtain the Boltzmann distribution

$$% \rho=A{\rm e}^{-\beta m\Phi},$$ (A.11)

where $\Phi({\vec r})=\int \rho({\vec r}')u({\vec r}-{\vec r}'){\rm d}^{3}{\vec r}'$ is the self-consistent gravitational potential. Adding the gaussian velocity factor, we obtain the Maxwell-Boltzmann distribution

$$% f=A'{\rm e}^{-\beta m({v^2\over 2}+\Phi)}.$$ (A.12)

As we have seen in Sect. 2, the distribution function (A.12) can also be obtained by minimizing the Boltzmann free energy  $F_{\rm B}[f]$ at fixed mass M and temperature T. This method provides a condition of thermodynamical stability $\delta^{2}F\ge 0$, which is not captured by the equilibrium BBGKY hierarchy; it just yields critical points of free energy. To get the condition of stability, we need to consider time-dependent solutions and kinetic equations, i.e. the non-equilibrium BBGKY hierarchy. Indeed, the thermodynamical stability is related to the dynamical stability with respect to the Fokker-Planck equation (Chavanis 2004c).

The equation for the two-body distribution function (A.8) is complicated because the one-body distribution function is spatially inhomogeneous. It may be of interest, however, to advocate the Jeans swindle and consider, formally, the case of an infinite homogeneous self-gravitating system (this can be made rigorous in a cosmological context; see Kandrup 1983). Making the drastic approximation $P_1=\rho/M$where $\rho$ is a constant, Eq. (A.8) simplifies into

$$% {\partial h\over\partial {\vec r}_{1}}({\vec r}_{1},{\vec r... ...tial u_{13}\over\partial {\vec r}_{1}}{\rm d}^{3}{\vec r}_{3},$$ (A.13)

where the second term in the right hand side of Eq. (A.8), of order 1/N², has been neglected. The correlation function h is defined by

$$% P_{2}(\vert{\vec x}\vert)={\rho^{2}\over M^{2}}\lbrack 1+h(\vert{\vec x}\vert)\rbrack,$$ (A.14)

where ${\vec x}={\vec r}_{1}-{\vec r}_{2}$. Taking the divergence of Eq. (A.13) and using $\Delta u=4\pi G m^{2}\delta({\vec r}_1-{\vec r}_2)$, we obtain

$$% \Delta h+k_{J}^{2} h=-4\pi G\beta m^{2}\delta({\vec x}),$$ (A.15)

where $k_{J}=(4\pi G m\beta\rho )^{1/2}$ is the inverse of the Jeans length. This equation is easily integrated to yield

$$% h(x)=\beta G m^{2}{\cos(k_{J}x)\over x}\cdot$$ (A.16)

This is the counterpart of the Debye-Hückel result in the gravitational case (Kandrup 1983). We emphasize that the above results are valid for other systems with long-range interactions (Chavanis 2004b). In particular, for the HMF model for which a homogeneous phase rigorously exists, we find by the same method that $h(\theta)={2\over N}{{\beta/\beta_{\rm c}}\over 1-\beta/\beta_{\rm c}}\cos\theta$, where $\beta_{\rm c}={4\pi\over kM}$ is the critical inverse temperature. In particular, the correlation function diverges close to the critical point where the homogeneous phase becomes unstable, so that the mean-field approximation ceases to be valid. We expect a similar behavior for inhomogenous self-gravitating systems close to $T_{\rm c}$.

Considering now an isolated Hamiltonian system, the N-body microcanonical distribution function is given by

 

$$% P_{N}({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N})=... ...ver g(E)}\delta (E-H({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N})).$$ (A.17)

From this expression it is easy to write the equilibrium BBGKY hierarchy (Chavanis 2004b). The first two equations of this hierarchy are

 

$$% {\partial P_{1}\over\partial {\vec r}_{1}}(1)\!=\!-(N-1)\int {\... ...l \over\partial E}\biggl\lbrack\! g(E)P_{2}(1,2)\!\biggr\rbrack {\rm d}^{3}(2),$$ (A.18)

 

$${\partial P_{2}\over\partial {\vec r}_{1}}({1},{2})=- {\partial u... ...ver\partial E}\biggl\lbrack g(E)P_{3}(1,2,3)\biggr\rbrack {\rm d}^{3}(3),\qquad$$ (A.19)

where we have written $(j)=({\vec r}_{j},{\vec v}_{j})$. Now,

$$% {1\over g(E)}{\partial \over\partial E}\biggl\lbrack g(E)P_{j}\biggr\rbrack=\beta P_{j}+{\partial P_{j}\over\partial E},$$ (A.20)

where $\beta=\partial S/\partial E$ and $S(E)=\ln g(E)$. The ratio of $\partial P_{j}/\partial E$ on $\beta P_{j}$ is of order ${1\over E\beta}={1\over \Lambda\eta N}$. Therefore, in the thermodynamic limit $N\rightarrow +\infty$ with $\Lambda$, $\eta$ fixed, the second term in the r.h.s. of Eq. (A.20) is always negligible with respect to the first. To leading order in N, we obtain the same equations as in the canonical ensemble. Therefore, the mean-field approximation is exact and leads to the Boltzmann distribution (A.11). Observing that

 

$$% {\partial P_{1}\over\partial {\vec v}_{1}}(1)=- {1\over g(E)}{\partial \over\partial E}\bigl \lbrack g(E)P_{1}(1)\bigr \rbrack {\vec v}_{1},$$ (A.21)

and taking the $N\rightarrow +\infty$ limit, we find that $P_{1}\sim {\rm e}~^{-\beta v_{1}^{2}/2}$. Combined with Eq. (A.11), this leads to the Maxwell-Boltzmann distribution (A.12). Therefore, the equilibrium BBGKY hierarchy in the microcanonical ensemble leads to the same result (A.12) as in the canonical ensemble. As indicated previously, the inequivalence of ensembles will appear by considering the non-equilibrium BBGKY hierarchy. The thermodynamical stability in the microcanonical ensemble is connected to the dynamical stability with respect to the Landau equation (see Chavanis 2004c) which can be deduced from the non-equilibrium BBGKY hierarchy (Balescu 1963) to order 1/N.

   
Appendix B: Derivation of Eq. (128)

The phase space hypervolume with energy less than E is defined by

 

$$% I(E)=\int_{H\le E}\prod_{i=1}^{N}{\rm d}^{3}{\vec r}_{i}{\rm d}^{3}{\vec v}_{i}.$$ (B.1)

Integrating over the velocities and using the fact that the kinetic term in the Hamiltonian is quadratic, a standard calculation yields

 

$$% I(E)=A\int \bigl \lbrack E-U({\vec r}_{1},...,{\vec r}_{N})\bigr \rbrack^{3N\over 2}\prod_{i=1}^{N}{\rm d}^{3}{\vec r}_{i},$$ (B.2)

where A=(2/m)^3N/2V_3N and V_n is the volume of a unit-hypersphere in a space of dimension n. The density of states $g(E)={\rm d}I/{\rm d}E$ is therefore given by

 

$$% g(E)={3N\over 2}A\int \bigl \lbrack E-U({\vec r}_{1},...,{\vec r}_{N})\bigr \rbrack^{{3N\over 2}-1}\prod_{i=1}^{N}{\rm d}^{3}{\vec r}_{i}.$$ (B.3)

Assuming now that $P_{N}({\vec r}_{1},{\vec v}_{1},...,{\vec r}_{N},{\vec v}_{N},t)\simeq P_{N}(E,t)$, and substituting this ansatz in the N-body Fokker-Planck equation Eq. (126), we obtain after simplification

 

$$% {\partial P_{N}\over\partial t}=2m\bigl\lbrack E-U({\vec r}_{1}... ...i P_{N}\biggr ) +3Nm\biggl (D{\partial P_{N}\over\partial E}+\xi P_{N}\biggr ),$$ (B.4)

where the term in bracket is $\sum_{i=1}^{N}v_{i}^{2}$. We note that P_N=P_N(E,t) is not an exact solution of (126), as expected. To get rid of the dependence in ${\vec r}_{1},...,{\vec r}_{N}$, we shall average Eq. (B.4) over the hypersurface of energy E using

 

$$% \overline{X}(E)={\int (E-U)^{{3N\over 2}-1}X({\vec r}_{1},...,{... ... r}_{i}\over \int (E-U)^{{3N\over 2}-1}\prod_{i=1}^{N}{\rm d}^{3}{\vec r}_{i}},$$ (B.5)

according to Eq. (B.3). This gives

 

$$% g(E){\partial P_{N}\over\partial t}=3M I(E){\partial\over\parti... ...{N}\biggr ) +3M g(E)\biggl (D{\partial P_{N}\over\partial E}+\xi P_{N}\biggr ).$$ (B.6)

Using $g(E)={\rm d}I/{\rm d}E(E)$, we can put this equation in the form (128).

   
Appendix C: Rotating self-gravitating systems

In this Appendix, we briefly consider the case of rotating self-gravitating systems. Introducing dimensionless variables as in Sect. 2.3, the conservation of angular momentum is equivalent to ${\setbox 0=\hbox{$\lambda$ }\kern-.025em\copy0\kern-\wd0 \kern-0.05em\copy0\ker... ...copy0\kern-\wd0 \kern-0.05em\copy0\kern-\wd0\kern-.025em\raise.0233em\box0}[f']$ with

$$% {\setbox 0=\hbox{$\lambda$ }\kern-.025em\copy0\kern-\wd0 \k... ...r}'{\times}{\vec v}' {\rm d}^{3}{\vec r}'{\rm d}^{3}{\vec v}'.$$ (C.1)

Now, repeating the argumentation of Sect. 3, the density of states

 

$$% g(E,{\vec L})=\int \delta (E-H)\delta \biggl ({\vec L}-\sum_{i=... ...ec v}_{i}\biggr )\prod_{i=1}^{N}{\rm d}^{3}{\vec r}_{i}{\rm d}^{3}{\vec v}_{i},$$ (C.2)

can be written formally as

 

$$% g(E,{\vec L})=\int {\cal D}f \ {\rm e}^{S[f]} \delta(E-E[f]) \delta(M-M[f])\delta({\vec L}-{\vec L}[f]).$$ (C.3)

Similarly, the partition function

$$% Z(\beta,{\bf\Omega})=\int {\rm e}^{-\beta H+\beta{\bf\Omega... ...}\prod_{i=1}^{N}{\rm d}^{3}{\vec r}_{i}{\rm d}^{3}{\vec v}_{i}$$ (C.4)

can be written as

$$% Z(\beta,{\bf\Omega})=\int {\rm e}^{-\beta E+\beta{\bf\Omega}\cdot{\vec L}}g(E,{\vec L}){\rm d}E{\rm d}^{3}{\vec L},$$ (C.5)

or

$$% Z(\beta,{\bf\Omega})=\int {\cal D}f \ {\rm e}^{J[f]}\delta(M-M[f]),$$ (C.6)

where $J[f]=S[f]-\beta E[f]+\beta {\setbox 0=\hbox{$\Omega$ }\kern-.025em\copy0\kern-\wd0 \kern-0.05em\copy0\kern-\wd0\kern-.025em\raise.0233em\box0}\cdot {\vec L}[f]$ is the free energy. In order to apply the saddle point approximation, we just need to impose that $\Lambda=-ER/GM^{2}$, $\eta=\beta GMm/R$, ${\setbox 0=\hbox{$\lambda$ }\kern-.025em\copy0\kern-\wd0 \kern-0.05em\copy0\kern-\wd0\kern-.025em\raise.0233em\box0}={\vec L}/(GM^{3}R)^{1/2}$ and ${\bf\omega}={\bf\Omega}(R^{3}/GM)^{1/2}$ remain of order unity in the limit $N\rightarrow +\infty$ (in the case of self-gravitating fermions, we also need to impose that $\mu=\eta_0\sqrt{G^3 M R^3}$ is fixed and in the case of a soften potential that $\epsilon=r_{0}/R$ is fixed). This defines the thermodynamic limit for rotating self-gravitating systems. The corresponding scalings are given in Chavanis & Rieutord (2003). In particular, $S\sim N$ and $J\sim N$. Therefore, in the $N\rightarrow \infty$ limit, we have to maximize S[f] at fixed E, M and ${\vec L}$ in the microcanonical ensemble and we have to maximize $J[f]=S[f]-\beta E[f]+\beta\ {\bf\Omega}\cdot {\vec L}[f]$ at fixed $\beta$, M and  ${\bf\Omega}$ in the canonical ensemble. Computation of rotating self-gravitating systems in relation with statistical mechanics have been performed by Votyakov et al. (2002) for a classical gas on a lattice with a large spacing and by Chavanis & Rieutord (2003) for fermions with arbitrary degree of degeneracy.

References

• Antonov, V. A. 1962, Vest. Leningr. Gos. Univ., 7, 135 (In the text)
• Aronson, E. B., & Hansen, C. J. 1972, ApJ, 177, 145 [NASA ADS] [CrossRef] (In the text)
• Balescu, R. 1963, Statistical Mechanics of Charged Particles (New-York: Interscience) (In the text)
• Bouchet, F., & Barré, J. 2005, J. Stat. Phys., submitted [cond-mat/0303307] (In the text)
• Chavanis, P. H. 2002a, A&A, 381, 340 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Chavanis, P. H. 2002b, PRE, 65, 056123 [NASA ADS] [MathSciNet] (In the text)
• Chavanis, P. H. 2002c, Statistical mechanics of two-dimensional vortices and stellar systems, in Dynamics and thermodynamics of systems with long range interactions, ed. T. Dauxois, S. Ruffo, E. Arimondo, & M. Wilkens, Lect. Notes Phys., 602 (Springer) [cond-mat/0212223] (In the text)
• Chavanis, P. H. 2002d, The self-gravitating Fermi gas, in Proceedings of the Fourth International Heidelberg Conference on Dark Matter in Astro and Particle Physics, ed. Klapdor-Kleingrothaus (New-York: Springer) [arXiv:astro-ph/0205426] (In the text)
• Chavanis, P. H. 2003, A&A, 401, 15 [EDP Sciences] [NASA ADS] (In the text)
• Chavanis, P. H. 2004a, PRE, 69, 066126 [NASA ADS] (In the text)
• Chavanis, P. H. 2004b [arXiv:cond-mat/0409641] (In the text)
• Chavanis, P. H. 2004c, Physica A, 332, 89 [NASA ADS] [MathSciNet] (In the text)
• Chavanis, P. H., & Ispolatov, I. 2002, PRE, 66, 036109 [NASA ADS] (In the text)
• Chavanis, P. H., & Rieutord, M. 2003, A&A, 412, 1 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Chavanis, P. H., Ribot, M., Rosier, C., & Sire, C. 2004 (Banach Center Publications), in press [cond-mat/0407386] (In the text)
• Chavanis, P. H., Rosier, C., & Sire, C. 2002, PRE, 66, 036105 [NASA ADS] (In the text)
• Chavanis, P. H., & Sommeria, J. 1998, MNRAS, 296, 569 [NASA ADS] [CrossRef] (In the text)
• de Vega, H. J., & Sanchez, N. 2002, Nucl. Phys. B, 625, 409 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Gross, D. H. E. 2003 [cond-mat/0307535] (In the text)
• Gross, D. H. E. 2004 [cond-mat/0403582] (In the text)
• Hertel, P., & Thirring, W. 1971, in Quanten und Felder, ed. H. P. Dürr (Vieweg: Braunschweig) (In the text)
• Horwitz, G., & Katz, J. 1978, ApJ, 222, 941 [NASA ADS] [CrossRef] (In the text)
• Huang, K. 1963, Statistical mechanics (John Wiley) (In the text)
• Ispolatov, I., & Karttunen, M. 2003, PRE, 68, 036117 [NASA ADS] (In the text)
• Kandrup, H. 1983, Astrophy. Space. Sci., 89, 143 [NASA ADS] [MathSciNet] (In the text)
• Katz, J. 1978, MNRAS, 183, 765 [NASA ADS] (In the text)
• Katz, J. 2003, Found. Phys., 33, 223 [CrossRef] [MathSciNet] (In the text)
• Katz, J., & Okamoto, I. 2000, MNRAS, 317, 163 [NASA ADS] [CrossRef] (In the text)
• Kiessling, M. K. H. 1989, J. Stat. Phys., 55, 203 [CrossRef] [MathSciNet] (In the text)
• Kramers, H. A. 1940, Physica, 7, 284 [CrossRef] [MathSciNet] (In the text)
• Lundgren, T. S., & Pointin, Y. B. 1977, J. Stat. Phys., 17, 323 [CrossRef] (In the text)
• Lynden-Bell, D. 1999, Physica A, 263, 293 [NASA ADS] (In the text)
• Lynden-Bell, D., & Lynden-Bell, R. M. 1977, MNRAS, 181, 405 [NASA ADS] (In the text)
• Lynden-Bell, D., & Wood, R. 1968, MNRAS, 138, 495 [NASA ADS] (In the text)
• Monaghan, J. J. 1978, MNRAS, 184, 25 [NASA ADS] (In the text)
• Ogorodnikov, K. F. 1965, Dynamics of stellar systems (Pergamon) (In the text)
• Padmanabhan, T. 1990, Phys. Rep., 188, 285 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Risken, H. 1989, The Fokker-Planck equation (Springer) (In the text)
• Sire, C., & Chavanis, P. H. 2002, PRE, 66, 046133 [NASA ADS] (In the text)
• Stahl, B., Kiessling, M. K. H., & Schindler, K. 1994, Planet. Space Sci., 43, 271 [NASA ADS] (In the text)
• Thirring, W. 1970, Z. Phys., 235, 339 (In the text)
• Votyakov, E. V., De Martino, A., & Gross, D. H. E. 2002, Eur. Phys. J. B, 29, 593 [EDP Sciences] [NASA ADS] (In the text)
• Youngkins, V., & Miller, B. 2000, PRE, 62, 4583 [NASA ADS] (In the text)