text:
# STATISTICAL MECHANICS OF THE SELF-GRAVITATING GAS: I. THERMODYNAMIC LIMIT AND PHASE DIAGRAMS

## 1 Statistical Mechanics of the Self-Gravitating Gas

Physical systems at thermal equilibrium are usually homogeneous. This is the case for gases with short range intermolecular forces (and in absence of external fields). In such cases the entropy is maximum when the system homogenizes.

When long range interactions as the gravitational force are present, even the ground state is inhomogeneous. In this case, each element of the substance is acted on by very strong forces due to distant particles of the gas. Hence, regions near to and far from the boundary of the volume occupied by the gas will be in very different conditions, and, as a result, the homogeneity of the gas is destroyed . The state of maximal entropy for gravitational systems is inhomogeneous. This basic inhomogeneity suggested us that fractal structures can arise in a self-interacting gravitational gas.

The inhomogeneous character of the ground state for gravitational systems explains why the universe is not going towards a ‘thermal death’. A ‘thermal death’ would mean that the universe evolves towards more and more homogeneity. This can only happen if the entropy is maximal for an homogeneous state. Instead, it is the opposite what happens, structures are formed in the universe through the action of the gravitational forces as time evolves.

Usual theorems in statistical mechanics break down for inhomogeneous ground states. For example, the specific heat may be negative in the microcanonical ensemble (not in the canonical ensemble where it is always positive).

As is known, the thermodynamic limit for self-gravitating systems does not exist in its usual form ($`N\mathrm{},V\mathrm{},N/V=`$ fixed). The system collapses into a very dense phase which is determined by the short distance (non-gravitational) forces between the particles.

We instead find that the thermodynamic functions exist in the dilute limit

$$N\mathrm{},V\mathrm{},\frac{N}{V^{1/3}}=\text{fixed}$$
(1)
where $`V`$ stands for the volume of the box containing the gas. In such a limit, the energy $`E`$, the free energy and the entropy turns to be extensive. That is, we find that they take the form of $`N`$ times a function of

$$\eta =\frac{Gm^2N}{LT}\text{or}\xi =\frac{EL}{Gm^2N^2}$$

where $`\eta `$ and $`\xi `$ are intensive variables. Namely, $`\eta `$ and $`\xi `$ stay finite when $`N`$ and $`VL^3`$ tend to infinite. $`\eta `$ is appropriate for the canonical ensemble and $`\xi `$ for the microcanonical ensemble. Physical magnitudes as the specific heat, speed of sound, chemical potential and compressibility only depend on $`\eta `$ or $`\xi `$. $`\eta `$ and $`\xi `$ as well as the ratio $`N/L`$ are therefore intensive magnitudes. The energy, the free energy, the Gibbs free energy and the entropy are of the form $`N`$ times a function of $`\eta `$. These functions of $`\eta `$ have a finite $`N=\mathrm{}`$ limit for fixed $`\eta `$ (once the ideal gas contributions are subtracted). Moreover, the dependence on $`\eta `$ in all these magnitudes express through a single universal function $`f(\eta )`$.

We study here and in the companion paper (called paper II in what follows) the statistical mechanics of the self-gravitating gas. That is, our starting point is the partition function for non-relativistic particles interacting through their gravitational attraction in thermal equilibrium. We study the self-gravitating gas in the three ensembles: microcanonical (MCE), canonical (CE) and grand canonical (GCE). We performed calculations by three methods:

* By expanding the partition function through direct calculation in powers of $`1/\xi `$ and $`\eta `$ for the MCE and CE, respectively. These expressions apply in the dilute regime ($`\xi 1,\eta 1`$) and become identical for both ensembles for $`N\mathrm{}`$. At $`\eta =0=1/\xi `$ we recover the ideal gas behaviour.
* By performing Monte Carlo simulations both in the MCE and in the CE. We found in this way that the self-gravitating gas collapses at a critical point which depends on the ensemble considered. As shown in fig. 1 the collapse occurs first in the canonical ensemble (point T). The microcanonical ensemble exhibits a larger region of stability that ends at the point MC (fig. 1). Notice that the physical magnitudes are identical in the common region of validity of both ensembles within the statistical error. Beyond the critical point T the system becomes suddenly extremely compact with a large negative pressure in the CE. Beyond the point MC in the MCE the pressure and the temperature increase suddenly and the gas collapses. The phase transitions at T and at MC are of zeroth order since the Gibbs free energy has discontinuities in both cases.
* By using the mean field approach we evaluate the partition function for large $`N`$. We do this computation in the grand canonical, canonical and microcanonical ensembles. In the three cases the partition function is expressed as a functional integral over a statistical weight which depends on the (continuous) particle density. These statistical weights are of the form of the exponential of an ‘effective action’ proportional to $`N`$. Therefore, the $`N\mathrm{}`$ limit follows by the saddle point method. The saddle point is a space dependent mean field showing the inhomogeneous character of the ground state. Corrections to the mean field are of the order $`1/N`$ and can be safely ignored for $`N1`$ except near the critical points. These mean field results turned out to be in excellent agreement with the Monte Carlo results and with the low density expansion.

We calculate the saddle point (mean field) for spherical symmetry and we obtain from it the various physical magnitudes (pressure, energy, entropy, free energy, specific heats, compressibilities, speed of sound and particle density). Furthermore, we compute the determinants of small fluctuations around the saddle point solution for spherical symmetry for the three statistical ensembles in paper II.

When any small fluctuation around the saddle point decreases the statistical weight in the functional integral, the saddle point is dominating the integral and the mean field approach is fully valid. In that case the determinant of small fluctuations is positive. A negative determinant of small fluctuations indicates that some fluctuations around the saddle point are increasing the statistical weight in the functional integral and hence the saddle point does not dominate the partition function. The mean field approach cannot be used when the determinant of small fluctuations is negative.

The zeroes of the small fluctuations determinant determine the position of the critical points for the three statistical ensembles. The Monte Carlo simulations for the CE and the MCE show that the self-gravitating gas collapses near the critical points obtained from mean field.

The saddle point solution is identical for the three statistical ensembles. This is not the case for the fluctuations around it. The presence of constraints in the CE (on the number of particles) and in the MCE (on the energy and the number of particles) changes the functional integral over the quadratic fluctuations with respect to the GCE.

The saddle point of the partition function turns out to coincide with the hydrostatic treatment of the self-gravitating gas-. (Usually known as the ‘isothermal sphere’ in the spherically symmetric case).

Our Monte Carlo simulations are performed in a cubic geometry. The equilibrium configurations obtained in this manner can thus be called the ‘isothermal cube’.

We find for spherical symmetry: $`\eta _{GC}^R=0.797375\mathrm{},\eta _C^R=2.517551\mathrm{}`$ and $`\eta _{MC}^R=2.03085\mathrm{}`$. The variable $`\eta ^R`$ appropriate for a spherical symmetry is defined as $`\eta ^R\frac{Gm^2N}{RT}=\eta \left(\frac{4\pi }{3}\right)^{1/3}=1.61199\mathrm{}\eta `$.

The values of $`\eta _T`$ and $`\eta _C`$ change by a few percent with the geometry and with the number of particles (for large $`N>500`$). For spherical symmetry and $`N=\mathrm{}`$ (mean field) we obtain $`\eta _C=\left(\frac{3}{4\pi }\right)^{1/3}\eta _C^R=1.56176\mathrm{}`$. Our Monte Carlo simulations yield $`\eta _T1.515`$. We find from the mean field approach that the isothermal compressibility diverges at $`\eta =\eta _0=1.51024\mathrm{}\eta _T`$ for spherical symmetry.

The conclusion being that the MF correctly describes the thermodynamic limit except near the critical points (where the small fluctuations determinant vanishes); the MF is valid for $`N|\eta \eta _{crit}|1`$. The vicinity of the critical point should be studied in a double scaling limit $`N\mathrm{},\eta \eta _{crit}`$.

In summary, the picture we get from our calculations using these three methods show that the self-gravitating gas behaves as a perfect gas for $`\eta 0,\mathrm{\hspace{0.33em}1}/\xi 0`$. When $`\eta `$ and $`1/\xi `$ grow, the gas becomes denser till it suddenly condenses into a high density object at a critical point GC, C or MC depending upon the statistical ensemble chosen. In the Monte Carlo simulations for the CE the collapse takes place at the point T slightly before $`\eta _C`$. $`\eta `$ is related with the Jeans’ length $`d_J`$ of the gas through $`\eta =3(L/d_J)^2`$. Hence, when $`\eta `$ goes beyond $`\eta _T`$, the length of the system becomes larger than $`d_J/\sqrt{\eta _T/3}`$. The collapse at T in the CE is therefore a manifestation of the Jeans’ instability. The saddle point ceases to describe the physics at C since the determinant of fluctuations for the CE vanishes there.

In the MCE, the determinant of fluctuations vanishes at the point MC. The physical states beyond MC are collapsed configurations shown by the Monte Carlo simulations \[see fig. 4\]. Actually, the gas collapses in the Monte Carlo simulations slightly before the mean field prediction for the point MC. The phase transition at the critical point MC is the so called gravothermal catastrophe .

The gravitational interaction being attractive without lower bound, a short distance cut-off ($`A`$) must be introduced in order to give a meaning to the partition function. We take the gravitational force between particles as $`Gm^2/r^2`$ for $`r>A`$ and zero for $`r<A`$ where $`r`$ is the distance between the two particles. We show that the cut-off effects are negligible in the $`N=\mathrm{}`$ limit. That is, once we set $`N=\mathrm{}`$ with fixed $`\eta `$, all physical quantities are finite in the zero cut-off limit ($`A=0`$). The cut-off effects are of the order $`A^2/L^2`$ and can be safely ignored.

All physical quantities are expressed in terms of $`f(\eta )`$. Besides computing $`f(\eta )`$ numerically in the mean field approach, we obtain analytic results about it from the Abel’s equation. There is a square root branch point in $`f(\eta )`$ at $`\eta _C`$. The specific heat is positive in the first sheet and negative in the second sheet. This second sheet is only physically realized in the microcanonical ensemble (MCE). \[The specific heat is positive definite in the CE\]. $`f(\eta )`$ has infinitely many branches in the $`\eta `$ plane but only the first two are physically realized. Beyond MC the states described by the mean field saddle point are unstable. We plot and analyze the equation of state, the energy, the entropy, the free energy, $`c_V`$, $`c_P`$, the isothermal compressibility and the speed of sound \[figs. 10-15\]. Most of these physical magnitudes were not previously computed in the literature as functions of $`\eta `$.

We find analytically the behaviour of $`f(\eta )`$ near $`\eta _C`$ in mean field,

$$f_{MF}(\eta ^R)\stackrel{\eta ^R\eta _C^R}{=}\frac{1}{3}+0.213738\mathrm{}\sqrt{\eta _C^R\eta ^R}+0.172225\mathrm{}(\eta _C^R\eta ^R)+𝒪\left[(\eta _C^R\eta ^R)^{3/2}\right].$$

This shows that the specific heat at constant volume diverges as $`\left(\eta _C\eta \right)^{1/2}`$ for $`\eta ^R\eta _C^R`$. The specific heat at constant pressure and the isothermal compressibility diverge at $`\eta _0`$ as $`\left(\eta _0\eta \right)^1`$. These mean field results apply for $`|\eta \eta _C|1N|\eta \eta _C|`$. Fluctuations around mean field can be neglected in such a regime.

The Monte Carlo calculations permit us to obtain $`f(\eta )`$ in the collapsed phase. Such result (which is cutoff dependent) cannot be obtained in the mean field approach. The mean field only provides information (as $`f(\eta )`$) in the gas phase.

For the self-gravitating gas, we find that the Gibbs free energy $`\mathrm{\Phi }`$ is not equal to $`N`$ times the chemical potential and that the thermodynamic potential $`\mathrm{\Omega }`$ is not equal to $`PV`$ as usual. This is a consequence of the dilute thermodynamic limit $`N\mathrm{},L\mathrm{},N/L=`$fixed.

We compute local properties of the gas in paper II. That is, the local energy density $`ϵ(r)`$, local particle density and local pressure. Furthermore, we analyze the scaling behaviour of the particle distribution and its fractal (Haussdorf) dimension.

This paper is organized as follows. In section II we present the statistical mechanics of the self-gravitating gas in the microcanonical ensemble, in sec. III we do the analogous presentation for the canonical ensemble and in sec. IV we contrast the results for the CE and the MCE. Sec. V contains the results from Monte Carlo simulations and we develop in sec. VI the mean field approach. In sec. VII we present the results for intensive magnitudes. Discussion and remarks are presented in section VIII whereas appendixes A-C contain relevant mathematical developments.

## 2 Statistical Mechanics of the Self-Gravitating Gas: the microcanonical ensemble

We investigate in this section an isolated set of $`N`$ non-relativistic particles with mass $`m`$ interacting through Newtonian gravity with total energy $`E`$. That is, a self-gravitating gas in the microcanonical ensemble. We assume the system being on a cubic box of side $`L`$ just for simplicity. We consider spherical symmetry in sec. VI. Please notice that we never use periodic boundary conditions.

At short distances, the particle interaction for the self-gravitating gas in physical situations is not gravitational. Its exact nature depends on the problem under consideration (opacity limit, Van der Waals forces for molecules etc.). We shall just assume a repulsive short distance potential, that is,

