Validity conditions of the hydrostatic approach for self-gravitating systems: a microcanonical analysis

We consider a classical gravitational system $\mathcal{S}$, made up of N identical hard spheres with mass m and diameter σ, enclosed in a spherical box of volume Λ = 4πR³/3. The corresponding Hamiltonian reads

$$\begin{equation} H_{N} = \sum _{i=1}^N \frac{\mathbf {p}_i^2}{ 2m} + \frac{1}{2}\sum _{i \ne j} v(|\mathbf {r}_i - \mathbf {r}_j|) \end{equation} \tag{ 1 }$$

with the two-body interaction potential

$$\begin{equation} v(r)=\infty \;\;\; {\rm for} \;\;\; r<\sigma , \;\; v(r)=-Gm^2/r \;\;\; {\rm for} \;\;\; r>\sigma . \end{equation} \tag{ 2 }$$

We consider that the previous system is isolated and does not exchange energy with some thermostat. Thus, its energy is fixed and equal to some value E. We assume that the corresponding stationary state is described within the microcanonical ensemble. The corresponding distribution of positions and momenta of the N particles in the canonical phase space reads

$$\begin{equation} f_{\rm micro}(\mathbf {r}_1,\ldots ,\mathbf {r}_N,\mathbf {p}_1,\ldots ,\mathbf {p}_N)=A_N \; \delta (E-H_N) , \end{equation} \tag{ 3 }$$

while the total number of microstates is

$$\begin{equation} \Omega (E,N,\Lambda ) = A_N \int _{\Lambda ^N \times R^{3N}} \prod _i \mathrm{d}^3 \mathbf {r}_i \,\mathrm{d}^3 \mathbf {p}_i \,\delta (E-H_N) , \end{equation} \tag{ 4 }$$

where A_N is some normalization constant, which is not relevant for our purpose.

The microcanonical distribution f_micro remains stationary under the time evolution generated by Hamiltonian H_N. We stress that Ω(E, N, Λ) is finite for σ > 0, while it diverges for point particles with σ = 0 for N ⩾ 3 as shown in [12]. The corresponding stationary state of total mass M = Nm is controlled by three independent dimensionless parameters, namely N, ε (which is the energy in units of GM²/R), and the packing fraction η = πnσ³/6 with the particle density n = N/Λ. Notice that ε is identical to the parameter which entirely controls the stationary state of a system composed of point particles within Antonov's mean-field theory [2].

Since we are interested in the properties of a system with a large number of particles, it is useful to consider a sequence of larger and larger auxiliary hard-sphere gravitating systems $\mathcal{S}_{\rm a}$ with increasing particle numbers N_a and take some limit where N_a → ∞. For ordinary systems with short-range interactions, that limit is nothing but the usual thermodynamic limit where both the energy per particle E_a/N_a and the particle density n_a = N_a/Λ_a are kept fixed. In that limit, the physical parameters which describe a given particle are also kept fixed. Because of both the attractive and the long-range natures of the gravitational interaction, such a limit would provide a collapsed state with non-extensive properties. In order to describe other physical situations of interest, other limits have been introduced in the literature, like the one describing an infinitely diluted system [14].

Here, we want to build a SL when N_a → ∞, which describes an infinite continuous fluid $\mathcal{S}_{\infty }$ with the usual extensive properties and depending only on the parameters ε and η characterizing the actual system of interest $\mathcal{S}$. For that purpose, we consider that the parameters which define the particles of $\mathcal{S}_{\rm a}$ vary with N_a like power laws, while the gravitational constant G is not rescaled, in contrast to the case for some mean-field scalings introduced in the literature [7, 8]. Therefore, we set σ_a = (N_a/N)^ασ and m_a = (N_a/N)^δm, while the size R_a of the spherical box is chosen to diverge as (N_a/N)^γR with γ > 0 so that the system does indeed become infinitely extended. The particle density n_a = N_a/Λ_a behaves then as $N_{\rm a}^{1-3\gamma }$. By imposing that $\eta _{\rm a}=\pi n_{\rm a}\sigma _{\rm a}^3/6$ remains constant and equal to η, we find a first constraint

$$\begin{equation} 1-3\gamma + 3 \alpha = 0 . \end{equation} \tag{ 5 }$$

Since E_a is chosen extensive, while $\varepsilon _{\rm a}=E_{\rm a}/ (GM_{\rm a}^2/R_{\rm a})$ with M_a = N_am_a is kept fixed, the gravitational energy $GM_{\rm a}^2/R_{\rm a}$ also has to be extensive. This implies the second constraint

$$\begin{equation} 2+ 2\delta -\gamma = 1 . \end{equation} \tag{ 6 }$$

At this stage, we have the two constraints (5) and (6) for three exponents, so we can impose a third constraint. By analogy with the usual TL where the number density is fixed, we impose here that the mass density ρ_a = m_an_a remains constant and equal to ρ, which leads to

$$\begin{equation} 1-3\gamma + \delta = 0 . \end{equation} \tag{ 7 }$$

The three exponents are then readily determined as α = −2/15, δ = −2/5 and γ = 1/5. Thus, the power laws defining the required SL when N_a → ∞ are

$$\begin{eqnarray} &&\nonumber R_{\rm a}=(N_{\rm a}/N)^{1/5} R , \quad\; m_{\rm a}= (N_{\rm a}/N)^{-2/5} m , \;\; \\ &&\sigma _{\rm a}= (N_{\rm a}/N)^{-2/15} \sigma , \;\; E_{\rm a}=(N_{\rm a}/N) E . \end{eqnarray} \tag{ 8 }$$

As illustrated in figure 1, that limit clearly describes an infinite continuous medium $\mathcal{S}_{\infty }$, the state of which is controlled by the two independent parameters ε and η. Note that, within the scaling considered, the mean free path $\ell _{\rm a}=1/(\pi \sqrt{2} \; \sigma _{\rm a}^2 n_{\rm a})$ remains proportional to the diameter σ_a of the hard spheres, $\ell _{\rm a}=\sigma _{\rm a}/(6 \sqrt{2} \; \eta )$. Therefore, for high dilutions such that η ≪ 1, ℓ_a is much larger than σ_a.

Zoom In Zoom Out Reset image size

For any allowed configuration, the potential energy of the system $\mathcal{S}_{\rm a}$,

$$\begin{equation} V_{N_{\rm a}} = -\frac{1}{ 2}\sum _{i \ne j} \frac{Gm_{\rm a}^2 }{ |\mathbf {r}_i - \mathbf {r}_j|}, \end{equation} \tag{ 9 }$$