$$v_A(|\stackrel{}{q}_l\stackrel{}{q}_j|)=\frac{1}{|\stackrel{}{q}_l\stackrel{}{q}_j|_A}=\{\begin{array}{c}\frac{1}{|\stackrel{}{q}_l\stackrel{}{q}_j|}\text{for}|\stackrel{}{q}_l\stackrel{}{q}_j|A\\ \\ +\frac{1}{A}\text{for}|\stackrel{}{q}_l\stackrel{}{q}_j|A\end{array}$$
(2)
where $`A<<L`$ is the short distance cut-off.

The presence of the repulsive short-distance interaction prevents the collapse (here unphysical) of the self-gravitating gas. In the situations we are interested to describe (interstellar medium, galaxy distributions) the collapse situation is unphysical.

The entropy of the system can be written as

$$S(E,N)=\mathrm{log}\left\{\frac{1}{N!}\mathrm{}\underset{l=1}{\overset{N}{}}\frac{d^3p_ld^3q_l}{(2\pi )^3}\delta \left[E\underset{l=1}{\overset{N}{}}\frac{p_l^2}{2m}U(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]\right\}$$
(3)
where

$$U(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)=Gm^2\underset{1l<jN}{}\frac{1}{|\stackrel{}{q}_l\stackrel{}{q}_j|_A}$$
(4)
and $`G`$ is Newton’s gravitational constant.

In order to compute the integrals over the momenta $`p_l,(1lN)`$, we introduce the variables,

$$\stackrel{}{\rho }_i=\frac{1}{\sqrt{2m}}\stackrel{}{p}_i.$$

We can now integrate over the angles in $`3N`$ dimensions,

$`{\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \underset{l=1}{\overset{N}{}}}{\displaystyle \frac{d^3p_l}{(2\pi )^3}}\delta \left[E{\displaystyle \underset{l=1}{\overset{N}{}}}\stackrel{}{\rho }_l^{\mathrm{\hspace{0.17em}2}}U(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]`$ (5)
$`=`$ $`\left({\displaystyle \frac{\sqrt{2m}}{2\pi }}\right)^{3N}{\displaystyle \frac{2\pi ^{3N/2}}{\mathrm{\Gamma }\left(\frac{3N}{2}\right)}}{\displaystyle _0^{\mathrm{}}}\rho ^{3N1}𝑑\rho \delta \left[E\rho ^2U(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]`$ (7)
$`=`$ $`\left({\displaystyle \frac{m}{2\pi }}\right)^{3N/2}{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\frac{3N}{2}\right)}}\left[EU(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]^{3N/21}\theta \left[EU(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]`$ (9)
The delta function in the energy thus becomes the constraint of a positive kinetic energy $`EU(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)>0`$. We then get for the entropy,

$$S(E,N)=\mathrm{log}\left\{\frac{\left(\frac{m}{2\pi }\right)^{3N/2}}{N!\mathrm{\Gamma }\left(\frac{3N}{2}\right)}_0^L\mathrm{}_0^L\underset{l=1}{\overset{N}{}}d^3q_l\left[EU(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]^{3N/21}\theta \left[EU(\stackrel{}{q}_1,\mathrm{}\stackrel{}{q}_N)\right]\right\}$$
(10)
It is convenient to introduce the dimensionless variables $`\stackrel{}{r}_l,\mathrm{\hspace{0.33em}1}lN`$ making explicit the volume dependence as

$`\stackrel{}{q}_l`$ $`=`$ $`L\stackrel{}{r}_l,\stackrel{}{r}_l=(x_l,y_l,z_l),`$ (11)
$`0`$ $``$ $`x_l,y_l,z_l1.`$ (13)
That is, in the new coordinates the gas is inside a cube of unit volume.

The entropy then becomes

$`S(E,N)`$ $`=`$ $`\mathrm{log}\{{\displaystyle \frac{N^{3N2}m^{9N/22}L^{3N/2+1}G^{3N/21}}{N!\mathrm{\Gamma }\left(\frac{3N}{2}\right)(2\pi )^{3N/2}}}`$
$`{\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}{\displaystyle \underset{l=1}{\overset{N}{}}}d^3r_l[\xi +{\displaystyle \frac{1}{N}}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)]^{3N/21}\theta [\xi +{\displaystyle \frac{1}{N}}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)]\}`$
where we introduced the dimensionless variable $`\xi `$,

$$\xi \frac{EL}{Gm^2N^2}$$
(16)
and

$$u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\frac{1}{N}\underset{1l<jN}{}\frac{1}{|\stackrel{}{r}_l\stackrel{}{r}_j|_a}.$$
(17)
where $`aA/L1`$.

Let us define the coordinate partition function in the microcanonical ensemble as

$$w(\xi ,N)_0^1\mathrm{}_0^1\underset{l=1}{\overset{N}{}}d^3r_l\left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]^{3N/21}\theta \left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right].$$
(18)
Therefore,

$$S(E,N)=\mathrm{log}\left[\frac{N^{3N2}m^{9N/22}L^{3N/2+1}G^{3N/21}}{N!\mathrm{\Gamma }\left(\frac{3N}{2}\right)(2\pi )^{3N/2}}\right]+\mathrm{log}w(\xi ,N).$$

We can now compute the thermodynamic quantities, temperature and pressure through the standard thermodynamic relations

$$\frac{1}{T}=\left(\frac{S}{E}\right)_V\text{and}p=T\left(\frac{S}{V}\right)_E,$$
(19)
where $`VL^3`$ stands for the volume of the system and $`p`$ is the external pressure on the system.

We obtain the temperature as a function of $`E`$ and $`\xi `$ from eqs.(2) and (19)

$$\frac{1}{T}=\frac{\xi }{E}\frac{}{\xi }\mathrm{log}w(\xi ,N)=\frac{3N\xi }{2E}\left[1\frac{2}{3N}\right]<\frac{1}{\xi +\frac{1}{N}u(.)}>.$$
(20)
where

$$<\frac{1}{\xi +\frac{1}{N}u(.)}>\frac{_0^1\mathrm{}_0^1\underset{l=1}{\overset{N}{}}d^3r_l\left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]^{3N/22}\theta \left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]}{_0^1\mathrm{}_0^1_{l=1}^Nd^3r_l\left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]^{3N/21}\theta \left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]}$$
(21)
The equation of state follows from eqs.(2) and (19)

$$\frac{pV}{NT}=\frac{1}{2}+\frac{1}{3N}+\frac{\xi }{3N}\frac{}{\xi }\mathrm{log}w(\xi ,N)=\frac{1}{2}\left(1+\frac{2}{3N}\right)+\frac{\xi }{2}<\frac{1}{\xi +\frac{1}{N}u(.)}>\left[1\frac{2}{3N}\right].$$
(22)
We are interested in the large size limit where $`N\mathrm{},L\mathrm{}`$ and $`E\mathrm{}`$. We consider that $`\xi =\frac{EL}{Gm^2N^2}`$ stays fixed in such limit. That is, we assume $`E/N`$ and $`L/N`$ bounded and nonzero when $`E,L`$ and $`N\mathrm{}`$. We shall see below that such limit is meaningful.

It is possible to write the energy and the equation of state in terms of a single function

$$g(\xi )\frac{\xi }{N}\frac{}{\xi }\mathrm{log}w(\xi ,N)=\frac{3\xi }{2}<\frac{1}{\xi +\frac{1}{N}u(.)}>\left[1\frac{2}{3N}\right].$$
(23)
We find from eqs.(20), (22) and (23),

$`{\displaystyle \frac{pV}{NT}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{3}}g(\xi )+{\displaystyle \frac{1}{3N}}`$ (24)
$`{\displaystyle \frac{E}{NT}}`$ $`=`$ $`g(\xi ).`$ (26)
We obtain the virial theorem by eliminating $`g(\xi )`$ in the eqs.(24)

$$pV=\frac{NT}{2}+\frac{E+T}{3}.$$
(27)
where the term $`T/3`$ can be neglected for large $`N`$.

In the case of a perfect gas (no gravity) we have $`u(.)0,g(\xi )=\frac{3}{2},pV=NT`$ and $`E=\frac{3}{2}NT`$ as it must be.

The function $`g(\xi )`$ is computed by Monte Carlo simulations, mean field methods and, in the weak field limit $`\xi <<1`$, is calculated analytically in powers of $`1/\xi `$ in subsection II.A.

The specific heat per particle is given by

$$c_V=\frac{T}{N}\left(\frac{S}{T}\right)_V=\frac{1}{N\left(\frac{T}{E}\right)_V}$$

Hence, using eq.(24) yields

$$c_V=\frac{g(\xi )}{1\frac{\xi g^{}(\xi )}{g(\xi )}}\text{or}\frac{1}{c_V}=\frac{d}{d\xi }\left[\frac{\xi }{g(\xi )}\right]$$

We can relate the specific heat $`c_V`$ with the fluctuations as follows. We can express $`g(\xi )`$ as an average value using eqs.(21) and (23)

$$\frac{\xi }{g(\xi )}=\frac{2}{3\frac{2}{N}}\frac{_0^1\mathrm{}_0^1\underset{l=1}{\overset{N}{}}d^3r_l\left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]^{3N/21}\theta \left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]}{_0^1\mathrm{}_0^1_{l=1}^Nd^3r_l\left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]^{3N/22}\theta \left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]}$$

Computing the derivative with respect to $`\xi `$ yields,

$$\frac{1}{c_V}=\frac{2}{3}\frac{\left(\mathrm{\Delta }\frac{1}{\xi +\frac{1}{N}u(.)}\right)^2}{<\frac{1}{\xi +\frac{1}{N}u(.)}>^2}+𝒪\left(\frac{1}{N}\right)$$
(28)
where

$$\left(\mathrm{\Delta }\frac{1}{\xi +\frac{1}{N}u(.)}\right)^2N\left\{<\frac{1}{[\xi +\frac{1}{N}u(.)]^2}><\frac{1}{\xi +\frac{1}{N}u(.)}>^2\right\}$$

is of order $`N^0`$ for $`N\mathrm{}`$. \[Notice that in the calculation of the fluctuations we must keep the $`1/N`$ corrections till the end\].

We can express $`c_V`$ in terms of the fluctuations of the inverse temperature $`\beta 1/T`$ using eq.(20):

$$\frac{1}{c_V}=\frac{2}{3}\left(\frac{\mathrm{\Delta }\beta }{\beta }\right)^2$$
(29)
It must be noticed that in the microcanonical ensemble, $`c_V`$ may be positive as well as negative. In fact, it becomes negative when the fluctuations are large enough \[see sec. V and VI\].

We see that extensivity holds here in an specific way. $`T,S/N`$ and $`pV/N`$ are of order one for $`N\mathrm{}`$ provided $`\xi `$ stays fixed in such limit. That is, we must keep $`E/N`$ and $`L/N`$ fixed in the $`N\mathrm{}`$ limit.

### 2.1 The diluted regime: $`\xi >>1`$

We can obtain the thermodynamic quantities as a series in powers of $`1/\xi `$ just expanding the integrand in eq.(18).

We find

$`w(\xi ,N)`$ $`\stackrel{\xi \mathrm{}}{=}`$ $`\xi ^{3N/21}\{1+{\displaystyle \frac{9b_0N}{2\xi }}(1{\displaystyle \frac{2}{3N}})(1{\displaystyle \frac{1}{N}})`$ (30)
$`+`$ $`{\displaystyle \frac{9}{8\xi ^2}}(1{\displaystyle \frac{2}{3N}})(1{\displaystyle \frac{4}{3N}})[9N^2b_0^2(1{\displaystyle \frac{1}{N}})(1{\displaystyle \frac{2}{N}})(1{\displaystyle \frac{3}{N}})`$ (32)
$`+`$ $`Nb_1(1{\displaystyle \frac{1}{N}})(1{\displaystyle \frac{2}{N}})+{\displaystyle \frac{1}{2}}b_2(1{\displaystyle \frac{1}{N}})]+𝒪(\xi ^3)\}`$ (34)
where $`b_0,b_1`$ and $`b_2`$ are pure numbers,

$`b_0`$ $`=`$ $`{\displaystyle \frac{1}{6}}{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|}}`$ (35)
$`b_1`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3r_1d^3r_2d^3r_3}{|\stackrel{}{r}_1\stackrel{}{r}_2||\stackrel{}{r}_1\stackrel{}{r}_3|}}`$ (37)
$`b_2`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|^2}}`$ (39)
For the cubic geometry chosen, it takes the value

$$b_0^{cube}=\frac{4}{3}_0^1(1x)𝑑x_0^1(1y)𝑑y_0^1\frac{(1z)dz}{\sqrt{x^2+y^2+z^2}}=0.31372\mathrm{}.$$

For a sphere of unit volume we find

$`b_0^{sphere}`$ $`=`$ $`{\displaystyle \frac{1}{5}}\left({\displaystyle \frac{4\pi }{3}}\right)^{1/3}=0.32239839\mathrm{},`$ (40)
$`b_1^{sphere}`$ $`=`$ $`{\displaystyle \frac{51}{35}}\left({\displaystyle \frac{4\pi }{3}}\right)^{2/3}=3.786412026\mathrm{},`$ (42)
$`b_2^{sphere}`$ $`=`$ $`{\displaystyle \frac{9}{4}}\left({\displaystyle \frac{4\pi }{3}}\right)^{2/3}=5.846665629\mathrm{}.`$ (44)
We see that the coefficient $`b_0`$ for cubic and spherical geometries only differ by about $`3\%`$.

We thus find from eq.(30) in the $`N\mathrm{}`$ limit

$$\text{lim}_N\mathrm{}\frac{1}{N}\mathrm{log}w(\xi ,N)=\frac{3}{2}\mathrm{log}\xi +\frac{9b_0}{2\xi }+\frac{9}{8\xi ^2}\left(b_142b_0^2\right)+𝒪(\xi ^3)$$
(45)
We considered here these integrals in the zero cutoff limit since $`b_0,b_1`$ and $`b_2`$ have finite zero cutoff limits. It is easy to see that their finite cutoff values behave as

$$b_0(a)b_0=𝒪(a^2),b_1(a)b_1=𝒪(a^4),b_2(a)b_2=𝒪(a)$$
(46)
Inserting eq.(45) into eq.(23) yields,

$$g(\xi )=\frac{3}{2}\frac{9b_0}{2\xi }\frac{9}{4\xi ^2}\left(b_142b_0^2\right)+𝒪(\xi ^3)$$
(47)
and

$$\frac{pV}{NT}=1\frac{3b_0}{2\xi }\frac{3}{4\xi ^2}\left(b_142b_0^2\right)+𝒪(\xi ^3).$$
(48)
We see that after letting $`N\mathrm{}`$ the zero cutoff limit is finite. We further discuss this important issue in the next section.

## 3 Statistical Mechanics of the Self-Gravitating Gas: the canonical ensemble

We investigate in this section the self-gravitating gas considered in the previous section but in thermal equilibrium at temperature $`T\beta ^1`$. That is, we work in the canonical ensemble where the system of $`N`$ particles is not isolated but in contact with a thermal bath at temperature $`T`$. We keep assuming the gas being on a cubic box of side $`L`$.

The partition function of the system can be written as

$$𝒵_C(N,T)=\frac{1}{N!}\mathrm{}\underset{l=1}{\overset{N}{}}\frac{d^3p_ld^3q_l}{(2\pi )^3}e^{\beta H_N}$$
(49)
where

$$H_N=\underset{l=1}{\overset{N}{}}\frac{p_l^2}{2m}Gm^2\underset{1l<jN}{}\frac{1}{|\stackrel{}{q}_l\stackrel{}{q}_j|_A}$$
(50)
$`G`$ is Newton’s gravitational constant.

Computing the integrals over the momenta $`p_l,(1lN)`$

$$_{\mathrm{}}^+\mathrm{}\frac{d^3p}{(2\pi )^3}e^{\frac{\beta p^2}{2m}}=\left(\frac{m}{2\pi \beta }\right)^{3/2}$$

yields

$$𝒵_C(N,T)=\frac{1}{N!}\left(\frac{m}{2\pi \beta }\right)^{\frac{3N}{2}}_0^L\mathrm{}_0^L\underset{l=1}{\overset{N}{}}d^3q_le^{\beta Gm^2_{1l<jN}\frac{1}{|\stackrel{}{q}_l\stackrel{}{q}_j|_A}}.$$
(51)
We make now explicit the volume dependence introducing the variables $`\stackrel{}{r}_l,\mathrm{\hspace{0.33em}1}lN`$ defined in eq.(11). The partition function takes then the form,

$$𝒵_C(N,T)=\frac{1}{N!}\left(\frac{mTL^2}{2\pi }\right)^{\frac{3N}{2}}_0^1\mathrm{}_0^1\underset{l=1}{\overset{N}{}}d^3r_le^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)},$$
(52)
where we introduced the dimensionless variable $`\eta `$

$`\eta `$ $``$ $`{\displaystyle \frac{Gm^2N}{LT}}`$ (53)
and $`u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)`$ is defined by eq.(17). Recall that

$$U\frac{Gm^2N}{L}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N),$$
(54)
is the potential energy of the gas.

The free energy takes then the form,

$$F=T\mathrm{log}𝒵_C(N,T)=NT\mathrm{log}\left[\frac{eV}{N}\left(\frac{mT}{2\pi }\right)^{3/2}\right]T\mathrm{\Phi }_N(\eta ),$$
(55)
where

$$\mathrm{\Phi }_N(\eta )=\mathrm{log}_0^1\mathrm{}_0^1\underset{l=1}{\overset{N}{}}d^3r_le^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)},$$
(56)
The derivative of the function $`\mathrm{\Phi }_N(\eta )`$ will be computed by Monte Carlo simulations, mean field methods and, in the weak field limit $`\eta <<1`$, it will be calculated analytically.