is larger than that of the collapsed configuration where the N_a hard spheres make a single cluster with size $L_{\rm coll} \sim N_{\rm a}^{1/3} \sigma _{\rm a}$, which is of order $-Gm_{\rm a}^2N_{\rm a}^2/L_{\rm coll}$. In the SL, the collapse radius L_coll diverges as $N_{\rm a}^{1/5}$ and this provides the classical version of H-stability [13]:

$$\begin{equation} V_{N_{\rm a}} \ge -C_{\rm HS} \frac{G N^{2/3} m^2 }{ \sigma } N_{\rm a} , \end{equation} \tag{ 10 }$$

where C_HS is a positive real number entirely determined by the geometrical arrangement of the hard spheres at maximum packing. On the other hand, the potential energy should reach its maximum when all particles are homogeneously distributed on the spherical surface of the box. That maximum behaves as the gravitational energy of an empty sphere carrying the surface mass density $N_{\rm a}m_{\rm a}/(4\pi R_{\rm a}^2)$, so we expect the upper bound

$$\begin{equation} V_{N_{\rm a}} \le - \frac{G N m^2}{2 R} N_{\rm a} . \end{equation} \tag{ 11 }$$

According to bounds (10) and (11), $V_{N_{\rm a}}$ remains of order N_a for any configuration. This implies extensive lower and upper bounds for the average potential energy $\langle V_{N_{\rm a}} \rangle$, where 〈A〉 denotes the microcanonical average of any observable A weighted by the microcanonical distribution for system $\mathcal{S}_{\rm a}$, namely expression (3) with (N, A_N, H_N, E) replaced by $N_{\rm a},A_{N_{\rm a}},H_{N_{\rm a}},E_{\rm a}$. Therefore, the non-extensive behaviour of the average potential energy encountered in the usual thermodynamic limit does not occur in the present SL. Accordingly, the extensivity of $\langle V_{N_{\rm a}} \rangle$ should be ensured here. Since the typical potential energy $GM_{\rm a}^2/R_{\rm a}$ is also extensive by construction of the SL, we find

$$\begin{equation} \lim _{\rm SL}\frac{\langle V_{N_{\rm a}} \rangle }{(GM_{\rm a}^2/R_{\rm a})} = u(\varepsilon ,\eta ) , \end{equation} \tag{ 12 }$$

where u(ε, η) is the dimensionless potential energy per particle of the infinite system $\mathcal{S}_{\infty }$, which only depends on the intensive dimensionless parameters ε and η.

For further purposes, it is useful to introduce the n-body mass distributions. Let us consider first the one-body mass distribution, which reads

$$\begin{equation} \rho _{\rm a}(\mathbf {r}) = m_{\rm a}\left\langle \sum _{i=1}^{N_{\rm a}} \delta (\mathbf {r}_i - \mathbf {r}) \right\rangle . \end{equation} \tag{ 13 }$$

According to the scaling properties of the potential energy described above, we can reasonably expect ρ_a(r) to take a well-defined shape in the SL, which is the mass distribution ρ_∞(r) of the infinite system $\mathcal{S}_{\infty }$. Taking into account dimensional and spherical-symmetry properties, this leads to

$$\begin{equation} \lim _{\rm SL}\rho _{\rm a}(qR_{\rm a}) = \rho g(q;\varepsilon ,\eta ) , \end{equation} \tag{ 14 }$$

where q = r/R_a remains fixed in the SL. The function g(q; ε, η) is the dimensionless mass distribution, as a function of the rescaled distance q, for the equilibrium state of $\mathcal{S}_{\infty }$ determined by the intensive parameters (ε, η). The two-body mass distribution is

$$\begin{equation} \rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime }) = m_{\rm a}^2\left\langle \sum _{i \ne j} \delta (\mathbf {r}_i - \mathbf {r})\delta (\mathbf {r}_j - \mathbf {r}^{\prime }) \right\rangle , \end{equation} \tag{ 15 }$$

while similar definitions hold for higher order n-body mass distributions with n ⩾ 3. The scaling property (14) can be extended to such mass distributions, namely

$$\begin{equation} \lim _{\rm SL}\rho _{\rm a}^{(2)}(R_{\rm a}\mathbf {q},R_{\rm a}\mathbf {q}^{\prime }) = \rho ^2 g^{(2)}(\mathbf {q},\mathbf {q}^{\prime };\varepsilon ,\eta ) , \end{equation} \tag{ 16 }$$

and so on.

We can check that the expected extensive behaviour of the average potential energy of $\mathcal{S}_{\rm a}$ is consistent with the above scaling properties for the mass distributions. This is easily seen by starting from the integral expression

$$\begin{equation} \langle V_{N_{\rm a}} \rangle = -\frac{1}{ 2} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}\, \mathrm{d}^3 \mathbf {r}^{\prime } \,\rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime }) \frac{G}{|\mathbf {r}- \mathbf {r}^{\prime }|} . \end{equation} \tag{ 17 }$$

After making the change of variables r → Rq, r' → Rq', and using the scaling behaviour (16) for $\rho _{\rm a}^{(2)}(R\mathbf {q},R\mathbf {q}^{\prime })$, we find that $\langle V_{N_{\rm a}} \rangle$ is indeed an extensive quantity in the SL.

Let us now consider the fluctuations of the potential energy around its average. Such fluctuations can be rewritten, similarly to formula (17), as

$$\begin{eqnarray} &&\nonumber \langle V_{N_{\rm a}}^2 \rangle - [\langle V_{N_{\rm a}} \rangle ]^2 = \frac{1}{2} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}\, \mathrm{d}^3 \mathbf {r}^{\prime } \,\rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime })\frac{G^2 m_{\rm a}^2}{ |\mathbf {r}- \mathbf {r}^{\prime }|^2} \\ &&\nonumber + \int _{\Lambda _{\rm a}^3} \mathrm{d}^3 \mathbf {r}\,\mathrm{d}^3 \mathbf {r}^{\prime }\, \mathrm{d}^3 \mathbf {r}^{\prime \prime }\, \rho _{\rm a}^{(3)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime }) \frac{G^2 m_{\rm a}}{ |\mathbf {r}- \mathbf {r}^{\prime }||\mathbf {r}- \mathbf {r}^{\prime \prime }|} \\ &&+ \frac{1}{4} \int _{\Lambda _{\rm a}^4} \!\!\!\!\mathrm{d}^3 \mathbf {r}\,\mathrm{d}^3 \mathbf {r}^{\prime } \,\mathrm{d}^3 \mathbf {r}^{\prime \prime }\, \mathrm{d}^3 \mathbf {r}^{\prime \prime \prime }\, \big[\rho _{\rm a}^{(4)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime }) - \rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime })\rho _{\rm a}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })\big]\nonumber\\ &&\times \frac{G^2}{ |\mathbf {r}- \mathbf {r}^{\prime }||\mathbf {r}^{\prime \prime } - \mathbf {r}^{\prime \prime \prime }|} , \end{eqnarray} \tag{ 18 }$$