We get for the pressure of the gas,

$$p=\left(\frac{F}{V}\right)_T=\frac{NT}{V}\frac{\eta T}{3V}\mathrm{\Phi }_N^{}(\eta ).$$
(57)
\[Here, $`VL^3`$ stands for the volume of the box and $`p`$ is the external pressure on the system.\]. We see from eq.(56) that $`\mathrm{\Phi }_N(\eta )`$ increases with $`\eta `$ since $`u(.)`$ is positive. Therefore, the second term in eq.(57) is a negative correction to the perfect gas pressure $`\frac{NT}{V}`$.

The mean value of the potential energy $`U`$ can be written from eq.(54) as

$$<U>=T\eta \mathrm{\Phi }_N^{}(\eta )$$
(58)
Combining eqs.(57) and (58) yields the virial theorem,

$$\frac{pV}{NT}=1+\frac{<U>}{3NT}\text{or}\frac{pV}{NT}=\frac{1}{2}+\frac{E}{3NT},$$
(59)
where we use that the average value of the kinetic energy of the gas is $`\frac{3}{2}NT`$.

A more explicit form of the equation of state is

$$\frac{pV}{NT}=1\frac{\eta }{3N}\mathrm{\Phi }_N^{}(\eta ),$$
(60)
where

$`\mathrm{\Phi }_N^{}(\eta )`$ $`=`$ $`e^{\mathrm{\Phi }_N(\eta )}{\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}{\displaystyle \underset{l=1}{\overset{N}{}}}d^3r_lu(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)e^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)}`$ (61)
$`=`$ $`{\displaystyle \frac{1}{2}}(N1)e^{\mathrm{\Phi }_N(\eta )}{\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}{\displaystyle \underset{l=1}{\overset{N}{}}}d^3r_l{\displaystyle \frac{1}{|\stackrel{}{r}_1\stackrel{}{r}_2|_a}}e^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)}.`$ (63)
This formula indicates that $`\mathrm{\Phi }_N^{}(\eta )`$ is of order $`N`$ for large $`N`$. Monte Carlo simulations as well as analytic calculations for small $`\eta `$ show that this is indeed the case. In conclusion, we can write the equation of state of the self-gravitating gas as

$$\frac{pV}{NT}=f(\eta ),$$
(64)
where the function

$$f(\eta )1\frac{\eta }{3N}\mathrm{\Phi }_N^{}(\eta ),$$

is independent of $`N`$ for large $`N`$ and fixed $`\eta `$. \[In practice, Monte Carlo simulations show that $`f(\eta )`$ is independent of $`N`$ for $`N>100`$\].

We get in addition,

$$<U>=3NT[1f(\eta )].$$
(65)
In the dilute limit, $`\eta 0`$ and we find the perfect gas value

$$f(0)=1.$$

Equating eqs.(60) and (64) yields,

$$\mathrm{\Phi }_N(\eta )=3N_0^\eta 𝑑x\frac{1f(x)}{x}.$$

Relevant thermodynamic quantities can be expressed in terms of the function $`f(\eta )`$. We find for the free energy from eq.(55),

$$F=F_03NT_0^\eta 𝑑x\frac{1f(x)}{x}.$$
(66)
where

$$F_0=NT\mathrm{log}\left[\frac{eV}{N}\left(\frac{mT}{2\pi }\right)^{3/2}\right]$$
(67)
is the free energy for an ideal gas.

We find for the total energy $`E`$, chemical potential $`\mu `$ and entropy $`S`$ the following expressions,

$`E`$ $`=`$ $`3NT\left[f(\eta ){\displaystyle \frac{1}{2}}\right],`$ (68)
$`\mu `$ $`=`$ $`\left({\displaystyle \frac{F}{N}}\right)_{T,V}=T\mathrm{log}\left[{\displaystyle \frac{V}{N}}\left({\displaystyle \frac{mT}{2\pi }}\right)^{3/2}\right]3T[1f(\eta )]3T{\displaystyle _0^\eta }𝑑x{\displaystyle \frac{1f(x)}{x}},`$ (70)
$`S`$ $`=`$ $`\left({\displaystyle \frac{F}{T}}\right)_V`$ (72)
$`=`$ $`S_0+3N\left[{\displaystyle _0^\eta }𝑑x{\displaystyle \frac{1f(x)}{x}}+f(\eta )1\right],`$ (74)
where

$$S_0=\frac{F_0}{T}+\frac{3}{2}N.$$

is the entropy of the ideal gas.

Notice that here the Gibbs free energy

$$\mathrm{\Phi }=F+pV=F_0+NT\left[f(\eta )3_0^\eta 𝑑x\frac{1f(x)}{x}\right],$$
(75)
is not proportional to the chemical potential. That is, here $`\mathrm{\Phi }\mu N`$ and we have instead,

$$\mathrm{\Phi }\mu N=2NT\left[1f(\eta )\right].$$
(76)
This relationship differs from the customary one (see ) due to the fact that the dilute scaling relation $`NL`$ holds here instead of the usual one $`NL^3`$. The usual relationship is only recovered in the ideal gas limit $`\eta =0`$.

The specific heat at constant volume takes the form,

$`c_V`$ $`=`$ $`{\displaystyle \frac{T}{N}}\left({\displaystyle \frac{S}{T}}\right)_V=3\left[f(\eta )\eta f^{}(\eta ){\displaystyle \frac{1}{2}}\right].`$ (77)
where we used eq.(72). This quantity is also related to the fluctuations of the potential energy $`(\mathrm{\Delta }U)^2`$ and it is positive defined in the canonical ensemble,

$$c_V=\frac{3}{2}+(\mathrm{\Delta }U)^2.$$
(78)
Here,

$$(\mathrm{\Delta }U)^2\frac{<U^2><U>^2}{NT^2}=3\left[f(\eta )\eta f^{}(\eta )1\right].$$
(79)
The specific heat at constant pressure is given by

$$c_P=c_V\frac{T}{N}\frac{\left(\frac{p}{T}\right)_V^2}{\left(\frac{p}{V}\right)_T}.$$
(80)
and then,

$`c_P`$ $`=`$ $`c_V+{\displaystyle \frac{\left[f(\eta )\eta f^{}(\eta )\right]^2}{f(\eta )+\frac{1}{3}\eta f^{}(\eta )}}`$ (81)
$`=`$ $`{\displaystyle \frac{3}{2}}+{\displaystyle \frac{4f(\eta )\left[f(\eta )\eta f^{}(\eta )\right]}{f(\eta )+\frac{1}{3}\eta f^{}(\eta )}}.`$ (83)
The isothermal ($`K_T`$) and adiabatic ($`K_S`$) compressibilities take the form

$`K_T`$ $`=`$ $`{\displaystyle \frac{1}{V}}\left({\displaystyle \frac{V}{p}}\right)_T={\displaystyle \frac{V}{NT}}{\displaystyle \frac{1}{f(\eta )+\frac{1}{3}\eta f^{}(\eta )}},`$ (84)
$`K_S`$ $`=`$ $`{\displaystyle \frac{1}{V}}\left({\displaystyle \frac{V}{p}}\right)_S={\displaystyle \frac{c_V}{c_P}}K_T.`$ (86)
It is then convenient to introduce the compressibilities

$`\kappa _T{\displaystyle \frac{NT}{V}}K_T={\displaystyle \frac{1}{f(\eta )+\frac{1}{3}\eta f^{}(\eta )}}\text{and}\kappa _S{\displaystyle \frac{NT}{V}}K_S={\displaystyle \frac{c_V}{c_P}}\kappa _T,`$ (87)
which are both of order one (intensive) in the $`N,L\mathrm{}`$ limit with $`N/L`$ fixed.

The speed of sound $`v_s`$ can be written as

$$v_s^2=\frac{c_PV^2}{c_VN}\left(\frac{p}{V}\right)_T=\frac{V^2}{N}\left[\frac{T}{Nc_V}\left(\frac{p}{T}\right)_V^2\left(\frac{p}{V}\right)_T\right].$$
(88)
where we used eq.(80) in the last step. Therefore,

$$\frac{v_s^2}{T}=\frac{\left[f(\eta )\eta f^{}(\eta )\right]^2}{3\left[f(\eta )\eta f^{}(\eta )\frac{1}{2}\right]}+f(\eta )+\frac{1}{3}\eta f^{}(\eta ).$$
(89)
The pressure $`p`$ used in this calculation corresponds to the pressure on the surface of the system. Hence, this is the speed of sound on the surface of the system, this is different from the speed of sound inside the volume since the ground state is inhomogeneous. We compute the speed of sound as a function of the point in paper II.

We see that the large $`N`$ limit of the self-gravitating gas is special. Energy, free energy and entropy are extensive magnitudes in the sense that they are proportional to the number of particles $`N`$ (for fixed $`\eta `$). They all depend on the variable $`\eta =\frac{Gm^2N}{LT}`$ which is to be kept fixed for the thermodynamic limit ($`N\mathrm{}`$ and $`V\mathrm{}`$) to exist. Notice that $`\eta `$ contains the ratio $`N/L=NV^{1/3}`$ which must be considered here an intensive variable. Here, the presence of long-range gravitational situations calls for this new intensive variable in the thermodynamic limit.

In addition, all physical magnitudes can be expressed in terms of a single function of one variable: $`f(\eta )`$.

### 3.1 The diluted regime: $`\eta <<1`$

We can obtain the thermodynamic quantities as a series in powers of $`\eta `$ just expanding the exponent in the integrand of $`\mathrm{\Phi }_N(\eta )`$ \[eq.(56)\].

To first order in $`\eta `$ we get,

$`\mathrm{\Phi }_N(\eta )`$ $`=`$ $`\eta {\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}{\displaystyle \underset{l=1}{\overset{N}{}}}d^3r_lu(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)+𝒪(\eta ^2)`$ (90)
$`=`$ $`{\displaystyle \frac{1}{2}}\eta (N1){\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|_a}}+𝒪(a^2)+𝒪(\eta ^2)`$ (92)
$`=`$ $`3(N1)b_0\eta +𝒪(\eta a^2)+𝒪(\eta ^2).`$ (94)
where the coefficient $`b_0`$ is defined by eq.(35).

To first order in $`\eta `$ we see that the cutoff effect is negligible $`𝒪(a^2)`$ \[see (46)\].

To second order in $`\eta `$ we find from eq.(56),

$`e^{\mathrm{\Phi }_N(\eta )}`$ $`=`$ $`{\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}{\displaystyle \underset{l=1}{\overset{N}{}}}d^3r_le^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)}`$ (95)
$`=`$ $`1+3(N1)b_0\eta `$ (97)
$`+`$ $`{\displaystyle \frac{\eta ^2}{2N^2}}[{\displaystyle \frac{N(N1)(N2)(N3)}{4}}{\displaystyle }{\displaystyle \frac{d^3r_1d^3r_2d^3r_3d^3r_4}{|\stackrel{}{r}_1\stackrel{}{r}_2||\stackrel{}{r}_3\stackrel{}{r}_4|}}`$ (99)
$`+`$ $`N(N1)(N2){\displaystyle \frac{d^3r_1d^3r_2d^3r_3}{|\stackrel{}{r}_1\stackrel{}{r}_2||\stackrel{}{r}_1\stackrel{}{r}_3|}}`$ (101)
$`+`$ $`{\displaystyle \frac{N(N1)}{2}}{\displaystyle }{\displaystyle \frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|^2}}]+𝒪(\eta ^3,\eta a^2,\eta ^2a).`$ (103)
where the coefficients in front of the integrals count the number of combinations of particles yielding the same contribution. Using the notation defined by eqs.(35) we get

$`e^{\mathrm{\Phi }_N(\eta )}`$ $`=`$ $`1+3(N1)b_0\eta +\eta ^2[{\displaystyle \frac{9(N1)(N2)(N3)}{2N}}b_0^2`$ (104)
$`+`$ $`{\displaystyle \frac{(N1)(N2)}{2N}}b_1+{\displaystyle \frac{(N1)}{4N}}b_2]+𝒪(\eta ^3,\eta a^2,\eta ^2a).`$ (106)
Taking the log we get in the infinite $`N`$ limit;

$$\text{lim}_N\mathrm{}\frac{1}{N}\mathrm{\Phi }_N(\eta )=3b_0\eta +\eta ^2\left[\frac{1}{2}b_118b_0^2\right]+𝒪(\eta ^3).$$

where we have now set $`a=0`$.

The cutoff effect is here again of order $`𝒪(a^2)`$. It must be noticed that the coefficient $`b_2`$ which has the stronger dependence on the cutoff \[see (46)\] cancels out in the $`N=\mathrm{}`$ limit.

We therefore find in the low density and the large $`N`$ limit using eqs.(60), (64) and (90):

$$\frac{pV}{NT}=f(\eta )=1b_0\eta \eta ^2\left[\frac{1}{3}b_112b_0^2\right]+𝒪(\eta ^3).$$
(107)
Furthermore, the speed of sound approaches for $`\eta 0`$ its perfect gas value,

$$\frac{v_s^2}{T}\stackrel{\eta 0}{=}\frac{5}{3}\frac{4}{3}b_0\eta \frac{5}{9}\eta ^2\left[b_136b_0^2\right]+𝒪(\eta ^3).$$

where we used eqs.(89) and (107).

As we see, there are no divergent contributions in $`\mathrm{\Phi }_N(\eta )`$ in the zero cutoff limit to the second order in $`\eta `$.

At order three a logarithmically divergent integral appears in $`e^{\mathrm{\Phi }_N(\eta )}`$. Namely,

$$\frac{\eta ^3}{3!N^3}\frac{1}{2}N(N1)\frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|_a^3}\frac{\eta ^3}{N}\mathrm{log}a$$

This integral gives to $`f(\eta )`$ and the other physical magnitudes a contribution of the order

$$\frac{\eta ^3}{N^2}\mathrm{log}a$$

Therefore, such quantities can be safely neglected for $`N\mathrm{}`$ and fixed (small) $`a`$ since $`f(\eta )`$ is of order $`N^0`$ for $`N\mathrm{}`$.

More generally, to the nth. order in $`\eta `$ and $`n>3`$ the leading divergent contribution to $`e^{\mathrm{\Phi }_N(\eta )}`$ for $`a0`$ is of the form

$$\frac{\eta ^n}{n!N^n}\frac{1}{2}N(N1)\frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|_a^n}\frac{\eta ^n}{n!N^{n2}}a^{3n}$$

This gives to $`f(\eta )`$ and the other physical magnitudes a contribution of the order

$$\frac{\eta ^n}{n!N^{n1}}a^{3n}$$

As in the $`n=3`$ case, such contributions are negligible in the $`N\mathrm{}`$ limit since we take it at fixed (small) $`a`$.

## 4 Microcanonical vs. Canonical Ensembles

Let us compare the thermodynamical quantities computed in the microcanonical and canonical ensembles in the $`N\mathrm{}`$ limit keeping $`\xi `$ and $`\eta `$ fixed, respectively.

We consider here the dilute limit where we dispose of analytic expressions. The Monte Carlo and mean field results for the two ensembles will be compared in the next sections and in paper II.

In the dilute limit, we have the expressions (48) and (107) for the equation of state in the microcanonical and canonical ensembles, respectively. We want to know whether they are or not equivalent.

Let us start from the microcanonical equation of state (48). We have to express $`\eta `$ in terms of $`\xi `$ in order to compare with the canonical equation of state (107).

It follows from eqs.(16), (24) and (53) that

$$\eta =\frac{g(\xi )}{\xi }$$

Hence, for large $`\xi `$ and small $`\eta `$,

$$\eta =\frac{3}{2\xi }\frac{9b_0}{2\xi ^2}\frac{9}{4\xi ^3}\left(b_142b_0^2\right)+𝒪(\xi ^4)$$
(108)
and then,

$$\frac{1}{\xi }=\frac{2}{3}\eta \left[1+2b_0\eta 2\left(10b_0^2\frac{1}{3}b_1\right)\eta ^2+𝒪(\eta ^3)\right].$$
(109)
One easily sees that inserting eq.(109) in the microcanonical equation of state (48) yields the canonical equation of state (107) \[up to $`𝒪(\eta ^3)=𝒪(\xi ^3)`$\].

Conversely, starting from the canonical ensemble, it follows from eqs.(24), (68) and (107) that

$$\frac{E}{NT}=g(\xi )=3\left[f(\eta )\frac{1}{2}\right]=\frac{3}{2}3b_0\eta \eta ^2\left[b_136b_0^2\right]+𝒪(\eta ^3).$$
(110)
and

$$\xi =\frac{3}{\eta }\left[f(\eta )\frac{1}{2}\right]=\frac{3}{2\eta }\left[12b_0\eta \frac{2}{3}\eta ^2\left(b_136b_0^2\right)+𝒪(\eta ^3)\right].$$

We see that this relation is identical to eqs.(108) and (109) obtained in the microcanonical ensemble \[up to $`𝒪(\eta ^3)=𝒪(\xi ^3)`$\].

Inserting now eq.(108) into the canonical equation of state (107) yields the microcanonical equation of state (48) \[up to $`𝒪(\eta ^3)=𝒪(\xi ^3)`$\].

One verifies in the same way that all thermodynamical quantities coincide to the same order in both ensembles.

In summary, the microcanonical and canonical ensembles yield the same results for $`N\mathrm{}`$ to the orders $`\eta ^0,\eta `$ and $`\eta ^2`$ (or equivalently $`\xi ^0,\xi ^1`$ and $`\xi ^3`$).

## 5 Monte Carlo calculations

We have applied first the standard Metropolis algorithm to the self-gravitating gas in a cube of size $`L`$ in the canonical ensemble at temperature $`T`$. We computed in this way the pressure, the energy, the average density, the potential energy fluctuations, the average particle distance and the average squared particle distance as functions of $`\eta `$. We implement the Metropolis algorithm changing at random the position of one particle chosen at random. The energy of the configurations is calculated performing the exact sums as in eq.(17). We used as statistical weight for the Metropolis algorithm in the canonical ensemble,

$$e^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)},$$

which appears in the coordinate partition function (56).

The number of particles $`N`$ went up to $`2000`$. We introduced a small short distance cutoff $`A=10^4L10^8L`$ in the attractive Newton’s potential according to eq.(2). All results in the gaseous phase were insensitive to the cutoff value. The partition function calculation turns to be much less sensible to the short distance singularities of the gravitational force than Newton’s equations of motion for $`N`$ particles. That is, solving the classical dynamics for $`N`$ particles interacting through gravitational forces as well as solving the Boltzman equation including the $`N`$-body gravitational interaction requires sophisticated algorithms to avoid excessively long computer times . As is clear, solving the $`N`$-body classical evolution or the kinetic equations provides the time-dependent dynamics and out of thermal equilibrium effects which are out of the scope of our approach.

In the CE, two different phases show up: for $`\eta <\eta _T`$ we have a non-perfect gas and for $`\eta >\eta _T`$ it is a condensed system with negative pressure. The transition between the two phases is very sharp. This phase transition is associated with the Jeans instability.

A negative pressure indicates that the free energy grows for increasing volume at constant temperature \[see eq.(57)\]. Therefore, the system wants to contract sucking on the walls.

We plot in figs. 1 and 2, $`f(\eta )=pV/[NT]`$ and $`(\mathrm{\Delta }U)^2`$ as functions of $`\eta `$, respectively.

We find that for small $`\eta `$, the Monte Carlo results for $`pV/[NT]`$ well reproduce the analytical formula (107). $`pV/[NT]`$ monotonically decreases with $`\eta `$.

In the Monte Carlo simulations the phase transition to the condensed phase happens for $`\eta =\eta _T`$ slightly below $`\eta _C`$. For $`N=2000`$ we find $`\eta _T1.515`$. For $`\eta _T<\eta <\eta _C`$, the gaseous phase can only exist as a metastable state.

The average distance between particles $`<r>`$ and the average squared distance between particles $`<r^2>`$ monotonically decrease with $`\eta `$. When the gas collapses at $`\eta _T,<r>`$ and $`<r^2>`$ exhibit a sharp decrease.

The values of $`pV/[NT],<r>`$ and $`<r^2>`$ in the condensed phase are independent of the cutoff for $`a<10^5`$. The Monte Carlo results in this condensed phase can be approximated for $`\eta >2`$ as

$$\frac{pV}{NT}=f(\eta )1K\eta ,<r>0.016.$$
(111)
where $`K14`$.

Since $`f(\eta )`$ has a jump at the transition, the Gibbs free energy $`\mathrm{\Phi }`$ is discontinuous and we have a phase transition of the zeroth order. We find from eq.(75)

$$\frac{\mathrm{\Phi }(\text{collapse})\mathrm{\Phi }(\eta _T)}{NT}=f(\text{collapse})f(\eta _T)21<0.$$
(112)
We can easily compute the latent heat of the transition per particle ($`q`$) using the fact that the volume $`V`$ stays constant. Hence, $`q=\mathrm{\Delta }E/N`$ and we obtain from eq.(68)

$$\frac{q}{T}=\frac{E(\text{collapse})E(\eta _T)}{NT}=3\left[f(\text{collapse})f(\eta _T)\right]23K\eta _T62<0.$$
(113)
This phase transition is different from the usual phase transitions since the two phases cannot coexist in equilibrium as their pressures are different.

Eq.(111) can be understood from the general treatment in sec. III as follows. We have from eqs.(60)-(61)

$$f(\eta )=1\frac{\eta }{3}<\frac{1}{r}>.$$
(114)
The Monte Carlo results indicate that $`<\frac{1}{r}>42`$ is approximately constant in the collapsed region as well as $`<r>`$ and $`<r^2>`$. Eq.(111) thus follows from eq.(114) using such value of $`<\frac{1}{r}>`$.

The behaviour of $`pV/[NT]`$ near $`\eta _C`$ in the gaseous phase can be well reproduced by

$$\frac{pV}{NT}=f(\eta )\stackrel{\eta \eta _C}{=}f_C+A\sqrt{\eta _C\eta }$$
(115)
where $`f_C0.316,A0.414`$ and $`\eta _C1.540`$.

In addition, the behaviour of $`(\mathrm{\Delta }U)^2`$ in the same region is well reproduced by

$$(\mathrm{\Delta }U)^2\stackrel{\eta \eta _C}{=}C+\frac{D}{\sqrt{\eta _C\eta }}$$
(116)
with $`C1.64`$ and $`D0.901`$. \[Notice that for finite $`N,(\mathrm{\Delta }U)^2`$ will be finite albeit very large at the phase transition\]. Eq.(79) relating $`f(\eta )`$ and $`(\mathrm{\Delta }U)^2`$ is satisfied with reasonable approximation.

We thus find a critical region just below $`\eta _C`$ where the energy fluctuations tend to infinity as $`\eta \eta _C`$.

The point $`\eta _T`$ where the phase transition actually takes place in the Monte Carlo simulations is at $`\eta _T1.51<\eta _C`$. This value for $`\eta _T`$ is close to the point $`\eta _0`$ where the isothermal compressibility $`\kappa _T`$ diverges (see sec. VII). They are probably the same point.

Since Monte Carlo simulations are like real experiments, we conclude that the gaseous phase extends from $`\eta =0`$ till $`\eta =\eta _T`$ in the CE and not till $`\eta =\eta _C`$. Notice that in the literature based on the hydrostatic description of the self-gravitating gas , only the instability at $`\eta =\eta _C`$ is discussed whereas the singularities at $`\eta =\eta _0`$ are not considered.

We then performed Monte Carlo calculations in the microcanonical ensemble where the coordinate partition function is given by eq.(18). We thus used

$$\left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right]^{3N/21}\theta \left[\xi +\frac{1}{N}u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)\right],$$

as the statistical weight for the Metropolis algorithm.

The MCE and CE Monte Carlo results coincide up to the statistical error for $`0<\eta <\eta _T`$, that is for $`\mathrm{}>\xi >\xi _T0.19`$. In the MCE the gas does not clump at $`\eta =\eta _C`$ (point $`C`$ in fig. 1) and the specific heat becomes negative between the points $`C`$ and $`MC`$. In the MCE the gas does clump at $`\xi 0.52,\eta _{MC}^T1.33`$ (point $`MC`$ in fig. 1) increasing both its temperature and pressure discontinuously. We find from the Monte Carlo data that the temperature increases by a factor $`2.4`$ whereas the pressure increases by a factor $`3.6`$ when the gas clumps. The transition point $`\eta _{MC}^T`$ in the Monte Carlo simulations is slightly to the right of the critical point $`\eta _{MC}`$ predicted by mean field theory. The mean field yields for the sphere $`\eta _{MC}=1.2598\mathrm{}`$.

In ref. finite $`N`$ corrections to the critical point $`\eta _{MC}`$ are computed in mean field for the sphere. This finite $`N`$ corrections shift $`\eta _{MC}`$ by $`+3.3\%`$ for $`N=2000`$. Since, $`\eta _{MC}^T`$ differs from $`\eta _{MC}`$ by $`+5.6\%`$, $`\eta _{MC}^T`$ and $`\eta _{MC}`$ are probably different critical points.

As is clear, the domain between $`C`$ and $`MC`$ cannot be reached in the CE since $`c_V>0`$ in the CE as shown by eq.(78).

We find an excellent agreement between the Monte Carlo and Mean Field (MF) results (both in the MCE and CE). (This happens although the geometry for the MC calculation is cubic while it is spherical for the MF). The points where the collapse phase transition occurs ($`\eta _T`$ and $`\eta _{MC}^T`$) slowly increase with the number of particles $`N`$.

We verified that the Monte Carlo results in the gaseous phase ($`\eta <\eta _T`$) are cutoff independent for $`10^3a10^7`$.

As for the CE, the Gibbs free energy is discontinuous at the transition in the MCE. The transition is then of the zeroth order. We find from eq.(75)

$$\frac{\mathrm{\Phi }(\text{collapse})\mathrm{\Phi }(\eta _T)}{NT_{gas}}=\frac{T_{coll}}{T_{gas}}f(\text{collapse})f(\eta _T)0.7>0.$$

where we used the numerical values from the Monte Carlo simulations. Notice that the Gibbs free energy increases at the MC transition whereas it decreases at the C transition \[see eq.(112)\].

Here again the two phases cannot coexist in equilibrium since their pressures and temperatures are different.

We display in figs. 3-4 the average particle distribution from Monte Carlo simulations with $`2000`$ particles in the microcanonical ensemble at both sides of the gravothermal catastrophe, i. e. $`\eta =\eta _{MC}`$. Fig. 3 corresponds to the gaseous phase and fig. 4 to the collapsed phase. The inhomogeneous particle distribution is clear in fig. 3 whereas fig. 4 shows a dense collapsed core surrounded by a halo of particles.

The different nature of the collapse in the CE and in the MCE can be explained using the virial theorem \[see eq.(59)\]

$$\frac{pV}{NT}=1+\frac{U}{NT}.$$

When the gas collapses in the CE the particles get very close and $`U`$ becomes large and negative while $`T`$ is fixed. Therefore, $`\frac{pV}{NT}`$ may become large and negative as it does.

We can write the virial theorem also as,

$$pV\frac{1}{2}NT=\frac{1}{3}E.$$

When the gas is near the point MC, $`E<0`$ is fixed and we have $`T>0`$. Therefore, $`\frac{pV}{NT}`$ as well as $`U=E3NT/2`$ cannot become large and negative as in the CE collapse. This prevents the distance between the particles to decrease. Actually, the Monte Carlo simulations show that $`<r>`$ increases by $`18\%`$ when the gas collapses in the MCE.

Figs. 5 and 6 depict the average particle distribution from Monte Carlo simulations with $`2000`$ particles in the canonical ensemble at both sides of the collapse critical point, i. e. $`\eta =\eta _C`$. Fig. 5 corresponds to the gaseous phase and fig. 6 to the collapsed phase. The inhomogeneous particle distribution is clear in fig. 5 whereas fig. 6 shows a dense collapsed core surrounded by a very little halo of particles.

Notice that the collapsed phases are of different nature in the CE and MCE. The core is much tighter and the halo much smaller in the CE than in the MCE.

Figs. 3 and 5 depict the average particle distribution for the gaseous phase in the MCE and the CE, respectively. In this phase, the MC simulations give identical descriptions for large $`N`$ in both ensembles. \[This important point will be further demostrated in sec. VI by functional integral methods\]. The average configurations in figs. 3 and 5 describe a self-gravitating gas in thermal equilibrium within a cube. We may call it the isothermal cube by analogy with the well known isothermal sphere-.

## 6 Mean Field Approach

Both in the microcanonical and the canonical ensembles the coordinate partition functions are given by $`3N`$-uple integrals \[eqs.(18) and (56), respectively\]. In the $`N\mathrm{}`$ limit both $`3N`$-uple integrals can be recasted as functional integrals over the continuous particle density as we see below.

### 6.1 The Canonical Ensemble

We now recast the coordinate partition function $`e^{\mathrm{\Phi }_N(\eta )}`$ in the canonical ensemble given by eq.(56) as a functional integral in the thermodynamic limit.

$`e^{\mathrm{\Phi }_N(\eta )}`$ $`\stackrel{N>>1}{=}`$ $`{\displaystyle D\rho 𝑑\widehat{a}e^{Ns_C[\rho (.),\widehat{a},\eta ]}}`$ (117)
$`s_C[\rho (.),\widehat{a},\eta ]`$ $`=`$ $`{\displaystyle \frac{\eta }{2}}{\displaystyle \frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}\rho (\stackrel{}{x})\rho (\stackrel{}{y})}+{\displaystyle d^3x\rho (\stackrel{}{x})\mathrm{log}\rho (\stackrel{}{x})}i\widehat{a}\left({\displaystyle d^3x\rho (\stackrel{}{x})}1\right).`$
where we used the coordinates $`\stackrel{}{x}`$ in the unit volume. The first term is the potential energy, the second term is the functional integration measure for this case (see appendix A). Here $`N\rho (\stackrel{}{x})`$ stands for the density of particles.

The integration over $`\widehat{a}`$ enforces the number of particles to be exactly $`N`$:

$$d^3x\rho (\stackrel{}{x})=1$$
(119)
That is, in the coordinates $`\stackrel{}{q}`$ (running from $`0`$ to $`L`$), the density of particles is

$$\frac{N}{L^3}\rho (\stackrel{}{q})\mathrm{with}d^3q\frac{N}{L^3}\rho (\stackrel{}{q})=N.$$

### 6.2 The Microcanonical Ensemble

Let us express the coordinate partition function in the microcanonical ensemble $`w(\xi ,N)`$ defined by eq.(18) in terms of the coordinate partition function in the canonical ensemble $`e^{\mathrm{\Phi }_N(\eta )}`$ defined by eq.(56). In order to do that we use the Fourier expansion