where $\rho _{\rm a}^{(3)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime })$ and $\rho _{\rm a}^{(4)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })$ are the three- and four-body mass distributions. In the SL, the first two terms in expression (18) can be estimated by replacing $\rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime })$ and $\rho _{\rm a}^{(3)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime })$ by their scaling behaviours in terms of g⁽²⁾(q, q'; ε, η) and g⁽³⁾(q, q', q''; ε, η). The first term is then found to be of order $N_{\rm a}^0$, while the second one behaves as N_a.

If we estimate the third term in expression (18) by using only the scaling behaviours of the mass distributions involved, we obtain not surprisingly an $N_{\rm a}^2$-behaviour. However such an estimation should overestimate the exact behaviour, because the correlations involved

$$\begin{equation} \big[\rho _{\rm a}^{(4)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime }) - \rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime })\rho _{\rm a}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })\big] \end{equation} \tag{ 19 }$$

can reasonably be expected to become rather small compared to ρ⁴, for spatial configurations (r, r', r'', r''') describing a large part of the integration domain $\Lambda _{\rm a}^4$. A precise estimation cannot be performed at this level, so here we only assume that the third term in expression (18) should grow slower than $N_{\rm a}^2$.

The previous simple arguments suggest that the square of fluctuations (18) should become small compared to $N_{\rm a}^2$ in the SL, i.e.

$$\begin{equation} \langle V_{N_{\rm a}}^2 \rangle - [\langle V_{N_{\rm a}} \rangle ]^2 \sim {\rm o}(N_{\rm a}^2) . \end{equation} \tag{ 20 }$$

Note that the corresponding fluctuations for an ordinary system with short-range interactions at thermodynamical equilibrium do satisfy the behaviour (20) since they are proportional to N_a.

In the following, we will have to take averages of quantities involving the potential energy $V_{N_{\rm a}}$. Taking into account the extensivity of its average, as well as the behaviour (20) of its fluctuations, we propose to decompose $V_{N_{\rm a}}$ as

$$\begin{equation} V_{N_{\rm a}} \rightarrow \langle V_{N_{\rm a}} \rangle + W_{N_{\rm a}} , \end{equation} \tag{ 21 }$$

with $\langle V_{N_{\rm a}} \rangle =\mathord {\mathrm{O}}(N_{\rm a})$ and $W_{N_{\rm a}}={\rm o}(N_{\rm a})$ in the SL. We assume that such a decomposition ansatz holds for most probable configurations, which mainly determine the averages of interest, so it should provide the exact behaviours in the cases studied further.

Note that both the extensivity of $\langle V_{N_{\rm a}} \rangle$ and the fluctuation behaviour (20), although quite plausible, are not rigorously proved from a mathematical point of view. In fact, we will show a posteriori that the mass distributions and correlations inferred from the decomposition ansatz are such that the a priori assumptions about $\langle V_{N_{\rm a}} \rangle$ and $(\langle V_{N_{\rm a}}^2 \rangle - [\langle V_{N_{\rm a}} \rangle ]^2)$ are indeed satisfied. This provides some kind of consistency check of the derivations.

As a starting point, we compute the average (13) with the microcanonical distribution for $\mathcal{S}_{\rm a}$. Then, the standard integration over the momenta of the N_a particles leads to

$$\begin{eqnarray} &&\nonumber \rho _{\rm a}(\mathbf {r})=B(E_{\rm a},N_{\rm a},\Lambda _{\rm a}) \int _{\Lambda _{\rm a}^{N_{\rm a}-1},|\mathbf {r}_i - \mathbf {r}_j| > \sigma _{\rm a}} \prod _{i=2}^{N_{\rm a}} \mathrm{d}^3 \mathbf {r}_i [E_{\rm a}-V_{N_{\rm a}}(\mathbf {r},\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})]^{3N_{\rm a}/2-1} \\ &&\times \theta (E_{\rm a}-V_{N_{\rm a}}(\mathbf {r},\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})) . \end{eqnarray} \tag{ 22 }$$

In that expression, θ(ξ) is the usual Heaviside function such that θ(ξ) = 1 for ξ > 0 and θ(ξ) = 0 for ξ < 0, while the normalization constant B(E_a, N_a, Λ_a) ensures that the spatial integral of ρ_a(r) over the box does provide the total mass M_a of the system. The conditions |r_i − r_j| > σ_a, arising from the hard-sphere interaction, apply to any pair of particles, and in particular those including particle 1 fixed at r₁ = r.

For ε > −1/2, it turns out that $(E_{\rm a}-V_{N_{\rm a}}(\mathbf {r},\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}}))$ is positive for any spatial configuration thanks to the upper bound (11) for the potential energy. Then, the Heaviside function can be exactly replaced by 1 in the expression (22). For ε ⩽ −1/2, there are spatial configurations which are forbidden. However, in the SL, the average of the kinetic energy $K_{N_{\rm a}}=E_{\rm a}-V_{N_{\rm a}}(\mathbf {r},\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})$ is necessarily strictly positive and extensive. Thus, for the most probable configurations which provide relatively small fluctuations as a consequence of the decomposition ansatz (21), that quantity is still positive. Thus, it is legitimate to discard the condition $(E_{\rm a}-V_{N_{\rm a}}(\mathbf {r},\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})) > 0$ in the SL for any value of ε, except that corresponding to the collapse into a single ball of course.

The previous simplification is crucial for further transformations. In particular, let us introduce the gravitational potential $\Phi (\mathbf {r}|\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})=\Phi _{N_{\rm a}-1}(\mathbf {r})$ at r created by the (N_a − 1) particles located at $\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}}$. According to the decomposition

$$\begin{equation} V_{N_{\rm a}}(\mathbf {r},\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}}) = V_{N_{\rm a}-1}(\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}}) + m_{\rm a} \Phi (\mathbf {r}|\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}}) , \end{equation} \tag{ 23 }$$

we can rewrite formula (22) in the SL as

$$\begin{eqnarray} &&\nonumber \rho _{\rm a}(\mathbf {r})\sim B(E_{\rm a},N_{\rm a},\Lambda _{\rm a}) \int _{\Lambda ^{N_{\rm a}-1}} \mathrm{d}\mu _{N_{\rm a}-1} \prod _{i=2}^{N_{\rm a}} \theta (|\mathbf {r}_i - \mathbf {r}|/\sigma _{\rm a} -1) \\ &&\times \left[E_{\rm a}-V_{N_{\rm a}-1}\right]^{3/2} \left[1-\frac{m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r}) }{ E_{\rm a}-V_{N_{\rm a}-1} }\right]^{3N_{\rm a}/2-1} , \end{eqnarray} \tag{ 24 }$$