$$x^\lambda \theta (x)=\frac{\mathrm{\Gamma }(\lambda +1)}{2\pi }_{\mathrm{}}^+\mathrm{}e^{i\omega x}\frac{d\omega }{(i\omega )^{\lambda +1}}$$
(120)
We thus find from eqs.(18), (56) and (120) that

$`w(\xi ,N)`$ $`=`$ $`\mathrm{\Gamma }\left({\displaystyle \frac{3N}{2}}\right){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\omega }{2\pi }}e^{i\omega \xi +\mathrm{\Phi }_N(i\omega /N)\frac{3N}{2}\mathrm{log}(i\omega )}`$ (121)
$`=`$ $`N\mathrm{\Gamma }\left({\displaystyle \frac{3N}{2}}\right){\displaystyle _\gamma }{\displaystyle \frac{d\eta }{2\pi i}}e^{N\eta \xi +\mathrm{\Phi }_N(\eta )\frac{3N}{2}\mathrm{log}(N\eta )}`$ (123)
where we introduced the integration variable $`\eta i\omega /N`$ and where $`\gamma `$ is an upward integration contour parallel to the imaginary $`\eta `$ axis. Using Stirling’s approximation for the $`\mathrm{\Gamma }`$ function we find for $`N1`$ up to irrelevant constants

$$w(\xi ,N)\stackrel{N>>1}{=}_\gamma \frac{d\eta }{2\pi i}e^{N\eta \xi +\mathrm{\Phi }_N(\eta )\frac{3N}{2}\mathrm{log}\eta }$$

Now, inserting the functional integral representation (117) for the coordinate canonical partition function yields,

$$w(\xi ,N)\stackrel{N>>1}{=}D\rho 𝑑\widehat{a}\frac{d\eta }{2\pi i}e^{N[\eta \xi \frac{3}{2}\mathrm{log}\eta s_C[\rho (.),\widehat{a},\eta ]]}$$
(124)
We thus find a functional integral representation in the microcanonical ensemble analogous to the canonical representation eq.(117) but with an extra integration (over $`\eta `$) that constrains the value of the energy.

The ‘effective action’ in the microcanonical ensemble takes thus the form,

$$s_{MC}[\rho (.),\widehat{a},\eta ]=\frac{3}{2}\mathrm{log}\eta \eta \xi \frac{\eta }{2}\frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}\rho (\stackrel{}{x})\rho (\stackrel{}{y})+d^3x\rho (\stackrel{}{x})\mathrm{log}\rho (\stackrel{}{x})i\widehat{a}(d^3x\rho (\stackrel{}{x})1)$$
(125)
### 6.3 The Grand Canonical Ensemble

The partition function in the grand canonical ensemble can be written as

$$𝒵_{GC}(z,T)=\underset{N=0}{\overset{\mathrm{}}{}}z^N𝒵(N,T),$$
(126)
where $`z=e^{\frac{\mu }{T}}`$ stands for the fugacity and $`𝒵(N,T)`$ is the partition function in the canonical ensemble given by eqs.(49) and (52).

As shown in ref., this grand canonical partition function can be expressed as a functional integral

$$𝒵_{GC}(z,T)=𝒟\mathrm{\Phi }e^{\frac{1}{T_{eff}}_Vd^3q\left[\frac{1}{2}\mathrm{\Phi }^2\mathrm{\Phi }+M^2e^{\mathrm{\Phi }(\stackrel{}{q})}\right]},$$
(127)
where

$$M^2=\sqrt{\frac{2T}{\pi }}zGm^{7/2},T_{eff}=4\pi \frac{Gm^2}{T}.$$
(128)
Notice that the representation (127) is exact while the functional integral representations in the microcanonical and canonical ensembles only apply for large number of particles.

Rewriting eq.(127) in terms of the dimensionless variables (11) yields for the exponent

$$\frac{1}{T_{eff}}_Vd^3q\left[\frac{1}{2}\mathrm{\Phi }_q^2\mathrm{\Phi }+\mu ^2e^{\mathrm{\Phi }(\stackrel{}{q})}\right]=\frac{L}{T_{eff}}_0^1d^3x\left[\frac{1}{2}\mathrm{\Phi }_r^2\mathrm{\Phi }+\zeta ^2e^{\mathrm{\Phi }(\stackrel{}{x})}\right]$$

where $`\zeta ML`$ is of the order one ($`L^0`$), since $`M^2z=e^{\frac{\mu }{T}}L^2`$ \[see eq.(70)\].

Since the exponent in the functional integral (127) is proportional to $`L`$, the large volume limit is dominated by the stationary points (mean field approximation)

$$_r^2\mathrm{\Phi }_s(\stackrel{}{x})+\zeta ^2e^{\mathrm{\Phi }_s(\stackrel{}{x})}=0.$$
(129)
We expand around the saddle point $`\mathrm{\Phi }_s(\stackrel{}{x})`$ changing to a new functional integration $`Y(\stackrel{}{x})`$ variable as follows,

$$\mathrm{\Phi }(\stackrel{}{x})=\mathrm{\Phi }_s(\stackrel{}{x})+Y(\stackrel{}{x}).$$
(130)
Keeping in eq.(127) quadratic terms in $`Y(.)`$ yields,

$$𝒵_{GC}(z,T)=e^{\frac{L}{T_{eff}}_0^1d^3x\left[\frac{1}{2}\mathrm{\Phi }_s_r^2\mathrm{\Phi }_s+\zeta ^2e^{\mathrm{\Phi }_s(\stackrel{}{x})}\right]}𝒟Ye^{\frac{L}{2T_{eff}}_0^1d^3x\left[Y^2Y+\zeta ^2Y^2e^{\mathrm{\Phi }_s(\stackrel{}{x})}\right]}\left[1+𝒪\left(\frac{1}{L}\right)\right]$$
(131)
where the Gaussian integral over $`Y(.)`$ gives a factor of order $`L^0`$ \[see paper II\].

We recall that the saddle point method applies while all eigenvalues of the quadratic form in the exponent of eq.(131) are positive. Therefore, the determinant of the quadratic fluctuations is positive. The determinant vanishing or changing sign indicates the presence of zero or negative eigenvalues. In such a case the system is no more on a stable phase but on a metastable or unstable phase. The free energy gets an imaginary part in such metastable or unstable situations.

The average number of particles in the grand canonical ensemble is given by

$$\overline{N}=\frac{1}{𝒵_{GC}}\underset{N=0}{\overset{\mathrm{}}{}}Nz^N𝒵(N,T)=\frac{\mathrm{log}𝒵_{GC}}{\mathrm{log}z}|_{V,T}.$$

We thus find in the mean field approximation,

$$\overline{N}=\frac{L\zeta ^2}{T_{eff}}_0^1d^3xe^{\mathrm{\Phi }_s(\stackrel{}{x})}$$

Therefore, using this and eq.(128) we can express $`\zeta ^2`$ in terms of $`\eta `$ where we denote $`\overline{N}`$ as $`N`$ to avoid cluttering of notation,

$$\zeta ^2=\frac{4\pi \eta }{_0^1d^3xe^{\mathrm{\Phi }_s(\stackrel{}{x})}},$$
(132)
and the fugacity results

$$z=\frac{N}{L^3}\left(\frac{2\pi }{mT}\right)^{3/2}\frac{1}{_0^1d^3xe^{\mathrm{\Phi }_s(\stackrel{}{x})}}.$$
(133)
We again see that $`zL^2`$ in the GCE.

Integrating eq.(129) over the unit volume yields

$$\stackrel{}{}\mathrm{\Phi }_s(\stackrel{}{x})𝑑\stackrel{}{s}=4\pi \eta $$
(134)
where we used eq.(132).

We find for the free energy,

$$F=T\mathrm{log}𝒵_{GC}+NT\mathrm{log}z=F_0+\frac{NT}{2}K(\eta )NT\mathrm{log}C(\eta )+𝒪(N^0),$$
(135)
where we used the grand canonical partition function (131) evaluated at the stationary point,

$$\mathrm{log}𝒵_{GC}=N\left[1\frac{1}{2}K(\eta )\right],$$
(136)
and $`z`$ is given by eq.(133) with

$$K(\eta )\frac{_0^1d^3x\mathrm{\Phi }_s(\stackrel{}{x})e^{\mathrm{\Phi }_s(\stackrel{}{x})}}{C(\eta )}\text{and}C(\eta )_0^1d^3xe^{\mathrm{\Phi }_s(\stackrel{}{x})}.$$
(137)
$`F_0`$ is given by eq.(67).

We easily calculate the mean value of the potential energy in the mean field approximation

$$<U>=T\frac{\mathrm{log}𝒵_{GC}}{\mathrm{log}G}=\frac{NT}{2}K(\eta )$$
(138)
Combining the two expressions for the entropy

$$S=\frac{EF}{T}\text{and}S=\left(\frac{F}{T}\right)_V,$$
(139)
yields

$$S=S_0N\left[K(\eta )\mathrm{log}C(\eta )\right]$$
(140)
and the first order differential equation

$$\eta K^{}(\eta )+K(\eta )=2\eta \frac{d}{d\eta }\mathrm{log}C(\eta ).$$
(141)
The boundary conditions $`K(0)=0,C(0)=1`$ ensure the ideal gas limit $`\eta =0`$.

The pressure takes the form,

$$P=\left(\frac{F}{V}\right)_T=\frac{NT}{V}\left[1+\frac{\eta }{3}\left(\frac{1}{2}K^{}(\eta )\frac{d}{d\eta }\mathrm{log}C(\eta )\right)\right]+𝒪(N^0).$$
(142)
These equations guarantee in addition that the virial theorem (59) holds.

### 6.4 Saddle point evaluation in the canonical ensemble

The functional integral in eq.(117) is dominated for large $`N`$ by the extrema of the ‘effective action’ $`s_C[\rho (.),\widehat{a},\eta ]`$, that is, the solutions of the stationary point equation

$$\mathrm{log}\rho _s(\stackrel{}{x})\eta \frac{d^3y\rho _s(\stackrel{}{y})}{|\stackrel{}{x}\stackrel{}{y}|}=a_s,$$
(143)
$`a=i\widehat{a}`$ is a Lagrange multiplier enforcing the constraint (119).

Applying the Laplacian and setting $`\varphi (\stackrel{}{x})\mathrm{log}\rho _s(\stackrel{}{x})`$ yields,

$$^2\varphi (\stackrel{}{x})+4\pi \eta e^{\varphi (\stackrel{}{x})}=0,$$
(144)
This equation is scale covariant . That is, if $`\varphi (\stackrel{}{x})`$ is a solution of eq.(144), then

$$\varphi _\lambda (\stackrel{}{x})\varphi (\lambda \stackrel{}{x})+\mathrm{log}\lambda ^2$$
(145)
where $`\lambda `$ is an arbitrary constant is also a solution of eq.(144). For spherically symmetric solutions this property can be found in ref..

Integrating eq.(144) over the unit volume and using the constraint (119) yields

$$\stackrel{}{}\varphi (\stackrel{}{x})𝑑\stackrel{}{s}=4\pi \eta $$
(146)
where the surface integral is over the boundary of the unit volume.

Comparing eqs.(129)-(134) with (144) and (146) shows that the grand canonical and canonical stationary points are related by

$$\mathrm{\Phi }_s(\stackrel{}{x})=\varphi (\stackrel{}{x})+\mathrm{log}C(\eta ).$$
(147)
Eq.(131) can then be written as

$`𝒵_{GC}(z,T)=e^{\frac{N}{4\pi \eta }\left\{_0^1d^3x\left[\frac{1}{2}\varphi _r^2\varphi +4\pi \eta e^{\varphi (\stackrel{}{x})}\right]2\pi \eta \mathrm{log}C(\eta )\right\}}\times `$ (148)
(149)
$`{\displaystyle 𝒟Ye^{\frac{N}{8\pi \eta }_0^1d^3x\left[Y^2Y+4\pi \eta Y^2e^{\varphi (\stackrel{}{x})}\right]}\left[1+𝒪\left(\frac{1}{N}\right)\right]}`$ (150)
where we used eqs.(53), (128), (132) , (146) and (147).

We have taken the zero cutoff limit in eqs.(143)-(144). The mean field equations turn to be finite with regular solutions in such limit. This can be understood from our perturbative calculation in sec. III.A. All potentially divergent contributions at zero cutoff are suppressed by factors $`1/N^2`$ and therefore disappear in the $`N=\mathrm{}`$ limit. Hence one can set the cutoff to zero in the mean field approximation.

In order to evaluate the functional integral in eq.(117) by the saddle point method we change the functional integration variable as follows,

$$\rho (\stackrel{}{x})=\rho _s(\stackrel{}{x})+Y(\stackrel{}{x}),a=a_s+y_0$$
(151)
where $`\rho _s(\stackrel{}{x})`$ and $`a_s`$ obey eq.(143). We can expand the exponent to second order as

$$s_C[\rho (.),a,\eta ]s_C[\rho _s(.),a_s,\eta ]=s_C^{(2)}[Y(.),y_0]+𝒪(Y^3,Y^2y_0)$$
(152)
where we use that

$$\frac{\delta s_C}{\delta \rho (\stackrel{}{x})}|_{\rho =\rho _s,a=a_s}=0,\frac{s_C}{a}|_{\rho =\rho _s,a=a_s}=0$$

and

$$s_C^{(2)}[Y(.),y_0]\frac{1}{2}d^3xd^3yY(\stackrel{}{x})Y(\stackrel{}{y})\frac{\delta ^2s_C}{\delta \rho (\stackrel{}{x})\delta \rho (\stackrel{}{y})}|_{\rho =\rho _s,a=a_s}+y_0d^3xY(\stackrel{}{x})\frac{\delta ^2s_C}{\delta \rho (\stackrel{}{x})a}|_{\rho =\rho _s,a=a_s}.$$

Notice that

$$\frac{^2s_C}{a^2}=0.$$

We evaluate explicitly the second derivatives from eq.(117) with the result,

$$\frac{\delta ^2s_C}{\delta \rho (\stackrel{}{x})\delta \rho (\stackrel{}{y})}=\frac{\delta (\stackrel{}{x}\stackrel{}{y})}{\rho (\stackrel{}{x})}\frac{\eta }{|\stackrel{}{x}\stackrel{}{y}|},\frac{\delta ^2s_C}{\delta \rho (\stackrel{}{x})a}=1$$

Therefore,

$$s_C^{(2)}[Y(.),y_0]=\frac{1}{2}d^3x\frac{Y^2(\stackrel{}{x})}{\rho (\stackrel{}{x})}\frac{\eta }{2}\frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}Y(\stackrel{}{x})Y(\stackrel{}{y})y_0d^3xY(\stackrel{}{x})$$
(153)
Inserting eqs.(151) and (152) into eq.(117) yields

$$e^{\mathrm{\Phi }_N(\eta )}\stackrel{N>>1}{=}e^{Ns(\eta )}DY𝑑y_0e^{Ns_C^{(2)}[Y(.),y_0]}\left[1+𝒪\left(\frac{1}{N}\right)\right]$$
(154)
where $`s(\eta )s_C[\rho _s(.),a_s,\eta ]`$ stands for the value of the exponent at the saddle point. Terms of order higher than quadratic in $`s_C[\rho (.),a,\eta ]`$ contribute to the $`1/N`$ corrections.

The Gaussian functional integral (154) can be exactly computed in terms of the functional determinant of the quadratic form (153) \[see paper II\]. It gives a result of order one ($`N^0`$).

In the mean field approximation we only keep the dominant order for large $`N`$. Therefore, only the exponent at the saddle point accounts and according to eq.(55) we find for the free energy

$`F`$ $`=`$ $`F_0+NTs(\eta )+𝒪(N^0)`$ (155)
$`{\displaystyle \frac{pV}{NT}}`$ $`=`$ $`1+{\displaystyle \frac{\eta }{3}}{\displaystyle \frac{ds}{d\eta }}+𝒪(N^1)`$ (157)
Hence, in the mean field approximation, the function $`f(\eta )`$ is given by

$$f_{MF}(\eta )1+\frac{\eta }{3}\frac{ds}{d\eta },$$
(158)
From eq.(117) we can compute $`s(\eta )`$ in terms of the saddle point solution as follows

$$s(\eta )s_C[\rho _s(.),a_s,\eta ]=\frac{\eta }{2}\frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}\rho _s(\stackrel{}{x})\rho _s(\stackrel{}{y})+d^3x\rho _s(\stackrel{}{x})\mathrm{log}\rho _s(\stackrel{}{x}).$$
(159)
Using eq.(143) we find an equivalent expression that will be useful in paper II,

$$s(\eta )=\frac{a_s}{2}+\frac{1}{2}\varphi (\stackrel{}{x})e^{\varphi (\stackrel{}{x})}d^3x.$$
(160)
### 6.5 Saddle point evaluation in the microcanonical ensemble

The extrema of the ‘effective action’ (125) dominate the microcanonical partition function (124) in the large $`N`$ limit. Extremizing eq.(125) with respect to $`\rho (.)`$ and $`\widehat{a}`$ gives again eqs.(143) and (119), respectively.

An extra equation follows by extremizing the ‘effective action’ (125) with respect to $`\eta `$:

$$\xi =\frac{3}{2\eta _s}\frac{1}{2}\frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}\rho _s(\stackrel{}{x})\rho _s(\stackrel{}{y})$$
(161)
Going back to dimensionful variables this equation takes the familiar form

$$E=\frac{3}{2}NT\frac{Gm^2}{2}\frac{d^3qd^3q^{}}{|\stackrel{}{q}\stackrel{}{q}^{}|}N\rho (\stackrel{}{q})N\rho (\stackrel{}{q}^{})$$

That is, eq.(161) enforces the fixed value of the energy in the microcanonical ensemble.

Therefore, the stationary point equations in the microcanonical and canonical ensembles are identical. Both ensembles yield the same results in the $`N\mathrm{}`$ limit in their common region of validity. We derive the domain of validity of the mean field approach for each of the three statistical ensembles in paper II. That is, the regions where all fluctuations around it decrease its statistical weight within their common region of validity.

In order to evaluate the functional integral for the microcanonical partition function (124)

$$w(\xi ,N)\stackrel{N>>1}{=}D\rho 𝑑\widehat{a}\frac{d\eta }{2\pi i}e^{Ns_{MC}[\rho (.),\widehat{a},\eta ]}$$
(162)
we expand the ‘effective action’ $`s_{MC}[\rho (.),\widehat{a},\eta ]`$ around the stationary point $`\rho _s(.),\widehat{a}_s,\eta _s`$ to second order. This gives

$$w(\xi ,N)\stackrel{N>>1}{=}e^{Ns(\eta )}DY𝑑y_0\frac{d\stackrel{~}{\eta }}{2\pi i}e^{Ns_{MC}^{(2)}[Y(.),y_0,\stackrel{~}{\eta }]}\left[1+𝒪\left(\frac{1}{N}\right)\right]$$
(163)
where $`Y(.)`$ and $`y_0`$ are defined by eq.(151) and we set $`\eta =\eta _s+\stackrel{~}{\eta }`$. The second order piece of the ‘effective action’ takes now the form

$$s_{MC}^{(2)}[Y(.),y_0,\stackrel{~}{\eta }]=s_C^{(2)}[Y(.),y_0]\stackrel{~}{\eta }\frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}\rho _s(\stackrel{}{x})Y(\stackrel{}{y})\frac{3}{4\eta _s^2}\stackrel{~}{\eta }^2.$$
(164)
The Gaussian functional integral in eq.(163) yields a contribution of order one ($`N^0`$) \[see paper II\]. The dominant (mean field) contribution, $`e^{Ns(\eta )}`$, exactly coincides with the mean field result in the canonical ensemble \[eq.(154)\] Therefore, the canonical and microcanonical ensembles yields identical physical magnitudes and the same equation of state in the mean field limit.

### 6.6 Spherically symmetric case

We shall consider here the spherically symmetric case where eq.(144) takes the form

$$\frac{d^2\varphi }{dR^2}+\frac{2}{R}\frac{d\varphi }{dR}+4\pi \eta e^{\varphi (R)}=0.$$
(165)
where we work on an unit volume sphere instead of an unit volume cube as in eq.(11). Therefore, the radial variable runs in the interval

$$0R\left(\frac{3}{4\pi }\right)^{1/3}$$

It is more convenient to introduce a new radial variable

$$rR\left(\frac{4\pi }{3}\right)^{1/3}$$

such that $`0r1`$.

The saddle point equation (165) takes then the form

$$\frac{d^2\varphi }{dr^2}+\frac{2}{r}\frac{d\varphi }{dr}+4\pi \eta ^Re^{\varphi (r)}=0.$$
(166)
where

$$\eta ^R\eta \left(\frac{4\pi }{3}\right)^{1/3}=1.61199\mathrm{}\eta \text{and}e^{\varphi (r)}=e^{\varphi (R)}\frac{3}{4\pi }.$$
(167)
In order to have a regular solution at $`r=0`$ one has to impose

$$\varphi ^{}(0)=0.$$
(168)
Otherwise, the second term in eq.(165) diverges for $`r0`$.

In the spherically symmetric case, the constraint (146) becomes

$$\varphi ^{}(1)=\eta ^R.$$
(169)
Using the scale covariance (145) we can express $`\varphi (r)`$ as

$$\varphi (r)=\mathrm{log}\left(\frac{\lambda ^2}{4\pi \eta ^R}\right)+\chi (\lambda r)$$
(170)
where

$$\chi ^{\prime \prime }(\lambda )+\frac{2}{\lambda }\chi ^{}(\lambda )+e^{\chi (\lambda )}=0,\chi ^{}(0)=0$$
(171)
This equation is invariant under the transformation:

$$\lambda \lambda e^\alpha \chi (\lambda )\chi (\lambda )2\alpha ,$$
(172)
where $`\alpha `$ is a real number. Hence, we can set $`\chi (0)0`$ without loosing generality.

$`\chi (x)`$ is independent of $`\eta ^R`$, and $`\lambda `$ is related to $`\eta ^R`$ through eq.(169)

$$\lambda \chi ^{}(\lambda )=\eta ^R.$$
(173)
Since $`\lambda `$ and $`\eta ^R`$ are always positive, $`\chi (\lambda )`$ is a monotonically decreasing function of $`\lambda `$.

Eq.(171) can be easily solved for small arguments as

$$\chi (x)=\frac{x^2}{6}+\frac{x^4}{120}+𝒪(x^6)$$

Hence, in the dilute limit eq.(173) relating $`\eta ^R`$ with $`\lambda `$ gives

$$\eta ^R=\frac{\lambda ^2}{3}\frac{\lambda ^4}{30}+𝒪(\lambda ^6).$$
(174)
For large argument, the solution of eq.(171) takes the asymptotic form

$$\chi (x)=\mathrm{log}\frac{2}{x^2}+\frac{A}{\sqrt{x}}\mathrm{cos}\left(\frac{\sqrt{7}}{2}\mathrm{log}x+B\right)\left[1+𝒪\left(\frac{1}{x}\right)\right]$$
(175)
where $`A`$ and $`B`$ are numerical constants. Using eq.(173) this gives for $`\eta ^R`$

$$\eta ^R=2+\frac{C}{\sqrt{\lambda }}\mathrm{cos}\left(\frac{\sqrt{7}}{2}\mathrm{log}\lambda +D\right)\left[1+𝒪\left(\frac{1}{\lambda }\right)\right]$$
(176)
where $`C`$ and $`D`$ are constants related to $`A`$ and $`B`$. By numerically solving eq.(171) we find

$$C=1.667\mathrm{}.$$

It must be noticed, however, that the mean field solution is unphysical for $`\lambda >\lambda _{MC}=34.36361\mathrm{}`$ as we shall see in paper II. Anyway, we see from fig. 7 that $`\eta ^R`$ approaches very fast its asymptotic behaviour (176) for $`\mathrm{log}\lambda >2`$.

We plot in fig. 8 $`\chi (\lambda (\eta ^R))`$ as a function of $`\eta ^R`$.

In the spherically symmetric case the integral over the angles in eq.(143) is immediate with the result,

$$\varphi (r)=a_s+4\pi \eta ^R\left[\frac{1}{r}_0^rr^{\mathrm{\hspace{0.17em}2}}𝑑r^{}e^{\varphi (r^{})}+_r^1r^{}𝑑r^{}e^{\varphi (r^{})}\right]$$
(177)
Deriving with respect to $`r`$ yields,

$$\frac{d\varphi (r)}{dr}=\frac{4\pi \eta ^R}{r^2}_0^rr_{}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}𝑑r^{}e^{\varphi (r^{})}$$

This again shows that $`\varphi (r)`$ is a monotonically decreasing function of $`r`$ \[see above, eq.(173)\].

Setting $`r=1`$ in eq.(177) leads to the relation

$$\varphi (1)=a_s+4\pi \eta ^R_0^1r^2𝑑re^{\varphi (r)}$$

Using now the constraint (119) allows us to compute the Lagrange multiplier $`a`$ at the saddle point

$$a_s=\varphi (1)\eta ^R$$
(178)
The particle density in MF is given by

$$\rho (r)=e^{\varphi (r)}=\frac{\lambda ^2}{4\pi \eta ^R}e^{\chi (\lambda r)},\mathrm{\hspace{0.33em}0}r1.$$

Since $`\chi (\lambda )`$ monotonically decreases with $`\lambda `$, the particle density monotonically decreases with $`r`$ for fixed $`\eta ^R`$.

Let us now compute $`s(\eta ^R)`$ \[the exponent in eq.(117) at the saddle point\] for the spherically symmetric case. We find from eq.(160)

$`s(\eta ^R)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\varphi (1)\eta ^R\right]+2\pi {\displaystyle _0^1}r^2𝑑r\varphi (r)e^{\varphi (r)}`$ (179)
$`=`$ $`\mathrm{log}\left({\displaystyle \frac{\lambda ^2}{4\pi \eta ^R}}\right)+\chi (\lambda ){\displaystyle \frac{\eta ^R}{2}}+{\displaystyle \frac{1}{2\lambda \eta ^R}}{\displaystyle _0^\lambda }x^2𝑑x[\chi ^{}(x)]^2`$ (181)
where we integrated by parts and used eqs.(170)-(173).

The integral in the r.h.s. of eq.(179) can be computed in closed form \[see appendix B\] with the result,

$$s(\eta ^R)=\mathrm{log}\left(\frac{\lambda ^2}{4\pi \eta ^R}\right)+\chi (\lambda )+3\eta ^R\frac{\lambda ^2}{\eta ^R}e^{\chi (\lambda )}$$

Inserting now $`s(\eta ^R)`$ into eq.(158) and using eqs.(171)-(173) yields after calculation

$`f_{MF}(\eta ^R)`$ $`=`$ $`{\displaystyle \frac{\lambda ^2}{3\eta ^R}}e^{\chi (\lambda )},`$ (182)
$`s(\eta ^R)`$ $`=`$ $`3[1f_{MF}(\eta ^R)]\eta ^R+\mathrm{log}\left[{\displaystyle \frac{3f_{MF}(\eta ^R)}{4\pi }}\right]`$
Notice that $`f_{MF}(\eta ^R)`$ as well as the other physical quantities are invariant under the transformation (172) as it must be.

It follows from eqs.(171), (173) and (182) that $`f_{MF}(\eta ^R)`$ obeys the first order non-linear differential equation

$$\eta ^R(3f_{MF}1)f_{MF}^{}(\eta ^R)+(3f_{MF}3+\eta ^R)f_{MF}=0.$$
(184)
which reduces to an Abel equation of first kind.

We thus find that in the mean field approximation all thermodynamic quantities follow from the resolution of the single first order non-linear differential eq.(184) with the initial condition $`f_{MF}(0)=1`$.

Integrating eq.(184) with respect to $`\eta ^R`$ yields,

$$3_0^{\eta ^R}\frac{dx}{x}[1f_{MF}(x)]=3[f_{MF}(\eta ^R)1]+\eta ^R\mathrm{log}f_{MF}(\eta ^R)$$

Further useful relations follow from eqs.(170) and (182)

$$\varphi (1)=\mathrm{log}\left[\frac{3f_{MF}(\eta ^R)}{4\pi }\right],\rho (1)=\frac{3}{4\pi }f_{MF}(\eta ^R).$$
(185)
That is, the particle density at the surface ($`r=1`$) is proportional to $`f_{MF}(\eta ^R)`$.

We can then write the different physical magnitudes in the MF approximation as

$`{\displaystyle \frac{pV}{NT}}`$ $`=`$ $`f_{MF}(\eta ^R)`$ (186)
$`{\displaystyle \frac{FF_0}{NT}}`$ $`=`$ $`3[1f_{MF}(\eta ^R)]\eta ^R+\mathrm{log}f_{MF}(\eta ^R)`$ (188)
$`{\displaystyle \frac{SS_0}{N}}`$ $`=`$ $`6[f_{MF}(\eta ^R)1]+\eta ^R\mathrm{log}f_{MF}(\eta ^R)`$ (190)
$`{\displaystyle \frac{E}{NT}}`$ $`=`$ $`3[f_{MF}(\eta ^R){\displaystyle \frac{1}{2}}]`$
where we used eqs.(64), (66), (68) and (72).

We derive in appendix C the properties of the function $`f_{MF}(\eta ^R)`$ from the differential equation (184). One easily obtains for small $`\eta ^R`$ (dilute regime),

$$f_{MF}(\eta ^R)=1\frac{\eta ^R}{5}\frac{(\eta ^R)^2}{175}+𝒪([\eta ^R]^3).$$

These terms exactly coincide with the perturbative calculation in the dilute regime for spherical symmetry \[see eq.(40), (107) and (167)\].

We plot in fig. 1 $`f_{MF}(\eta ^R)`$ as a function of $`\eta ^R`$ obtained by solving eq.(184) by the Runge-Kutta method. We see that $`f_{MF}(\eta ^R)`$ is a monotonically decreasing function of $`\eta ^R`$ for $`0<\eta ^R<\eta _C^R`$. At the point $`\eta ^R=\eta _C^R`$, the derivative $`f_{MF}^{}(\eta ^R)`$ takes the value $`\mathrm{}`$. It then follows from eq.(184) that

$$f_{MF}(\eta _C^R)=\frac{1}{3}.$$

At the point $`\eta _C^R`$ the series expansion for $`f_{MF}(\eta ^R)`$ in powers of $`\eta ^R`$ diverges. Both, from the ratio test on its coefficients and from the Runge-Kutta solution, we find that

$$\eta _C^R=2.517551\mathrm{}$$
(192)
From eq.(184) we find that $`f_{MF}(\eta ^R)\frac{1}{3}`$ has a square root behaviour around $`\eta ^R=\eta _C^R`$:

$$f_{MF}(\eta ^R)\stackrel{\eta ^R\eta _C^R}{=}\frac{1}{3}+\sqrt{\frac{2(\eta _C^R2)}{9\eta _C^R}}\sqrt{\eta _C^R\eta ^R}+\frac{2(\eta _C^R1)}{7\eta _C^R}(\eta _C^R\eta ^R)+𝒪\left[(\eta _C^R\eta ^R)^{3/2}\right]$$