where $\mathrm{d}\mu _{N_{\rm a}-1}$ denotes the unnormalized microcanonical measure for the spatial configurations $(\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})$ of a system made up of (N_a − 1) particles, enclosed in the same spherical box with volume Λ_a and the same energy E_a as the system $\mathcal{S}_{\rm a}$ with N_a particles,

$$\begin{eqnarray} &&\mathrm{d}\mu _{N_{\rm a}-1} = \prod _{i=2}^{N_{\rm a}} \mathrm{d}^3 \mathbf {r}_i \prod _{i<j} \theta (|\mathbf {r}_i - \mathbf {r}_j|/\sigma _{\rm a} -1) [E_{\rm a}-V_{N_{\rm a}-1}(\mathbf {r}_2,\ldots ,\mathbf {r}_{N_{\rm a}})]^{3(N_{\rm a}-1)/2-1} . \end{eqnarray} \tag{ 25 }$$

Let us consider the identity

$$\begin{eqnarray} &&\left[1-\frac{m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r}) }{ E_{\rm a}-V_{N_{\rm a}-1} }\right]^{3N_{\rm a}/2-1} =\exp \left\lbrace (3N_{\rm a}/2-1) \ln \left[1-\frac{m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r}) }{ E_{\rm a}-V_{N_{\rm a}-1} }\right] \right\rbrace . \end{eqnarray} \tag{ 26 }$$

Since $m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r})=\mathord {\mathrm{O}}(1)$ and $E_{\rm a}-V_{N_{\rm a}-1}=\mathord {\mathrm{O}}(N_{\rm a})$, we can expand the logarithm in powers of $m_{\rm a}\Phi /(E_{\rm a}-V_{N_{\rm a}-1})$. After multiplication by the factor (3N_a/2 − 1), we see that the linear term provides a contribution of order $\mathord {\mathrm{O}}(1)$, while higher powers provide vanishing contributions when N_a → ∞. Accordingly, we obtain for any spatial configuration

$$\begin{equation} \left[1-\frac{m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r})}{ E_{\rm a}-V_{N_{\rm a}-1} }\right]^{3N_{\rm a}/2-1} \sim \exp \left\lbrace -\frac{3N_{\rm a}m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r}) }{ 2(E_{\rm a}-V_{N_{\rm a}-1}) } \right\rbrace \end{equation} \tag{ 27 }$$

in the SL. Inserting that asymptotic behaviour inside the rhs of formula (24), we infer the SL behaviour

$$\begin{eqnarray} &&\nonumber \rho _{\rm a}(\mathbf {r}) \sim B(E_{\rm a},N_{\rm a},\Lambda _{\rm a}) \int _{\Lambda ^{N_{\rm a}-1}} \mathrm{d}\mu _{N_{\rm a}-1} \prod _{i=2}^{N_{\rm a}} \theta (|\mathbf {r}_i - \mathbf {r}|/\sigma _{\rm a} -1) \\ &&\times [E_{\rm a}-V_{N_{\rm a}-1}]^{3/2} \exp \left\lbrace - \frac{3N_{\rm a}m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r}) }{ 2(E_{\rm a}-V_{N_{\rm a}-1}) } \right\rbrace . \end{eqnarray} \tag{ 28 }$$

We stress that the extensivity of the potential energy for any spatial configuration, which follows from the particular scaling defined here, plays a crucial role in the derivation of the asymptotic expression (28). Within other scalings, which do not preserve that remarkable extensivity property, the asymptotic behaviour (28) would break down.

If we introduce the kinetic energy $K_{N_{\rm a}-1}$ of the (N_a − 1) particles for a given spatial configuration, the exponential factor on the r.h.s. of the asymptotic expression (28) can be seen as some kind of Boltzmann factor. However, at this level, the corresponding temperature fluctuates. Here, we show that such fluctuations can be neglected on applying the fluctuation ansatz described in section 3.

The integral with the measure $\mathrm{d}\mu _{N_{\rm a}-1}$ in the rhs of formula (28) is proportional to the microcanonical average of the quantity

$$\begin{equation} \prod _{i=2}^{N_{\rm a}} \theta (|\mathbf {r}_i - \mathbf {r}|/\sigma _{\rm a} -1) [E_{\rm a}-V_{N_{\rm a}-1}]^{3/2} \exp \left\lbrace - \frac{3N_{\rm a}m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r}) }{ 2(E_{\rm a}-V_{N_{\rm a}-1}) } \right\rbrace \; \end{equation} \tag{ 29 }$$

for the system with (N_a − 1) particles. For estimating that average in the SL, we can apply to $V_{N_{\rm a}-1}$ the decomposition ansatz (21) with N_a → (N_a − 1). The resulting contribution of fluctuations $W_{N_{\rm a}-1}$ can be omitted at leading order. Moreover, since the corresponding average of $V_{N_{\rm a}-1}$ differs from $\langle V_{N_{\rm a}} \rangle$ by a term of order $\mathord {\mathrm{O}}(1)$ in the SL, all factors $(E_{\rm a}-V_{N_{\rm a}-1})$ can be replaced by $(E_{\rm a}-\langle V_{N_{\rm a}} \rangle )$ at leading order.

Similarly to the potential energy $V_{N_{\rm a}-1}$, the gravitational potential energy $m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r})$ can be reasonably expected to display small relative fluctuations around its average, which is of order $\mathord {\mathrm{O}}(1)$ in the SL. Apart from terms of order o(1), that average reduces to the mean gravitational potential

$$\begin{equation} \phi _{\rm a}(\mathbf {r}) = -\int _{\Lambda _{\rm a}} \mathrm{d}^3 \mathbf {r}^{\prime } \,\rho _{\rm a}(\mathbf {r}^{\prime }) \frac{G}{ |\mathbf {r}^{\prime } - \mathbf {r}|} , \end{equation} \tag{ 30 }$$