Inserting the numerical value (192) for $`\eta _C^R`$ yields,

$$f_{MF}(\eta ^R)\stackrel{\eta ^R\eta _C^R}{=}\frac{1}{3}+0.213738\mathrm{}\sqrt{\eta _C^R\eta ^R}+0.172225\mathrm{}(\eta _C^R\eta ^R)+𝒪\left[(\eta _C^R\eta ^R)^{3/2}\right]$$
(193)
We see that $`f_{MF}(\eta ^R)`$ becomes complex for $`\eta ^R>\eta _C^R`$. Recall that in the Monte Carlo simulations the gas phase collapses at the point $`\eta _T^R<\eta _C^R`$.

From eqs.(186), we plot $`pV/[NT],S/N`$ and $`\frac{FF_0}{NT}`$ as a function of $`\eta ^R`$ in figs. 9, 10 and 11, respectively.

The points GC, C and MC correspond to the collapse phase transition in the grand canonical, canonical and microcanonical ensembles, respectively. Their positions are determined by the breakdown of the mean field approximation through the analysis of the small fluctuations \[see paper II\].

$`f_{MF}(\eta ^R)`$ is a multivalued function of $`\eta ^R`$ as well as all physical magnitudes \[see eq.(186)\].

As noticed before, the CE only describes the region between the ideal gas point, $`\eta ^R=0`$ and $`C`$ in fig. 1. The MCE goes beyond the point $`C`$ (till the point $`MC`$) with the physical magnitudes described by the second sheet of the square root in eqs.(193) (minus sign). We have near $`C`$ between $`C`$ and $`MC`$,

$`f_{MF}(\eta ^R)\stackrel{\eta ^R\eta _C^R}{=}{\displaystyle \frac{1}{3}}0.213738\mathrm{}\sqrt{\eta _C^R\eta ^R}+0.172225\mathrm{}(\eta _C^R\eta ^R)+𝒪\left[(\eta _C^R\eta ^R)^{3/2}\right]`$
The function $`f_{MF}(\eta ^R)`$ takes its absolute minimum at $`\eta ^R=\eta _{min}^R=2.20731\mathrm{}`$ in the second sheet where $`f_{MF}(\eta _{min}^R)=0.264230\mathrm{}`$.

Since $`f_{MF}(\eta ^R)<\frac{1}{2}`$ implies that the total energy is negative \[see eq.(186)\], the gas is in a ‘bounded state’ for $`\eta ^R`$ beyond $`\eta _2^R=2.18348\mathrm{}`$ in the first sheet.

Since $`\chi (\lambda )`$ and $`\eta (\lambda )`$ are single-valued functions of $`\lambda ,f_{MF}(\eta ^R(\lambda ))`$ defined by eq.(182) is also a single-valued function of $`\lambda `$. That is, $`\lambda `$ is the uniformization  variable. All physical magnitudes are single-valued functions of $`\lambda `$. On the other hand, $`\lambda `$ is an infinite-valued function of $`\eta ^R`$ as one sees from fig. 7 and eq.(176). That is, $`f_{MF}(\eta ^R)`$ has an infinite number of Riemann sheets. However, only the first two sheets are physically realized. The rest are unphysical. A plot of $`f_{MF}(\eta ^R)`$ including all sheets produces a nice spiral converging towards $`\eta ^R=2,f_{MF}(\eta ^R)=1/3`$ for $`\lambda =\mathrm{}`$ as follows from eqs. (175), (176) and (182).

$`\lambda `$ induces a scale transformation in coordinate space as we see in eq.(170) whereas $`\eta ^R`$ plays the coupling constant \[Recall that $`\eta ^R`$ is proportional to Newton’s gravitational constant\].

The variation of $`\eta ^R`$ with respect to $`\lambda `$ yields the renormalization group equation

$$\lambda \frac{d\eta ^R}{d\lambda }=\eta ^R[3f_{MF}(\eta ^R)1]$$

where we used eqs.(171), (173) and (182). Here $`\eta ^R[3f_{MF}(\eta ^R)1]`$ plays the role of the renormalization group beta function. We see that it has two fixed points at $`\eta ^R=0`$ and at $`\eta ^R=\eta _C^R`$. \[See fig. 7 where the running of $`\eta ^R`$ with $`\lambda `$ is exhibited\].

We find from eqs.(174) and (193) near these fixed points

$$\eta ^R\stackrel{\lambda 0}{=}\frac{\lambda ^2}{3}$$

$$\eta ^R\stackrel{\lambda \lambda _C}{=}\eta _C^R\frac{\eta _C^R(\eta _C^R2)}{2\lambda _C^2}(\lambda \lambda _C)^2$$

where the coefficient has the numerical value $`\frac{\eta _C^R(\eta _C^R2)}{2\lambda _C^2}=0.0085515\mathrm{}`$.

### 6.7 Canonical vs. Grand Canonical Ensembles in the Mean Field Approximation

We have seen that the stationary point equations and their respective solutions are closely related in the canonical and grand canonical ensembles \[eqs.(129)-(134) and (146)-(147)\].

Let us now show that physical quantities obtained from both ensembles do coincide in the mean field approximation.

From eq.(137) and (147) we find that

$$K(\eta )=\varphi (\stackrel{}{x})e^{\varphi (\stackrel{}{x})}d^3x+\mathrm{log}C(\eta ).$$
(194)
\[Recall that $`e^{\varphi (\stackrel{}{x})}d^3x=1`$\].

In the spherically symmetric case this integral takes the form

$$4\pi _0^1r^2𝑑r\varphi (r)e^{\varphi (r)}=\varphi (1)+\frac{1}{\eta ^R}_0^1r^2𝑑r\left(\frac{d\varphi }{dr}\right)^2=6[1f_{MF}(\eta ^R)]\eta ^R+\mathrm{log}\left[\frac{3f_{MF}(\eta ^R)}{4\pi }\right]$$
(195)
where we integrated by parts and used eqs.(170), (173) and appendix B.

From eqs.(194) and (195) we find

$$K(\eta ^R)\mathrm{log}C(\eta ^R)=6[1f_{MF}(\eta ^R)]\eta ^R+\mathrm{log}\left[\frac{3f_{MF}(\eta ^R)}{4\pi }\right].$$

Inserting this result into the linear differential equations (141) leads to the solution,

$$C(\eta ^R)=\frac{4\pi }{3}\frac{\mathrm{exp}[\eta ^R]}{f_{MF}(\eta ^R)}\text{and}K(\eta ^R)=6[1f_{MF}(\eta ^R)]$$
(196)
We then find from eqs.(147), (178) and (185) that

$$\mathrm{log}C(\eta )=a_s.$$
(197)
Combining eq.(196) with eqs.(138), (135)-(140) and (142) shows that the canonical and the grand canonical ensembles yields identical physical magnitudes (pressure, energy, entropy, free energy, specific heats, compressibilities, speed of sound) and the same equation of state in the mean field approximation.

The thermodynamical potential,

$$\mathrm{\Omega }T\mathrm{log}𝒵_{GC}=N[3f_{MF}(\eta ^R)2]$$

is not equal to $`PV`$. That is, here $`\mathrm{\Omega }PV`$ and we have instead

$$\mathrm{\Omega }+PV=2NT[1f_{MF}(\eta ^R)]$$

This relation is analogous to eq.(76). $`\mathrm{\Omega }`$ differs here from $`PV`$ since for the self-gravitating gas we have $`NL`$ instead of the usual relation $`NL^3`$.

## 7 Specific Heats, Speed of Sound and Compressibility

The specific heat at constant volume in the mean field approximation takes the form

$$(c_V)_{MF}=6f_{MF}(\eta ^R)\frac{7}{2}+\eta ^R+\frac{\eta ^R2}{3f_{MF}(\eta ^R)1}$$
(198)
where we used eqs.(77) and (184).

We plot in Fig. 12 eq.(198) for $`(c_V)_{MF}`$ as a function of $`\eta `$. We see that $`(c_V)_{MF}`$ increases with $`\eta `$ till it tends to $`+\mathrm{}`$ for $`\eta ^R\eta _C^R`$. It has a square-root branch point at the point C. In the stretch C-MC (only physically realized in the microcanonical ensemble), $`(c_V)_{MF}`$ becomes negative. We shall not discuss here the peculiar properties of systems with negative $`C_V`$ as they can be find in refs.

From eqs.(193) and (198) we obtain the following behaviour near the point $`C`$ in the positive (first) branch

$`(c_V)_{MF}`$ $`\stackrel{\eta ^R\eta _C^R}{=}`$ $`0.80714\mathrm{}(\eta _C^R\eta ^R)^{1/2}0.19924\mathrm{}+𝒪(\sqrt{\eta _C^R\eta ^R})`$ (199)
and between $`C`$ and $`MC`$ in the negative (second) branch

$`(c_V)_{MF}`$ $`\stackrel{\eta ^R\eta _C^R}{=}`$ $`0.80714\mathrm{}(\eta _C^R\eta ^R)^{1/2}0.19924\mathrm{}+𝒪(\sqrt{\eta _C^R\eta ^R})`$
Finally, $`(c_V)_{MF}`$ vanishes at the point MC $`\eta _{MC}^R=2.03085\mathrm{}`$.

The isothermal compressibility in mean field follows from eqs.(84) and (184)

$$(\kappa _T)_{MF}=\frac{3}{2f_{MF}(\eta ^R)}\left[1+\frac{\eta ^R2}{6f_{MF}(\eta ^R)\eta ^R}\right].$$
(200)
We plot $`(\kappa _T)_{MF}`$ in fig. 13. We see that $`(\kappa _T)_{MF}`$ is positive for $`0\eta ^R<\eta _0^R=2.43450\mathrm{}`$ where $`(\kappa _T)_{MF}`$ diverges. The point $`\eta _0^R`$ is defined by the equation

$$6f_{MF}(\eta _0^R)\eta _0^R=0.$$
(201)
We find from eqs.(184) and (201) that

$$f_{MF}^{}(\eta _0^R)=\frac{1}{2}$$
(202)
$`(\kappa _T)_{MF}`$ diverges for $`\eta ^R\eta _0^R`$ as

$$(\kappa _T)_{MF}\stackrel{\eta ^R\eta _0^R}{=}\frac{9(\eta _0^R2)}{4\eta _0^R(\eta _0^R\eta ^R)}+𝒪(1)=\frac{0.40157\mathrm{}}{\eta _0^R\eta ^R}+𝒪(1).$$

$`(\kappa _T)_{MF}`$ is negative for $`\eta _0^R<\eta ^R<\eta _C^R`$ and exactly vanishes at the point $`C`$. $`(\kappa _T)_{MF}`$ then becomes positive in the stretch between $`C`$ and $`MC`$ only physically realized in the microcanonical ensemble.

Notice that the singularity of $`(\kappa _T)_{MF}`$ at $`\eta ^R=\eta _0^R=2.43450\mathrm{}`$ is before but near the point $`C`$. It appears as a preliminary signal of the phase transition at $`C`$. $`\eta _0^R`$ is probably the transition point $`\eta _T`$ seen with the Monte Carlo simulations (see fig. 1). (Recall that $`\eta _T1.515`$ corresponds to $`\eta _T^R2.44`$).

It is easy to understand the meaning of a large compressibility. From the definition (84)

$$\frac{\delta V}{V}=K_T\delta p=\kappa _T\frac{V\delta p}{NT}.$$
(203)
A large compressibility implies that a small increase in the pressure ($`\delta pNT/V`$) produces a large change in the density of the gas. That means a very soft fluid.

For negative compressibility, eq.(203) tells us that the gas increases its volume when the external pressure on it increases. This is clearly an unusual behaviour that leads to instabilities as we shall see below.

The specific heat at constant pressure in the mean field approximation takes the form

$$(c_P)_{MF}=12f_{MF}(\eta ^R)\frac{3}{2}+\frac{24\left(\eta ^R2\right)f_{MF}(\eta ^R)}{6f_{MF}(\eta ^R)\eta ^R}$$
(204)
where we used eqs.(81) and (184). We plot $`(c_P)_{MF}`$ in fig. 14. We see that $`(c_P)_{MF}`$ is positive and grows with $`\eta ^R`$ till it diverges at the same point where $`(\kappa _T)_{MF}`$ diverges $`\eta ^R=\eta _0^R=2.43450\mathrm{}`$. That is,

$$(c_P)_{MF}\stackrel{\eta ^R\eta _0^R}{=}\frac{\eta _0^R(\eta _0^R2)}{\eta _0^R\eta ^R}+𝒪(1)=\frac{1.05779\mathrm{}}{\eta _0^R\eta ^R}+𝒪(1).$$

$`(c_P)_{MF}`$ becomes negative for $`\eta _0^R<\eta ^R<\eta _C^R`$. It keeps negative in the C-MC section till the point $`\eta ^R=\eta _1^R=2.14675\mathrm{}`$ where it becomes positive. The point $`\eta _1^R`$ is defined by the equation

$$24f_{MF}^2(\eta _1^R)+\left(4\eta _1^R19\right)f_{MF}(\eta _1^R)+\frac{\eta _1^R}{2}=0.$$
(205)
The speed of sound squared at the surface in the mean field approximation takes the form

$$\frac{v_s^2}{T}=\frac{f_{MF}(\eta ^R)}{3}\left[4+\frac{3f_{MF}(\eta ^R)+\frac{\eta ^R}{2}2}{6f_{MF}^2(\eta ^R)+\left(\eta ^R\frac{11}{2}\right)f_{MF}(\eta ^R)+\frac{1}{2}}\right],$$
(206)
where we used eqs.(89) and (184). We plot $`\frac{v_s^2}{T}`$ as a function of $`\eta ^R`$ in fig. 15. We see that $`\frac{v_s^2}{T}(\eta ^R)`$ is positive and decreasing with $`\eta ^R`$ in the whole interval between $`\eta ^R=0`$ and $`C`$. At the point $`C`$ it takes the value $`\frac{v_s^2}{T}(\eta _C^R)=11/18`$. Then, $`\frac{v_s^2}{T}(\eta ^R)`$ decreases between $`C`$ and $`MC`$ becoming negative at $`\eta _1^R=2.14675\mathrm{}`$ in the second sheet where it vanishes. Notice that $`\frac{v_s^2}{T}(\eta ^R)`$ and $`(c_P)_{MF}`$ vanish at the same point $`\eta _1^R`$ defined by eq.(205).

$`v_s^2<0`$ indicates an instability. That is, small density fluctuations grow exponentially in time instead of propagating harmonically. It is remarkable that $`v_s^2`$ becomes negative at $`\eta _1^R=2.14675\mathrm{}`$ in the second sheet before but near the $`MC`$ critical point $`\eta _{MC}^R=2.03085\mathrm{}`$ in the second sheet. Somehow, the change of sign in $`v_s^2`$ announces the $`MC`$ critical point.

$`\frac{v_s^2}{T}(\eta ^R)`$ tends to $`\mathrm{}`$ for $`\eta ^R\eta _{MC}^R`$. Notice that the denominator in eq.(206) exactly vanishes at $`\eta ^R=\eta _{MC}^R`$ \[see Table 1\].

The adiabatic compressibility $`\kappa _S`$ is not here an independent quantity. We find from eqs.(87), (198), (200) and (204),