created by the whole mass distribution ρ_a(r'). Therefore, in the above average of quantity (29), we can replace $m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r})$ by m_aϕ_a(r) at leading order in the SL.

Since all contributions from fluctuations of both $V_{N_{\rm a}-1}$ and $m_{\rm a}\Phi _{N_{\rm a}-1}(\mathbf {r})$ can be neglected at leading order, we eventually recast the asymptotic behaviour (28) as

$$\begin{eqnarray} &&\rho _{\rm a}(\mathbf {r}) \sim C^{(1)}(E_{\rm a},N_{\rm a},\Lambda _{\rm a}) \left[ \int _{\Lambda _{\rm a}^{N-1}} \mathrm{d}\mu _{N_{\rm a}-1} \prod _{i=2}^{N_{\rm a}} \theta (|\mathbf {r}_i - \mathbf {r}|/\sigma _{\rm a} -1) \right] \exp \left\lbrace -\frac{m_{\rm a}\phi _{\rm a}(\mathbf {r}) }{ T(E_{\rm a},N_{\rm a},\Lambda _{\rm a}) } \right\rbrace ,\nonumber\\ \end{eqnarray} \tag{ 31 }$$

where the microcanonical temperature T(E_a, N_a, Λ_a) is defined through the usual relation

$$\begin{equation} E_{\rm a}-\langle V_{N_{\rm a}} \rangle = \langle K_{N_{\rm a}}\rangle = \frac{3N_{\rm a}T(E_{\rm a},N_{\rm a},\Lambda _{\rm a})}{2} , \end{equation} \tag{ 32 }$$

while C⁽¹⁾(E_a, N_a, Λ_a) is a normalization constant.

According to formula (31), the particle fixed at r interacts with the other particles only via hard-core interactions, while the average gravitational potential plays the role of an external potential V_ext(r) = lim_SLm_aϕ_a(r). On the other hand, the occurrence of the Boltzmann factor exp ( − m_aϕ_a(r)/T(E_a, N_a, Λ_a)) shows that a local thermodynamical equilibrium takes place at temperature T(E_a, N_a, Λ_a). Hence, a local equilibrium in $\mathcal{S}_{\infty }$ emerges, and it is identical to that of thermalized pure hard spheres subjected to an external potential V_ext(r) = lim_SLm_aϕ_a(r). That picture is confirmed by the calculation of the one-body distribution $f_{\rm a} ^{(1)}(\mathbf {r},\mathbf {p})$, which behaves in the SL as

$$\begin{eqnarray} &&f_{\rm a} ^{(1)}(\mathbf {r},\mathbf {p}) \sim \rho _{\rm a}(\mathbf {r}) (m_{\rm a}T(E_{\rm a},N_{\rm a},\Lambda _{\rm a})/\pi )^{3/2} \exp (-\mathbf {p}^2/(2m_{\rm a}T(E_{\rm a},N_{\rm a},\Lambda _{\rm a}))) , \end{eqnarray} \tag{ 33 }$$

while the n-body mass and phase-space distributions are given by simple extensions of formulae (31), (33).

The previous identification of a local thermodynamical equilibrium constitutes the central observation which will allow us to introduce the hydrostatic approach for studying the properties of $\mathcal{S}_{\infty }$, as described in the next section. Notice that the temperature defined through relation (32) remains intensive in the SL, and it reduces to some temperature T_∞ for the infinite system $\mathcal{S}_{\infty }$, where the rescaled dimensionless temperature T_∞/(GM²/NR) is a function T*(ε, η) of the intensive parameters ε and η.

In the SL limit, the local equilibrium in $\mathcal{S}_{\rm a}$ is entirely determined by the local packing fraction η_a(r) = ηρ_a(r)/ρ. In particular, the decay of particle correlations is controlled by the local hard-sphere correlation length [15, 16]

$$\begin{equation} \lambda _{\rm a}(\mathbf {r})=\sigma _{\rm a} \xi _{\rm HS}(\eta _{\rm a}(\mathbf {r})) \; \end{equation} \tag{ 34 }$$

where ξ_HS(η) is some dimensionless function of η. Let us assume a priori that the system remains in a fluid phase at the local level, i.e. the local packing fraction does not exceed some critical value η_LC, for which the system is expected to undergo a phase transition from a liquid to a crystalline phase [15, 16]. Then ξ_HS(η_a(r)) remains bounded everywhere, so λ_a(r) behaves as σ_a and vanishes as $N_{\rm a}^{-2/15}$ in the SL.

Since the mass density ρ_a(r) is expected to vary on a length scale of order R_a, V_ext(r) also varies on distances of order R_a. Thus, the external potential seen by the hard spheres is almost constant over their local correlation length λ_a(r). The corresponding density profile is then easily obtained through a straightforward application of density functional theory [17, 18], where the free energy functional can be merely replaced at leading order by its purely local expression built with the free energy density of a homogeneous system. Corrections due to spatial inhomogeneities of the density profile can be safely neglected here because they are of order $(\lambda _{\rm a}(\mathbf {r})/R_{\rm a})^2 \sim N_{\rm a}^{-2/3}$ in the SL, as estimated from a systematic expansion in powers of the density gradient [19–22]. The familiar hydrostatic equation for the density profile in the external potential V_ext(r) = m_aϕ_a(r) then follows on taking the gradient of the variational equation associated with the local functional. The local pressure comes out when differentiating the local free energy density and it reduces to that of a homogeneous gas of pure hard spheres without gravitational interactions, at temperature T and number density ρ_a(r)/m_a, i.e.

$$\begin{equation} P_{\rm a}(\mathbf {r}) = \frac{T \rho _{\rm a}(\mathbf {r})}{m_{\rm a}} p_{\rm HS}(\eta \rho _{\rm a}(\mathbf {r})/\rho ) , \end{equation} \tag{ 35 }$$

where p_HS is the dimensionless hard-sphere pressure which depends only on the local packing fraction η_a(r) = ηρ_a(r)/ρ.

According to the previous infinite separation of typical lengths in the SL limit, the mass density profile and the corresponding gravitational potential for the infinite system $\mathcal{S}_{\infty }$ are exactly determined within the hydrostatic approach. The corresponding coupled equations can be recast for the dimensionless rescaled quantities g(q; ε, η) = lim_SLρ_a(R_aq)/ρ and ψ(q; ε, η) = lim_SLϕ_a(R_aq)/(GM_a/R_a) of $\mathcal{S}_{\infty }$, as

$$\begin{equation} \nabla _{\mathbf {q}} \left[ g(\mathbf {q};\varepsilon ,\eta ) p_{\rm HS}(\eta g(\mathbf {q};\varepsilon ,\eta ))\right] = -g(\mathbf {q};\varepsilon ,\eta )\nabla _{\mathbf {q}} \psi (\mathbf {q};\varepsilon ,\eta )/T^*(\varepsilon ,\eta ) , \end{equation} \tag{ 36 }$$

and

$$\begin{equation} \psi (\mathbf {q};\varepsilon ,\eta ) = -\frac{3}{4\pi }\int _{q^{\prime } \le 1} \mathrm{d}^3 \mathbf {q}^{\prime } \, \frac{g(\mathbf {q}^{\prime };\varepsilon ,\eta )}{|\mathbf {q}^{\prime } - \mathbf {q}|} . \end{equation} \tag{ 37 }$$

Those equations have to be solved inside the unity sphere together with the rescaled-mass constraint

$$\begin{equation} \frac{3}{4\pi }\int _{q \le 1} \mathrm{d}^3 \mathbf {q}\, g(\mathbf {q};\varepsilon ,\eta ) = 1 , \end{equation} \tag{ 38 }$$

and the rescaled-energy constraint

$$\begin{equation} \frac{3T^*(\varepsilon ,\eta )}{2} + \frac{3}{8\pi }\int _{q \le 1} \mathrm{d}^3 \mathbf {q}\, g(\mathbf {q};\varepsilon ,\eta ) \; \psi (\mathbf {q};\varepsilon ,\eta ) = \varepsilon . \end{equation} \tag{ 39 }$$

We recall that the dimensionless temperature T*(ε, η) has to be determined through the resolution of the whole system of equations, where ε and η are the fixed control parameters.

The properties of the coupled equations (36), (37) have been widely studied in the literature. Once the SL has been taken, the limit of point particles corresponding to η = 0 is now well-behaved, and it is readily obtained by replacing the hard-sphere pressure by its ideal counterpart, namely by setting p_HS(ηg(q; ε, η)) = 1 in the dimensionless hydrostatic equation (36). We then recover Antonov's mean-field theory, for which analytical results are available [2, 24]. For η ≠ 0, various numerical studies have been performed by using approximate simple forms of the hard-sphere equation of state p_HS(η) [25, 26].

Eventually, let us quote that the system does not remain in a fluid phase everywhere if ε becomes too close to the collapse lower bound ε_col(η) related to the minimum (10) of the potential energy. The core region then undergoes crystallization⁴, so the hydrostatic equation should be replaced by some elasticity equation incorporating the stress induced by the gravitational field. Moreover, use of the hydrostatic approach for describing the interface with the fluid region outside the solid core would be rather questionable. Notice that close to the wall of the spherical box, i.e. for r almost equal to R_a, the mass density ρ_a(r) varies on a scale σ_a, so the hydrostatic approach fails in that region. Indeed, use of a purely local functional for the free energy is no longer valid and would lead to a violation of the contact theorem (see e.g. ref. [23]). This does not affect the shape of ρ_a(qR_a) in the SL for fixed q < 1.

The derivations in section 4, which ultimately lead to the hydrostatic equations, involve a priori assumptions about the relative smallness of fluctuations. Here we estimate the contributions of such fluctuations, based on simple modellizations of correlations which follow from the local structure of $\mathcal{S}_{\rm a}$ in the SL limit.

First, we can estimate the contributions of correlations to the average gravitational energy which, according to formula (17), can be recast as

$$\begin{equation} \langle V_{N_{\rm a}} \rangle = V_{\rm self}+ V_{\rm corr} . \end{equation} \tag{ 40 }$$

In that decomposition, V_self is the self-gravitational energy of a sphere with mass density ρ_a(r),

$$\begin{equation} V_{\rm self} = \frac{1}{2}\int _{\Lambda _{\rm a}} \mathrm{d}^3 \mathbf {r}\,\rho _{\rm a}(\mathbf {r})\phi _{\rm a}(\mathbf {r})= -\frac{1}{2} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}\,\mathrm{d}^3 \mathbf {r}^{\prime }\, \rho _{\rm a}(\mathbf {r})\rho _{\rm a}(\mathbf {r}^{\prime }) \frac{G}{ |\mathbf {r}- \mathbf {r}^{\prime }|} , \end{equation} \tag{ 41 }$$

while V_corr is the correlation energy

$$\begin{equation} V_{\rm corr} = -\frac{1}{2} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}\,\mathrm{d}^3 \mathbf {r}^{\prime } \,\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime }) \frac{G}{ |\mathbf {r}- \mathbf {r}^{\prime }|} \end{equation} \tag{ 42 }$$

with the truncated mass distribution or mass correlations,

$$\begin{equation} \rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime }) = \rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime }) - \rho _{\rm a}(\mathbf {r}) \rho _{\rm a}(\mathbf {r}^{\prime }) . \end{equation} \tag{ 43 }$$

In the SL, V_self behaves extensively, and it is determined at leading order by the solutions of the hydrostatic equations for $\mathcal{S}_{\infty }$. Taking into account the local structure of $\mathcal{S}_{\rm a}$, we find that $\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime })$ reduces to −ρ_a(r)ρ_a(r') for |r − r'| < σ_a. Moreover, we expect $\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime })$ to vanish over the hard-sphere local correlation length λ_a(r). As argued above, λ_a(R_aq) becomes proportional to σ_a in the SL. Thus, the contribution of the vicinity of point r = R_aq to

$$\begin{equation} \int _{\Lambda _{\rm a}} \mathrm{d}^3 \mathbf {r}^{\prime } \,\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime }) \frac{G}{ |\mathbf {r}- \mathbf {r}^{\prime }|} \end{equation} \tag{ 44 }$$

is of order $G \rho ^2 \sigma _{\rm a}^2=\mathord {\mathrm{O}}(N_{\rm a}^{-4/15})$ in the SL. The remaining contribution to the integral (44) of points r' such that |r − r'| ≫ σ_a can be roughly estimated by replacing $\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime })$ by a constant times m_aρ_a(r)/Λ_a which does not depend on r'. That spread homogeneous approximation is inspired by the sum rule

$$\begin{equation} \int _{\Lambda _{\rm a}} \mathrm{d}^3 \mathbf {r}^{\prime }\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime }) =-m_{\rm a} \rho _{\rm a}(\mathbf {r}) , \end{equation} \tag{ 45 }$$

which follows from particle conservation. The corresponding contribution to integral (44) is then of order $G m_{\rm a} \rho R_{\rm a}^2/\Lambda _{\rm a}= \mathord {\mathrm{O}}(N_{\rm a}^{-3/5})$ which becomes small compared to that of the region |r − r'| ∼ σ_a. Accordingly, we find that the correlation energy (42) is of order $G \rho ^2 R_{\rm a}^3 \sigma _{\rm a}^2=\mathord {\mathrm{O}}(N_{\rm a}^{1/3})$. Thus, $\langle V_{N_{\rm a}} \rangle$ is indeed extensive, and the rescaled potential energy per particle (12) is entirely given by the self-part V_self, namely