$$\kappa _S=\frac{c_V}{c_P}\kappa _T=\frac{3}{f_{MF}(\eta ^R)}\frac{12f_{MF}^2(\eta ^R)+(2\eta ^R11)f_{MF}(\eta ^R)+1}{48f_{MF}^2(\eta ^R)+(8\eta ^R38)f_{MF}(\eta ^R)+\eta ^R}$$

That is,

$$\kappa _S=\frac{T}{v_s^2}.$$

|  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- |
| POINT | $`\lambda `$ | $`\eta ^R`$ | Defining Equation | $`f_{MF}(\eta ^R)`$ | PHYSICAL MEANING |
|  |  |  |  |  |  |
|  |  |  |  |  |  |
| GC | $`1.7772\mathrm{}`$ | $`0.797375\mathrm{}`$ | $`23\eta _{GC}^Rf_{MF}(\eta _{GC}^R)=0`$ | $`0.836076\mathrm{}`$ | Collapse in the GCE. |
|  |  |  |  |  |  |
|  |  |  |  |  | Energy density |
| $`\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}3}`$ | $`3.38626\mathrm{}`$ | $`1.73745\mathrm{}`$ | $`3\eta ^R+\chi (\lambda )=0`$ | $`0.622424\mathrm{}`$ | vanishes at $`r=0`$. |
|  |  |  |  |  |  |
|  |  |  |  |  |  |
| $`\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}2}`$ | $`4.73739\mathrm{}`$ | $`2.18348\mathrm{}`$ | $`2f_{MF}(\eta _2^R)1=0`$ | $`0.5`$ | Total Energy vanishes. |
|  |  |  |  |  |  |
|  |  |  |  |  |  |
| $`\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}`$ | $`6.45077\mathrm{}`$ | $`2.43450\mathrm{}`$ | $`6f_{MF}(\eta _0^R)\eta _0^R=0`$ | $`0.40575\mathrm{}`$ | $`\kappa _T`$ and $`c_P`$ diverge. |
|  |  |  |  |  |  |
|  |  |  |  |  | Collapse in the CE. |
| C | $`8.993195\mathrm{}`$ | $`2.517551\mathrm{}`$ | $`13f_{MF}(\eta _C^R)=0`$ | $`1/3`$ | $`c_V`$ diverges. |
|  |  |  |  |  |  |
|  |  |  |  |  | Minimum of |
| Min | $`22.5442\mathrm{}`$ | $`2.20731\mathrm{}`$ | $`f_{MF}^{}(\eta _{min}^R)=0`$ | $`0.264230\mathrm{}`$ | pV/\[NT\] |
|  |  |  |  |  | in the gas phase |
|  |  |  |  |  |  |
| $`\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}`$ | $`25.7991\mathrm{}`$ | $`2.14675\mathrm{}`$ | $`\mathrm{\hspace{0.33em}\hspace{0.33em}48}f_{MF}^2(\eta _1^R)(388\eta _1^R)f_{MF}(\eta _1^R)+`$ | $`0.265290\mathrm{}`$ | $`v_s^2`$ and $`c_P`$ vanish. |
|  |  |  | $`+\eta _1^R=0`$ |  |  |
|  |  |  |  |  | Collapse in the MCE. |
| MC | $`34.36361\mathrm{}`$ | $`2.03085\mathrm{}`$ | $`12f_{MF}^2(\eta _{MC}^R)(112\eta _{MC}^R)f_{MF}(\eta _{MC}^R)+`$ | $`0.273512\mathrm{}`$ | $`c_V`$ vanishes. |
|  |  |  | $`+1=0`$ |  |  |

TABLE 1. Values of the critical points in the three ensembles GC, C and MC (using mean field) and further characteristic points for spherical symmetry. $`pV/[NT],E`$ and $`S\mathrm{}`$ for $`\eta ^R\eta _{GC}^R`$ and $`\eta ^R\eta _C^R`$. Notice that $`\eta _{Min},\eta _1`$ and $`\eta _{MC}`$ are in the second Riemann sheet.

## 8 Discussion

We have presented here a set of new results for the self-gravitating thermal gas obtained by Monte Carlo and analytic methods. They provide a complete picture for the thermal self-gravitating gas.

Contrary to the usual hydrostatic treatments , we do not assume here an equation of state but we obtain the equation of state from the partition function \[see eq.(64)\]. We find at the same time that the relevant variable is here $`\eta ^R=Gm^2N/[V^{1/3}T]`$. The relevance of the ratio $`Gm^2/[V^{1/3}T]`$ has been noticed on dimensionality grounds . However, dimensionality arguments alone cannot single out the crucial factor $`N`$ in the variable $`\eta ^R`$.

The crucial point is that the thermodynamic limit exist if we let $`N\mathrm{}`$ and $`V\mathrm{}`$ keeping $`\eta ^R`$ fixed. Notice that $`\eta `$ contains the ratio $`NV^{1/3}`$ and not $`N/V`$. This means that in this thermodynamic limit $`V`$ grows as $`N^3`$ and thus the volume density $`\rho =N/V`$ decreases as $`N^2`$. $`\eta `$ is to be kept fixed for a thermodynamic limit to exist in the same way as the temperature. $`pV`$, the energy $`E`$, the free energy, the entropy are functions of $`\eta `$ and $`T`$ times $`N`$. The chemical potential, specific heat, etc. are just functions of $`\eta `$ and $`T`$.

We find collapse phase transitions both in the canonical and in the microcanonical ensembles. They take place at different values of the thermodynamic variables and are of different nature. In the CE the pressure becomes large and negative in the collapsed phase. The phase transition in the MCE is sometimes called ‘gravothermal catastrophe’. We find that the temperature and pressure increase discontinuously at the MCE transition. Both are zeroth order phase transitions (the Gibbs free energy is discontinuous). The two phases cannot coexist in equilibrium since the pressure has different values at each phase.

The parameter $`\eta ^R`$ \[introduced in eq.(53)\] can be related to the Jeans length of the system

$$d_J=\sqrt{\frac{3T}{m}}\frac{1}{\sqrt{Gm\rho }},$$
(207)
where $`\rho N/V`$ stands for the number volume density. Combining eqs.(53) and (207) yields

$$\eta ^R=3\left(\frac{L}{d_J}\right)^2.$$

We see that the phase transition in the canonical ensemble takes place for $`d_JL`$. \[The precise numerical value of the proportionality coefficient depends on the geometry\]. For $`d_J>L`$ we find the gaseous phase and for $`d_J<L`$ the system condenses as expected. Hence, the collapse phase transition in the canonical ensemble is related to the Jeans instability.

The latent heat of the transition ($`q`$) is negative in the CE transition indicating that the gas releases heat when it collapses \[see eq.(113)\]. The MCE transition exhibits an opposite behaviour. The Gibbs free energy increases at the MCE collapse phase transition (point MC in fig.1) whereas it decreases at the CE transition \[point T in fig. 1, see eq.(112)\]. Also, the average distance between particles increases at the MCE phase transition whereas it decreases dramatically in the CE phase transition. These differences are related to the MCE constraint keeping the energy fixed whereas in the CE the system exchanges energy with an external heat bath keeping fixed its temperature. The constant energy constraint in the MCE keeps the gas stable in a wider domain and makes the collapse transition softer than in the CE. Notice that the core is much tighter and the halo much smaller in the CE than in the MCE \[see figs. 4 and 6\].

## 9 Acknowledgements

One of us (H J de V) thanks M. Picco for useful discussions on Monte Carlo methods. We thank S. Bouquet for useful discussions and J. Katz for calling our attention on ref..

## Appendix A Functional integration Measure in the Mean Field Approach

We follow the derivation of ref. for the functional integral measure. We want to recast

$$e^{\mathrm{\Phi }_N(\eta )}=_0^1\mathrm{}_0^1\underset{l=1}{\overset{N}{}}d^3r_le^{\eta u(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_N)},$$
(208)
as a functional integral in the large $`N`$ limit.

We start by dividing the domain of integration (of unit volume) into $`M`$ cells. Each cell is of volume $`\omega _r`$ and contains $`k_r`$ particles with $`1rM`$. Therefore,

$$\underset{r=1}{\overset{M}{}}k_r=N,\underset{r=1}{\overset{M}{}}\omega _r=1.$$

We can thus rewrite the multiple integral (208) as follows:

$$e^{\mathrm{\Phi }_N(\eta )}=\underset{k_1,\mathrm{},k_M}{}\delta \left(N\underset{r=1}{\overset{M}{}}k_r\right)\frac{N!}{_{r=1}^Mk_r!}\underset{r=1}{\overset{M}{}}(\omega _r)^{k_r}e^J$$

where

$$J=\frac{1}{2}\underset{r,r^{}}{}k_rk_r^{}V_{r,r^{}}+\frac{1}{2}\underset{r}{}k_rV_{r,r}+\frac{1}{2}\underset{r,r^{},r^{\prime \prime }}{}k_rk_r^{}k_{r^{\prime \prime }}\left[<V_{r,r^{}}V_{r,r^{\prime \prime }}><V_{r,r^{}}><V_{r,r^{\prime \prime }}>\right]+\mathrm{}$$

and

$$V_{r,r^{}}=\frac{1}{\omega _r\omega _r^{}}\frac{\eta }{N}_0^1_0^1\frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|}|_{\stackrel{}{r}_1\omega _r,\stackrel{}{r}_2\omega _r^{}}.$$

Assuming $`1/N\omega _r<N^{2/3}`$ one can neglect in $`J`$ terms quadratic and higher in $`V_{r,r^{}}`$.

The particle density is defined as

$$N\rho (\stackrel{}{r})\underset{r=1}{\overset{M}{}}\frac{k_r}{\omega _r}\theta (\stackrel{}{r}\omega _r).$$

Therefore, we can write the sums over $`r`$ as integrals in the following way

$$\frac{1}{2}\underset{r,r^{}}{}k_rk_r^{}V_{r,r^{}}=\frac{\eta }{2N}_0^1_0^1\frac{d^3r_1d^3r_2}{|\stackrel{}{r}_1\stackrel{}{r}_2|}\rho (\stackrel{}{r}_1)\rho (\stackrel{}{r}_2).$$

Using Stirling’s’ formula one finds that

$$\underset{r=1}{\overset{M}{}}\frac{(\omega _r)^{k_r}}{k_r!}\stackrel{N\mathrm{}}{=}\frac{1}{N^N}\underset{r=1}{\overset{M}{}}\frac{1}{\sqrt{2\pi k_r}}e^{N{\scriptscriptstyle d^3x\rho (\stackrel{}{x})\mathrm{log}[\rho (\stackrel{}{x})/e]}}.$$

Collecting all terms yields,

$$N!\underset{r=1}{\overset{M}{}}\frac{(\omega _r)^{k_r}}{k_r!}e^J\stackrel{N\mathrm{}}{=}e^{\frac{N\eta }{2}{\scriptscriptstyle {\scriptscriptstyle \frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}}\rho (\stackrel{}{x})\rho (\stackrel{}{y})}N{\scriptscriptstyle d^3x\rho (\stackrel{}{x})\mathrm{log}[\rho (\stackrel{}{x})/e]}}$$

whereas the constraint in the number of particles takes the form

$$\delta \left(N\underset{r=1}{\overset{M}{}}k_r\right)=\frac{1}{N}\delta \left(d^3x\rho (\stackrel{}{x})1\right)$$

and finally,

$$e^{\mathrm{\Phi }_N(\eta )}\stackrel{N>>1}{=}\frac{1}{N}D\rho e^{\frac{N\eta }{2}{\scriptscriptstyle {\scriptscriptstyle \frac{d^3xd^3y}{|\stackrel{}{x}\stackrel{}{y}|}}\rho (\stackrel{}{x})\rho (\stackrel{}{y})}N{\scriptscriptstyle d^3x\rho (\stackrel{}{x})\mathrm{log}[\rho (\stackrel{}{x})/e]}}\delta \left(d^3x\rho (\stackrel{}{x})1\right)$$

Replacing the Dirac delta by its Fourier representation

$$\frac{1}{N}\delta \left(d^3x\rho (\stackrel{}{x})1\right)=\frac{d\widehat{a}}{2\pi }e^{iN\widehat{a}\left({\scriptscriptstyle d^3x\rho (\stackrel{}{x})}1\right)}$$

yields eq.(117).

## Appendix B Calculation of the saddle point

We prove in this Appendix that the integral

$$I(\lambda )_0^\lambda x^2𝑑x[\chi ^{}(x)]^2$$
(209)
takes the value

$$I(\lambda )=\lambda \eta ^R(6\eta ^R)2\lambda ^3e^{\chi (\lambda )}$$
(210)
Here $`\chi (x)`$ is a regular solution of eq.(171) in the interval $`0x\lambda `$ fulfilling the relation (173).

We start by computing the derivative of $`I(\lambda )`$ in two ways. According to the definition (209)

$$\frac{dI(\lambda )}{d\lambda }=\lambda ^2[\chi ^{}(\lambda )]^2$$

Then, we compute the derivative of eq.(210) with respect to $`\lambda `$ and use eqs.(171) and (173). We find after calculation that both results coincide.

Finally, we observe that both eqs.(209) and (210) vanish at $`\lambda =0`$. Therefore, eq.(210) is valid.

## Appendix C Abel’s equation of first kind for the equation of state

In the mean field approximation the equation of state for spherical symmetry satisfies the first order differential equation (184)

$$\eta ^R(3f_{MF}1)f_{MF}^{}(\eta ^R)+(3f_{MF}3+\eta ^R)f_{MF}=0.$$
(211)
with the boundary condition $`f_{MF}(0)=1`$.

We can solve eq.(211) in power series in $`\eta ^R`$ around the origin,

$$f_{MF}(\eta )=1+\underset{n=1}{\overset{\mathrm{}}{}}f_n\eta ^n$$
(212)
Inserting eq.(212) into eq.(211) yields the quadratic recurrence relation

$$f_n=\frac{1}{2n+3}\left[f_{n1}+3\underset{k=2}{\overset{n}{}}kf_{k1}f_{nk+1}\right]\text{for}n2.$$

where $`f_1=\frac{1}{5}`$.

We find from this recurrence relation,

$$f_2=\frac{1}{175},f_3=\frac{2}{1575},f_4=\frac{991}{3031875}$$

All coefficients $`f_n`$ are negative rational numbers for $`n1`$. They decrease very fast with $`n`$ as

$$f_n\stackrel{n1}{=}\frac{0.0956678\mathrm{}}{[\eta _C^R]^nn^{3/2}}\left[1+𝒪\left(\frac{1}{n}\right)\right]$$

This formula reproduces the large orders of the expansion of $`\sqrt{\eta _C^R\eta ^R}`$ describing the behaviour of $`f_{MF}(\eta )`$ near $`\eta _C^R`$ \[see eq.(193) and ref.\]

$$\sqrt{\eta _C^R\eta ^R}=\frac{1}{2}\sqrt{\frac{\eta _C^R}{\pi }}\underset{n=0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(n\frac{1}{2})}{n!}\left(\frac{\eta ^R}{\eta _C^R}\right)^n$$

Notice that

$$\frac{1}{2}\sqrt{\frac{\eta _C^R}{\pi }}\frac{\mathrm{\Gamma }(n\frac{1}{2})}{n!}\stackrel{n1}{=}\frac{0.447594\mathrm{}}{n^{3/2}}\left[1+𝒪\left(\frac{1}{n}\right)\right]$$

and that $`0.213738\mathrm{}\times 0.447594\mathrm{}=0.0956678\mathrm{}`$.

The power series (212) thus has a radius of convergence $`\eta _C^R=2.517551\mathrm{}`$. The singularity of $`f_{MF}(\eta )`$ nearest to the origin is thus the critical point.