$$\begin{equation} u(\varepsilon ,\eta ) = -\frac{9}{ 32 \pi ^2} \int _{q \le 1,q^{\prime } \le 1} \mathrm{d}^3 \mathbf {q}\,\mathrm{d}^3 \mathbf {q}^{\prime }\, g(\mathbf {q};\varepsilon ,\eta )g(\mathbf {q}^{\prime };\varepsilon ,\eta ) \frac{1}{ |\mathbf {q}- \mathbf {q}^{\prime }|} . \end{equation} \tag{ 46 }$$

Let us consider now the fluctuations $(\langle V_{N_{\rm a}}^2 \rangle - [\langle V_{N_{\rm a}} \rangle ]^2)$ of the potential energy, which are given by formula (18). We can rewrite the first two terms as mean-field-like terms

$$\begin{equation} \frac{1}{2} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}\, \mathrm{d}^3 \mathbf {r}^{\prime } \, \rho _{\rm a}(\mathbf {r})\rho _{\rm a}(\mathbf {r}^{\prime }) \frac{G^2 m_{\rm a}^2 }{ |\mathbf {r}- \mathbf {r}^{\prime }|^2} \end{equation} \tag{ 47 }$$

and

$$\begin{equation} \int _{\Lambda _{\rm a}^3} \mathrm{d}^3 \mathbf {r}\,\mathrm{d}^3 \mathbf {r}^{\prime }\, \mathrm{d}^3 \mathbf {r}^{\prime \prime }\, \rho _{\rm a}(\mathbf {r})\rho _{\rm a}(\mathbf {r}^{\prime })\rho _{\rm a}(\mathbf {r}^{\prime \prime }) \frac{G^2 m }{ |\mathbf {r}- \mathbf {r}^{\prime }||\mathbf {r}- \mathbf {r}^{\prime \prime }|} , \end{equation} \tag{ 48 }$$

plus the corresponding correlation terms associated with $[\rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime })-\rho _{\rm a}(\mathbf {r})\rho _{\rm a}(\mathbf {r}^{\prime })]$ and $[\rho _{\rm a}^{(3)}(\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime })-\rho _{\rm a}(\mathbf {r})\rho _{\rm a}(\mathbf {r}^{\prime })\rho _{\rm a}(\mathbf {r}^{\prime \prime })]$. Since ρ_a(r) does satisfy the scaling property (14) in the SL, the mean-field contributions (47) and (48) are of order $N^0_{\rm a}$ and N_a respectively. The corresponding two- and three-body correlation contributions are readily estimated within a simple modellization of the truncated distributions analogous to the one used above for analysing contribution (42) to the average $\langle V_{N_{\rm a}} \rangle$ itself. They are found to become small compared to their mean-field counterparts in the SL.

It remains to estimate the contribution of the third term in expression (18) of the fluctuations $(\langle V_{N_{\rm a}}^2 \rangle - [\langle V_{N_{\rm a}} \rangle ]^2)$. If we define

$$\begin{equation} n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })=\rho _{\rm a}^{(4)} (\mathbf {r},\mathbf {r}^{\prime },\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })/\rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime }) - \rho _{\rm a}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime }) , \end{equation} \tag{ 49 }$$

we can rewrite that four-body correlation term as

$$\begin{equation} \frac{1}{4} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}\,\mathrm{d}^3 \mathbf {r}^{\prime }\, \rho _{\rm a}^{(2)}(\mathbf {r},\mathbf {r}^{\prime }) \frac{G }{ |\mathbf {r}- \mathbf {r}^{\prime }|}\int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}^{\prime \prime } \,\mathrm{d}^3 \mathbf {r}^{\prime \prime \prime }\, n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime }) \frac{G}{ |\mathbf {r}^{\prime \prime } - \mathbf {r}^{\prime \prime \prime }|} . \end{equation} \tag{ 50 }$$

Like for the case of $\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime })$, we expect, for two given points r and r', $n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })$ to take finite values of order ρ² for spatial configurations such that one or more of the relative distances |r'' − r|, |r'' − r'|, |r''' − r'|, |r''' − r''| is of order σ_a. The corresponding contributions to the integral

$$\begin{equation} \int _{\Lambda _{\rm a}^2} \mathrm{d}^3 \mathbf {r}^{\prime \prime } \,\mathrm{d}^3 \mathbf {r}^{\prime \prime \prime }\, n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime }) \frac{G}{ |\mathbf {r}^{\prime \prime } - \mathbf {r}^{\prime \prime \prime }|} \end{equation} \tag{ 51 }$$

are then readily estimated along lines similar to the above when analysing the two-body correlation term (44). For each of the four regions where one of the points r'' or r''' is close to either r or r', we find a contribution to integral (51) of order $G\rho ^2\sigma _{\rm a}^3 R_{\rm a}^2=\mathord {\mathrm{O}}(1)$. After integration over r and r', the corresponding contribution to the term (50) is $\mathord {\mathrm{O}}(N_{\rm a})$. All the other spatial configurations (r'', r''') for which hard-sphere correlations determine $n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })$, ultimately provide contributions o(N_a) to (50). Eventually, it remains to determine the contributions of regions such that all distances |r'' − r|, |r'' − r'|, |r''' − r'|, |r''' − r''| are large compared to σ_a. As for the case of $\rho _{\rm a}^{(2,{\rm T})}(\mathbf {r},\mathbf {r}^{\prime })$, we again use a spread homogeneous approximation inspired by the sum rule

$$\begin{equation} \int _{\Lambda _{\rm a}} \mathrm{d}^3 \mathbf {r}^{\prime \prime } \,\mathrm{d}^3 \mathbf {r}^{\prime \prime \prime }\, n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime }) =2(3-2N_{\rm a})m_{\rm a}^2 \; ; \end{equation} \tag{ 52 }$$

namely, we replace $n_{\mathbf {r},\mathbf {r}^{\prime }}^{(2)}(\mathbf {r}^{\prime \prime },\mathbf {r}^{\prime \prime \prime })$ by a constant proportional to $N_{\rm a}m_{\rm a}^2/\Lambda _{\rm a}^2$. This provides a contribution of order 1 to integral (51), and ultimately a contribution $\mathord {\mathrm{O}}(N_{\rm a})$ to (50). Accordingly, the four-body correlation term (50) is of order N_a, so fluctuations $(\langle V_{N_{\rm a}}^2 \rangle - [\langle V_{N_{\rm a}} \rangle ]^2)$ do indeed grow more slowly than $N_{\rm a}^2$, in agreement with the assumption (20).

The previous estimations confirm that gravitational interactions can be treated at the mean-field level in the SL limit. In particular, the gravitational energy is indeed dominated by the mean-field self-term (41), as assumed in variational derivations of the hydrostatic equations [2]. Furthermore, fluctuations of that gravitational energy can be discarded. This constitutes a consistency check a posteriori for our derivation of the hydrostatic approach.

The exponentiation formula (27) plays a crucial role in the further emergence of Boltzmann factors. That asymptotic formula works in the SL thanks to extensivity properties of systems $\mathcal{S}_{\rm a}$. For the finite system $\mathcal{S}$, its validity requires

$$\begin{equation} N \left[ \frac{ m \Phi _{N-1}(\mathbf {r}) }{E-V_{N-1}} \right]^2 \ll 1. \end{equation} \tag{ 53 }$$

That condition has to be fulfilled by the most probable configurations which determine the equilibrium state of $\mathcal{S}$. If we assume a priori that the hydrostatic approach is accurate for $\mathcal{S}$, the gravitational potential and gravitational energy involved in condition (53) can be replaced by their average values determined by the rescaled solutions of the hydrostatic equations (36), (37), namely

$$\begin{equation} \Phi _{N-1}(\mathbf {r}) \rightarrow (GM/R) \psi (0;\varepsilon ,\eta ) \end{equation} \tag{ 54 }$$

and

$$\begin{equation} V_{N-1} \rightarrow (GM^2/R) u(\varepsilon ,\eta ) . \end{equation} \tag{ 55 }$$

In the replacement (54), we have set r = 0 because |ϕ(r)| is maximum at r = 0, so the condition (53) will indeed be fulfilled for any r. According to the above replacements, that condition can be recast as

$$\begin{equation} \frac{[\psi (0;\varepsilon ,\eta )]^2}{N [T^*(\varepsilon ,\eta )]^2 } \ll 1 . \end{equation} \tag{ 56 }$$

Not surprisingly, for given values of ε and η, this means that the number of particles N of $\mathcal{S}$ must be large enough. For a given value of N, the temperature T*(ε, η) has to be large enough, or in other words the mean kinetic energy has to be large enough.

Another crucial step in the introduction of Boltzmann factors associated with the average gravitational potential relies on the estimation of the a priori fluctuating exponential

$$\begin{equation} \exp \left[ -\frac{3N m \Phi _{N-1}(\mathbf {r}) }{ 2(E-V_{N-1}) } \right] . \end{equation} \tag{ 57 }$$

Fluctuations can be omitted if the corresponding variations of the exponent

$$\begin{equation} -\frac{3N m \Phi _{N-1}(\mathbf {r}) }{2(E-V_{N-1}) } \end{equation} \tag{ 58 }$$

are small compared to 1. A typical fluctuation of Φ_{N − 1}(r) can be readily estimated by considering the fluctuations of position of one particle j with j > 1 close to the particle fixed at r₁ = r. On average, the corresponding relative distance |r_j − r₁| is of order a = (3/4πn)^1/3, while that distance becomes σ when the two particles are in contact. The resulting variation of mΦ_{N − 1}(r) is of order

$$\begin{equation} \frac{G m^2}{\sigma }. \end{equation} \tag{ 59 }$$

In the corresponding variation of the exponent (58), the kinetic energy (E − V_{N − 1}) can be replaced by its average expressed in terms of temperature T. This leads to the simple condition

$$\begin{equation} \frac{G m^2}{\sigma T} \ll 1 , \end{equation} \tag{ 60 }$$

which can be recast in terms of the dimensionless parameters as

$$\begin{equation} \frac{1}{N^{2/3} \eta ^{1/3} T^*(\varepsilon ,\eta )} \ll 1 . \end{equation} \tag{ 61 }$$

Like the first condition (56), that condition is indeed fulfilled for sufficiently large systems with given values of ε and η. For N given, η cannot be too small and/or temperature T*(ε, η) cannot be too low.

If the two independent conditions (56), (61) are clearly necessary, we expect them to also be sufficient for ensuring the validity of the hydrostatic approach for $\mathcal{S}$. In particular, notice that the condition (61) implies that the system has to be weakly coupled at the local level, namely gravitational interactions between close particles can be neglected. This precisely guarantees that the local equilibrium is entirely governed by hard-core interactions, so the hard-sphere pressure can indeed be used in the hydrostatic balance.

Antonov's mean-field theory [2] can be seen as a particular case of the hydrostatic approach where the local pressure reduces to its ideal gas expression, namely for η = 0. Obviously, the weak-coupling condition (61) cannot be fulfilled if one sets η = 0, in relation to the collapse of the finite system made with point particles. Thus, the corresponding analysis of the validity of the theory relies on the introduction of some realistic finite size σ for the particles, i.e. a finite value of η. Then, if both conditions (56), (61) are satisfied, the validity of Antonov's approach still requires in addition that the corresponding ideal mass distribution has to be close to its hard-sphere counterpart. This implies that ε has to be larger than ε_min ≃ −0.335 [2], since otherwise the ideal hydrostatic equations do not have any solution. Furthermore η has to be sufficiently small. Notice that this necessary requirement might not be sufficient when ε lies in some negative range above ε_min where multiple solutions to the hard-sphere hydrostatic equations exist. Accordingly, if Antonov's theory obviously fails for ε < ε_min, as frequently quoted in the literature, its validity for ε > ε_min is not guaranteed.

The previous validity conditions can be readily tested for astrophysical systems, where the required parameters and/or relevant quantities are roughly determined from observations. We consider two examples, namely globular clusters and dust gases.

Globular clusters

Gas of dust

Here, we consider belts made with dust. Such systems are subjected to a given external gravitational potential created by a planet or a star. Then the conditions (53), (60) for the internal gravitational interactions remain both unchanged and necessary for the validity of the hydrostatic approach. For the extensivity condition (53), we can replace Φ_{N − 1}(r) by the gravitational potential GM/R created by the belt with radius R and mass M. The typical kinetic energy per particle is taken as of order m(δv)² where δv is the velocity dispersion. That condition can be then recast as

$$\begin{equation} R (\delta v)^2 \gg \frac{GM}{\sqrt{N}} . \end{equation} \tag{ 62 }$$

Similarly, the weak-coupling condition (60) can be rewritten as

$$\begin{equation} \delta v \gg \sqrt{\frac{Gm}{\sigma }} . \end{equation} \tag{ 63 }$$

For instance, let us consider Fomalhaut's dust belt with data found by [29]. The measurements give M ≈ 1.1 × 10²⁶ g, while objects have a typical size σ ≈ 10 μm and their inner mass density is of order 2.5 g cm⁻³. This leads to a number of particles N ≈ 1 × 10³⁴. Within those numerical values, the rhs of conditions (62) and (63) are respectively of order 10⁻⁵ m³ s⁻² and 10⁻⁶ m s⁻¹, so both conditions are largely fulfilled. More generally, those conditions would be satisfied for similar astrophysical systems made with dust, in particular those which lead to planet formation.