# Thermodynamics of self-gravitating systems with softened potentials

## I Introduction

The statistical mechanics of self-gravitating systems is amazing. It has been studied since long ago by Antonov , Lynden-Bell and Wood , Thirring , and Kiessling , among others . One reason for the interesting and peculiar behavior of these systems is that they are thermodynamically unstable. The usual thermodynamical limit exists only for those systems which are thermodynamically stable . For a system of $`N_p`$ classical particles interacting via a two body potential $`\varphi (r)`$, a sufficient condition for thermodynamical stability states that there must exist a positive constant $`E_0`$ such that, for each configuration $`\{𝒓_1,\mathrm{},𝒓_{N_p}\}`$, the following inequality is obeyed :

$$\mathrm{\Phi }(𝒓_1,\mathrm{},𝒓_{N_p})=\frac{1}{2}\underset{ij}{}\varphi (|𝒓_i𝒓_j|)N_pE_0.$$
(1)   
In contrast, self-gravitating systems do not possess a proper thermodynamical limit. Moreover, due to the short distance singularity of the gravitational potential, the entropy is not even well defined: it diverges for any value of the energy . To define the thermodynamics of these systems the potential must be regularized at short distances. This can be done in many different ways. Particles endowed with a hard core is one possibility . In this case, the potential is repulsive and singular at short distances. Other popular choices are the so called softened potentials, which are smooth at the origin. As shown by Thirring , the thermodynamical instability is caused neither by the singularity nor by the long range nature of the potential, but is due to the fact that the potential is always attractive<sup>*</sup><sup>*</sup>*In the case of hard core particles the potential is repulsive at short distances and the thermodynamical instability is actually due to the long range forces.. The essential common feature of these purely attractive potentials is the appearance of a phase transition separating a high energy homogeneous phase (HP) from a low energy collapsing phase (CP) . The phase transition takes place in an energy interval with negative microcanonical specific heat. From the dynamical point of view, both phases are also different: the single particle motion is superdiffusive in the CP and ballistic in the HP . The dynamics and statistics of simple low dimensional models with long range attractive forces has been studied in . Their conclusions support the idea of a collapsing phase transition as in the Thirring model. If angular momentum is conserved, the situation could be notably altered .

As mentioned, the usual thermodynamical limit does not exist for unstable systems. To have well defined thermodynamics when the number of particles $`N_p`$ is huge, the following scaling must be considered: the potential energy is rescaled by $`1/N_p`$, and then the energy and entropy scale with $`N_p`$. It has been proved for the canonical ensemble that this scaling reproduces mean field theory exactly in the limit $`N_p\mathrm{}`$ . This means that correlations among two or more particles vanish, and therefore the equilibrium state is characterized by a one particle density only, which minimizes the free energy functional. Although we are not aware of any rigorous proof, we shall assume here that the same holds for the microcanonical ensemble, changing minimization of the free energy by maximization of the entropy functional.

If the troubles caused by the short distance singularity are ignored, it is possible to write down a mean field entropy functional for self-gravitating systems, which depends only on the particle density. This functional is not upper bounded, and, therefore, has no absolute maximum, a reflection of the fact that the entropy is not defined in the finite system. For energies larger than $`E_c0.335GM^2/R`$ there is however a local maximum. Below this energy, no local maximum of the entropy exists . This fact was explained in terms of a transition from the homogeneous isothermal sphere behavior to the CP at $`E_c`$. The transition produces negative specific heat, and was called the gravo-thermal catastrophe . Very recently, it has been pointed out that the low energy phase might be described by a spherically non-symmetric deformation of the singular solution of the isothermal Lane-Emden equation .

The gravitational potential must be modified at short distances to make equilibrium statistical mechanics applicable to self-gravitating systems. As shown by Kiessling for the canonical ensemble , in the limit where the classical gravitational potential is recovered the equilibrium state approaches a particle distribution with all particles collapsing at a single point. The behavior of the system will depend on the scales at which the regularization is effective. There might be regularized potentials which, in the mean field limit, produce the global maximum of their associated entropy close to the solution of the isothermal Lane-Emden equation for energies $`EE_c`$. If this is the case, a collapsing transition should occur at some energy close to $`E_c`$. The CP is expected to be very sensitive to the details of the regularization at short distances and the HP almost insensitive to it. On the other hand, if the regularization is effective only at very short distances, the solutions of the isothermal Lane-Emden equation will be global maxima of the entropy only at very high energies, and therefore the collapsing transition will take place at some energy much larger than $`E_c`$. The critical energy will be higher the smaller the scale at which the regularized potential differs significantly from the unregularized one. A completely similar picture was rigorously established by Kiessling for the canonical ensemble .

In this paper we introduce a convenient new softening procedure for the regularization of the gravitational potential, and we investigate the consequences in the microcanonical thermodynamics of self-gravitating systems. The rest of the article is organized as follows: in Sec. II we introduce the family of potentials to be studied; in Sec. III we derive the mean field equation and its general solution in terms of a set of algebraic equations; Sec. IV is devoted to the discussion of the results and Sec. V to summarize the conclusions.

## II Softened potential

As mentioned in the introduction, the short distance singularity of the gravitational potential causes many troubles. What is more, for a real system such singularity is not physical, since at short distances new physics must be taken into account. Thus, the potential should be modified at short distances to avoid the singularity. In simulations of cosmological problems a widely used potential is the so called Plummer softened :

$$\varphi (r)=\frac{GM^2}{\sqrt{r^2+\sigma ^2}}.$$
(2)   
For $`r\sigma `$, (2) coincides with the gravitational potential. Other softened potentials are those known as spline softened . The equilibrium thermodynamics of systems with these softened potentials has been studied in , and the dynamical effects of softening were considered in . The form of the potential at short distances is arbitrary to a large extent, since we do not know how the new interactions modify it.

The problems with the singularity of the gravitational potential in statistical mechanics disappear if the equilibrium distribution is modified appropriately. For instance, if the approach to equilibrium is collisionless, via violent relaxation, the equilibrium state is described by the Lynden-Bell statistics , whose one-particle distribution is of Fermi-Dirac type, and produces an effective repulsion at short distances .

The regularized potential we propose, which, as we shall discuss, has several remarkable features making it very convenient for thermodynamical purposes, is based on the following identityEquation (3) follows immediately from the sine series expansion of the constant function, $`f(x)=1`$, in the interval $`(0,1)`$.:

$$\frac{1}{x}=\mathrm{\hspace{0.33em}4}\underset{k=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}[(2k1)\pi x]}{(2k1)\pi x},0<x<\mathrm{\hspace{0.17em}1}.$$
(3)   
The singularity at the origin is removed by truncating the series to a given order $`N`$. Let us consider a system of $`N_p`$ particles confined within a sphere of radius $`R`$ in 3 dimensions. The maximum distance between two particles is $`2R`$. Hence, our potential must represent $`1/r`$ for distances $`0r<2R`$. Therefore, we choose the following interaction energy between two particles of mass $`m`$ located at $`𝒓`$ and $`𝒓^{}`$:

$$\varphi (|𝒓𝒓^{}|)=\frac{Gm^2}{R}\mathrm{\hspace{0.17em}2}\underset{k=1}{\overset{N}{}}\varphi _k(|𝒓𝒓^{}|/R),$$
(4)   
where $`\varphi _k(x)=\mathrm{sin}(\omega _kx)/(\omega _kx)`$ are spherical Bessel functions of order zero and $`\omega _k=(2k1)\pi /2`$. Fig. 1 displays the singular potential and the regularized potential with $`N=10`$ and $`N=20`$. A similar expansion has been used to introduce simple models in low dimensions which allow to perform numerical simulations of systems with long range attractive forces with CPU time growing only as the number of particles .

What is remarkable of (4) is that each term obeys the following differential equation:

$$(^2+\frac{\omega _k^2}{R^2})\varphi _k(r/R)=\mathrm{\hspace{0.33em}0},$$
(5)   
so that the potential (4) verifies

$$𝒟_N\varphi (r)=\mathrm{\hspace{0.33em}0},$$
(6)   
where

$$𝒟_N=\underset{k=1}{\overset{N}{}}(^2+\frac{\omega _k^2}{R^2}).$$
(7)   
This relation will prove very useful in the mean field analysis of the next section.

## III Mean Field analysis

It is well known that long range forces suppress fluctuations, and thus in these cases a mean field analysis is accurate or even exact. We expect that, dealing with an unstable system in the scaling regime described in the introduction, the description of the thermodynamical state in terms of a one particle density, neglecting two or more particle correlations, gives the essential physical behavior . We will derive in this section the form of the mean field equation and of the corresponding thermodynamic quantities for a system whose dynamics is governed by a potential of the form (4).

### A Mean Field Equation

Let us consider a system of particles enclosed in a spherical region of radius $`R`$ and volume $`V=4/3\pi R^3`$, with a total mass $`M`$ distributed according to a smooth density $`\rho (𝒓)`$, normalized such that $`d^3r\rho (𝒓)=1`$, and interacting via a two body central potential $`\varphi (|𝒓𝒓^{}|)`$. If the potential is smooth, the entropy per particle in the microcanonical ensemble can be written in terms of the particle density as :

$$𝒮=d^3r\rho (𝒓)\left[\mathrm{ln}V\rho (𝒓)\mathrm{\hspace{0.17em}1}\right]+\frac{3}{2}\mathrm{ln}\left(E\mathrm{\Phi }\right),$$
(8)   
where $`E`$ is the total energy and $`\mathrm{\Phi }`$ is the potential energy:

$$\mathrm{\Phi }[\rho ]=\frac{1}{2}d^3rd^3r^{}\rho (𝒓)\varphi (|𝒓𝒓^{}|)\rho (𝒓^{}).$$
(9)   
The volume $`V`$ entering the first term of the r.h.s in (8) has been included to make the entropy look dimensionally correct, and plays no significant role since it only shifts the entropy by a constant. The physical density is the absolute maximum of (8) under the constraint $`\rho =1`$. Differentiating with respect to $`\rho `$ we arrive at the following integral equation:

$$\mathrm{ln}V\rho (𝒓)=\mu \frac{3}{2}\beta d^3r^{}\varphi (|𝒓𝒓^{}|)\rho (𝒓^{}),$$
(10)   
where $`\beta =1/(E\mathrm{\Phi })`$ and $`\mu `$ is the Lagrange multiplier for the constraint $`\rho =1`$. Defining $`\nu (𝒓)`$ by

$$\rho (𝒓)=\frac{1}{V}\mathrm{exp}[\mu +\nu (𝒓)],$$
(11)   
the constraint is solved by taking

$$e^\mu =\frac{V}{d^3re^{\nu (𝒓)}}.$$
(12)   
Substituting (11) and (12) in (10), we obtain for $`\nu (𝒓)`$:

$$\nu (𝒓)=\frac{3}{2}\beta \frac{d^3r^{}\varphi (|𝒓𝒓^{}|)e^{\nu (𝒓^{})}}{d^3r^{}e^{\nu (𝒓^{})}}.$$
(13)   
If we take for $`\varphi (r)`$ the Newtonian potential, we know that the entropy is not well defined. Nevertheless, it is still possible to start formally with the entropy functional (8), which gives a finite result for any smooth distribution $`\rho (𝒓)`$, but is unbounded (see IV C). There can still exist local maxima, which are then solutions of (13). By expanding the right-hand side of this equation in a series of spherical Bessel functions and truncating after N terms, one would obtain results equivalent to the ones we get using the softened potential.

If we now particularize (13) to the softened potential (4), we see that $`\nu (𝒓)`$ obeys the same differential equation (6) as the potential. Imposing rotational symmetry on $`\nu `$ , we obtain the following ordinary differential equation:

$$\underset{k=1}{\overset{N}{}}\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}+\frac{\omega _k^2}{R^2}\right)\nu (r)=\mathrm{\hspace{0.33em}0}.$$
(14)   
The general solution of this equation is a linear combination of $`\{\mathrm{sin}(\omega _kr)/r\}`$ and $`\{\mathrm{cos}(\omega _kr)/r\}`$. The cosines should be absent from the solution since (13) implies that $`\nu `$ is smooth. (Only at $`T=0`$, i.e., $`\beta =\mathrm{}`$, is $`\nu `$ singular.) Indeed, it is shown explicitly in appendix A that the solution of (13) can be written as

$$\nu (r)=\underset{k=1}{\overset{N}{}}\nu _k\varphi _k(r/R),$$
(15)   
where $`\nu _k`$ are $`N`$ numerical coefficients determined by the following set of equations:

$$\nu _k=\mathrm{\hspace{0.33em}3}\beta \frac{GM^2}{R}\frac{_0^1𝑑xx^2\varphi _k(x)\mathrm{exp}\{_k\nu _k\varphi _k(x)\}}{_0^1𝑑xx^2\mathrm{exp}\{_k\nu _k\varphi _k(x)\}}.$$
(16)   
The integral equation (13) has been reduced to a system of $`N`$ non-linear algebraic equations with $`N`$ unknowns. It can be solved by iteration, for instance with a Newton algorithm (see appendix B for a summary of the method used in this work).

### B Thermodynamical quantities

Using formula (A2) of appendix A, it is straightforward to verify that, for a spherically symmetric mass distribution $`\mathrm{exp}(\mu +\nu (r))`$, the potential energy is given by:

$$\mathrm{\Phi }=\frac{GM^2}{R}\underset{k=1}{\overset{N}{}}\left[\frac{_0^R𝑑rr^2\varphi _k(r/R)e^{\nu (r)}}{_0^R𝑑rr^2e^{\nu (r)}}\right]^2.$$
(17)   
For an equilibrium distribution of the form (15) , Eq. (16) implies

$$\mathrm{\Phi }=\frac{1}{9\beta ^2}\frac{R}{GM^2}\underset{k=1}{\overset{N}{}}\nu _k^2.$$
(18)   
Hence, the total energy is

$$E=\frac{1}{\beta }\frac{1}{9\beta ^2}\frac{R}{GM^2}\underset{k}{}\nu _k^2.$$
(19)   
From (8) we easily obtain the equilibrium entropy:

$$𝒮=\mathrm{ln}\left(_0^1𝑑xx^2\mathrm{exp}(\underset{k}{}\nu _k\varphi _k(x))\right)+\frac{R}{GM^2}\frac{_k\nu _k^2}{3\beta }\frac{3}{2}\mathrm{ln}\beta .$$
(20)   
Since the entropy is stationary under variations of the mass distribution, the inverse temperature is $`1/T𝒮/E=\beta =1/(E\mathrm{\Phi })`$.

## IV Results

In order to present specific numerical results, it is convenient to work with dimensionless quantities. We measure the energy in units of the characteristic energy, $`GM^2/R`$, where $`M`$ is the total mass and $`R`$ the radius of the confining sphere. The dimensionless energy is then $`ϵ=ER/(GM^2)`$. Any other quantity with dimensions of energy (like the temperature $`1/\beta `$ and the potential energy $`\mathrm{\Phi }`$) must be also understood to be expressed in units of $`GM^2/R`$, and, similarly, magnitudes with dimensions of length are given in units of $`R`$. As a matter of terminology, we shall call Newtonian potential (NP) to the unregularized potential, $`GM^2/r`$, Newtonian entropy (NE) to its corresponding entropy, regularized potential (RP) to the potential regularized by Eq. (4), and regularized entropy (RE) to its associated entropy.

### A  N = 10

For values of $`N`$ not too large, computations are easy. Let us describe the case $`N=10`$ in detail. (See appendix B for a summary of the numerical methods used in this work). The solution of Eq. (16) as a function of the energy $`ϵ`$ provides all thermodynamic functions. For each value of $`ϵ`$ we found only one solution, which should then be the absolute maximum of (8). We shall return to this point later on, in Sec. IV C. Fig. 2 displays the inverse temperature $`\beta =1/T`$ versus $`ϵ`$. In the thermodynamics of stable systems, this function must be monotonically decreasing, since the entropy is a convex function of the energy . In the present case, however, $`\beta `$ decreases with $`ϵ`$ in the low and high energy regimes, but it increases for $`ϵ`$ in $`(4.46,0.2)`$ and, consequently, the specific heat is negative in this energy interval. This is a consequence of the instability of the system.

As is usually the case with these systems , the negative specific heat region is associated with a transition to a collapsed phase. To investigate this, let us define an order parameter $`\kappa =R_0/R`$, where $`R_0`$ is the radius of the sphere centered at the origin which contains 95% of the mass (of course the value of 95% is arbitrary). Fig. 3 displays $`\kappa `$ versus $`ϵ`$. At $`ϵ=\mathrm{}`$ the mass is distributed homogeneously, and then $`\kappa =(0.95)^{1/3}0.9830`$. When the energy is reduced, $`\kappa `$ decreases monotonically and slowly. Notice the anomaly at $`ϵ0.335`$; we shall discuss it in Sec. IV B. The collapsing order parameter $`\kappa `$ varies abruptly in the region where the specific heat is negative. It decays from $`\kappa 0.95`$, corresponding to an homogeneous phase, to $`\kappa 0.1`$. In the later case, the mass distribution consists of a small dense core and an homogeneous tenuous halo.

The results of this section are similar to those found by regularizing the potential with hard core repulsions , and to those derived from the Lynden-Bell statistics applied to the unregularized potential . It is remarkable that different regularizations lead to similar results.

### B  Newtonian potential

It is interesting to compare the maximum of the RE with the local maximum of the NE. Substituting $`\varphi (|𝒓𝒓^{}|)=GM^2/|𝒓𝒓^{}|`$ in (13), and using the fact that this potential is a Green function of the Laplacian, we get the following differential equation

$$^2\nu (𝒓)=\frac{3}{2}\beta e^\mu \mathrm{exp}(\nu (𝒓)),$$
(21)   
which, for spherically symmetric $`\nu `$, is equivalent to the isothermal Lane-Emden equation :

$$\frac{d^2\nu (r)}{dr^2}+\frac{2}{r}\frac{d\nu (r)}{dr}+\frac{3}{2}\beta e^\mu e^{\nu (r)}=\mathrm{\hspace{0.33em}0}$$
(22)   
The proper solutions of Eq. (22), with $`\beta `$ and $`\mu `$ such that $`\beta =1/(ϵ\mathrm{\Phi })`$ and $`\mu =\mathrm{ln}𝑑rr^2\mathrm{exp}\nu (r)`$, give local maxima of the entropy if $`ϵ>0.335`$ . The high energy phase should depend only weakly on the form of the potential at short distances. Therefore, the maximum of the RE in the high energy phase might be an approximation to the local maximum of the NE given by Eq. (22). This is indeed the case.

To see how close the maximum of the RE is to the appropriate solution of eq. (22), we define a distance between functions by $`D=\mathrm{max}\{|\nu _{N=10}(r)\nu _{LE}(r)|\}`$, where the subscripts indicate the solutions of Eq. (16) with $`N=10`$, and of Eq. (22) respectively. For $`ϵ>0.335`$, i.e., when the Lane-Emden Eq. determines a local maximum of the NE, $`D<10^4`$. The absolute maximum of the RE is indeed a very good approximation to the local maximum of the NE.

Now, we can understand the anomaly in $`\kappa `$ around $`ϵ0.335`$, which was mentioned in Sec. IV A and which can be appreciated in Fig. 3. At this point, which is close to the energy at which the solutions of the Lane-Emden eq. cease to be local maxima of the NE, the nature of the maximum of the RE also changes, originating anomalies such as the pick in $`1/T`$ (Fig. 2) and the fissure in $`\kappa `$ (Fig. 3).

The effect of the regularization is to deform the entropy functional dramatically for mass distributions $`\rho (r)`$ which are very concentrated at the origin. These distributions get a huge amount of negative entropy after softening the potential, at least for $`N=10`$, in such a way that there are no maxima of the RE close to them. On the other hand, the entropy of smooth distributions which are not concentrated is sensitive to the global form of the potential rather than to the short distance details. Therefore, these distributions have similar NE and RE, and they essentially do not feel the regularization. The solutions of the Lane-Emden Eq. belong to this class, and, consequently, close to them there is a local maximum of the RE which is indeed the global maximum for not too large $`N`$, in particular for $`N=10`$.

### C  $`N`$ dependence

The NP can be arbitrarily well approximated at short distances by a RP with $`N`$ sufficiently large. Consequently, the maximum of the RE close to the local maximum of the NE will attain the later in the $`N\mathrm{}`$ limit, and, obviously, must cease to be the global maximum of the RE and become a local one for some value of $`N`$, which will be denoted by $`N_c`$. Since the maximum of the entropy depends on the energy, $`N_c`$ is a function of $`ϵ`$. In principle, we can compute $`N_c(ϵ)`$ by solving Eq. (16) for large values of $`N`$. In practice, however, this is very difficult and we must content ourselves with an estimate of $`N_c`$.

To get the estimate, let us first analyze how matter distributions with arbitrarily high NE can be built. There is an upper bound for the entropy functional (8) if the potential energy (9) is bounded from below ($`\mathrm{\Phi }\mathrm{\Phi }_{min}`$ for any $`\rho (r)`$):

$$𝒮\mathrm{\hspace{0.33em}1}+\frac{3}{2}\mathrm{ln}(ϵ\mathrm{\Phi }_{min}).$$
(23)   
In the case of the RP, $`\mathrm{\Phi }_{min}=N`$. Since the entropy has an upper bound, it is reasonable to assume that it has a global maximum given by a regular function $`\nu (r)`$ of the form (15), with coefficients $`\nu _k`$ verifying (16). The potential energy associated to the NP has no lower bound and therefore $`\mathrm{\Phi }_{min}=\mathrm{}`$. Hence, (23) does not provide an upper bound for the NE. Indeed, it is straightforward to verify that the distribution

$$\rho (r)=\{\begin{array}{cc}\frac{3\alpha }{4\pi r_0^3}\hfill & 0<r<r_0\hfill \\ \frac{3(1\alpha )}{4\pi (1r_0^3)}\hfill & r_0<r<\mathrm{\hspace{0.17em}1}\hfill \end{array}$$
(24)   
with $`0<\alpha <1`$, has arbitrarily large entropy when we take the limit $`r_00`$, while maintaining $`\alpha \mathrm{ln}r_0`$ constantNotice that in the limit $`r_00`$ with $`\alpha \mathrm{ln}r_0`$ constant only an infinitesimal amount of matter collapses, while the rest is homogeneously distributed. This reflects the fact that it is enough that two particles (hard binaries) become arbitrarily close to make the potential energy arbitrarily negative, and therefore the kinetic energy arbitrarily large. However, as they constitute only two degrees of freedom, their contribution to the purely configurational entropy term $`\rho (\mathrm{ln}V\rho 1)`$ of (8) is negligible.. This is true for any value of $`ϵ`$. We shall call these distributions, for any values of $`\alpha `$ and $`r_0`$, special distributions (SD). If $`N`$ is large enough, there are SD with larger RE than the maximum close to the solution of the Lane-Emden equation..

As already claimed, it is very difficult to get the solutions of eq. (16) for large values of $`N`$. To overcome this problem and obtain an estimate of $`N_c`$, we shall study the restriction of the RE to particle distributions of the form (24) (SD). In such a way, we have a RE which depends only on two parameters, $`\alpha `$ and $`r_0`$. Now, the maximization of this entropy with respect to $`\alpha `$ and $`r_0`$ is an easy task, even for very large values of $`N`$. Obviously, the smaller value of $`N`$ for which the maximum of the restricted RE is larger than the RE of the corresponding solution of the Lane-Emden equation will give an estimate of $`N_c`$. Strictly speaking, this estimate is an upper bound on $`N_c`$.

Before analyzing the behavior with $`N`$, let us look again to the $`N=10`$ case, for which computations are easy. Fig. 4 displays the maximum of the restricted RE and the RE of the corresponding solution of (16), both for $`N=10`$, as a function of $`ϵ`$. The later distribution has always larger RE than any SD. This, besides the fact that we did not find other solutions by varying the initial guess, confirms that for each $`ϵ`$ only a local maximum of the $`N=10`$ RE exists. It is, obviously, the global maximum.

To investigate the behavior with $`N`$, we computed the estimate of the critical $`N_c`$ for several values of $`ϵ`$. As it could have been anticipated, $`N_c`$ grows with $`ϵ`$. Table B shows the results. Column one displays $`ϵ`$, column two the entropy of the solutions of the Lane-Emden eq., column three the maximum of the restriction of the RE to SD for the estimated $`N_c`$, and column four the estimated $`N_c`$. It is apparent that, in the high energy phase ($`ϵ0.335`$), we must go to $`N`$ larger than 30 to see global maxima different from the solutions of the Lane-Emden equation. It is worth noting that a similar scenario was rigorously established by Kiessling for the equilibrium state of self-gravitating systems in contact with a thermal bath .

## V Conclusions

To define the thermodynamics of gravitational systems properly, the Newtonian potential must be regularized at short distances, removing its singularity. Only then is the entropy well defined, or, in a mean field approach, the entropy functional upper bounded. One way to introduce a regularization is by softening, i.e., by making the potential smooth at short distances while keeping it basically unchanged at long distances. There are infinitely many ways to achieve that. One interesting possibility is given by the truncation of the expansion of the gravitational potential in spherical Bessel functions to a given order $`N`$, as in Eq. (4). This regularization has the virtue of reducing the mean field integral equation to a system of $`N`$ algebraic equations with $`N`$ unknowns. This simplifies considerably the solution of the problem.

The result which emerges from this approach is the following: if the regularization is mild enough, $`N<30`$, the system undergoes a phase transition separating a high energy homogeneous phase from a low energy collapsed phase. In the high energy phase, the mass distribution and the thermodynamic quantities are those of an isothermal sphere. Quantitatively, they are very close to the solutions of the corresponding Lane-Emden equation. The low energy phase is characterized by a mass distribution consisting of a dense core surrounded by a tenuous halo. As usual in these cases , the transition from the HP to the CP takes place in an energy interval with negative specific heat, an indication of the thermodynamical instability of the system. These results are remarkably similar to those found with a different regularization (hard core spheres) , and with those derived from the Lynden-Bell statistics applied to the unregularized potential . We can then conclude that the thermodynamics is not very sensitive to the form of the regularization.

The effect of a mild regularization is to deform the entropy functional in such a way that, in the high energy phase, the global maximum of the entropy with the regularized potential is very close to the local maximum with the unregularized potential. The analysis based on the Lane-Emden equation is therefore very accurate and the conclusions extracted from it hold. If the potential is too sharp at short distances, $`N>30`$, we expect also a collapsing transition, which however will take place at a much higher energy than the one predicted by the analysis of the stability of the solutions of the Lane-Emden equation. Below the critical energy the global maximum of the RE will not be in the vicinity of the solution of the Lane-Emden equation, where, nevertheless, there will be a local maximum of the RE. Besides describing such metastable states, the Lane-Emden equation might be physically relevant in diluted self-gravitating systems with an interparticle distance high enough to be insensitive to a truncation of the expansion of the gravitational potential in spherical Bessel functions, Eq. (4), to 30 terms.

Finally, let us comment on the structure of the low energy microcanonical equilibrium state when the regularized potential is very sharp ($`N>30`$) at short distances. In this regime there are high entropy mass distributions consisting of a small amount (infinitesimal when $`N\mathrm{}`$) of matter condensed and the rest homogeneously distributed. This might indicate that at these scales the system is not well described by a smooth density, and granularity is playing a major role.

## A

Let us show that the solutions of Eq. (13) for $`\nu `$ rotationally symmetric are of the form (15) if the potential is $`\varphi (r)=_{k=1}^N\varphi _k(r)`$. Introducing spherical coordinates, eq. (13) can be written

$$\nu (r)=\frac{3\beta GM^2}{R}\underset{k=1}{\overset{N}{}}\frac{_0^1𝑑r^{}r^{}e^{\nu (r^{})}_0^\pi 𝑑\theta \mathrm{sin}\theta \varphi _k(\sqrt{r^2+r^22rr^{}\mathrm{cos}\theta })}{2_0^1𝑑r^{}r^2e^{\nu (r^{})}}.$$
(A1)   
The integral in $`\theta `$ can be readily performed, and gives

$$_0^\pi 𝑑\theta \mathrm{sin}\theta \frac{\mathrm{sin}[\omega _k(r^2+r^22rr^{}\mathrm{cos}\theta )^{1/2}]}{\omega _k(r^2+r^22rr^{}\mathrm{cos}\theta )^{1/2}}=\mathrm{\hspace{0.33em}2}\varphi _k(r)\varphi _k(r^{}).$$
(A2)   
Eqs. (A1) and (A2) imply $`\nu (r)=_{k=1}^N\nu _k\varphi _k(r)`$, with the coefficients $`\nu _k`$ determined by (16), QED.

## B

Since the numerical solution of Eq. (16) is central in this work, we shall outline in this appendix the method used to solve it. The problem is to find the roots of a vector function defined in a multidimensional space. Eq. (16) can be written as

$$F_i(\nu _1,\mathrm{},\nu _N)=\mathrm{\hspace{0.33em}0}$$
(B1)   
with $`i=1,\mathrm{},N`$. We shall use matrix notation and denote by $`𝝂`$ the complete set $`(\nu _1,\mathrm{},\nu _N)`$ and by $`𝑭`$ the vector $`(F_1,\mathrm{},F_N)`$.

To solve systems of equations like (B1) we choose the Newton-Raphson method, which works as follows : provided we have an initial guess $`𝝂`$, which is close to the solution of (B1), we can expand $`F_i`$ in Taylor series in a neighbourhood of $`𝝂`$:

$$F_i(𝝂+\delta 𝝂)=F_i(𝝂)+\underset{j=1}{\overset{N}{}}J_{ij}\delta \nu _j+O(\delta 𝝂^2)$$
(B2)   
where $`J_{ij}=\frac{F_i}{\nu _j}`$ is the Jacobian matrix. In matrix notation we have:

$$𝑭(𝝂+\delta 𝝂)=𝑭(𝝂)+𝑱\delta 𝝂+O(\delta 𝝂^2)$$
(B3)   
By neglecting terms of order higher than linear in $`\delta 𝝂`$ and by setting $`𝑭(𝝂+\delta 𝝂)=0`$, we obtain a set of linear equations for the corrections $`\delta 𝝂`$ that move each function $`F_i`$ closer to zero simultaneously:

$$𝑱\delta 𝝂=𝑭$$
(B4)   
This linear equation is a standar problem in numerical linear algebra and can be solved by LU decomposition. The corrections are then added to the initial guess,

$$𝝂_{new}=𝝂+\delta 𝝂$$
(B5)   
and the process is iterated to convergence. It is possible to show that the method converges always provided the initial guess is close enough to the root. It can also spectacularly fail to converge indicating (though not proving) that the putative root does not exist nearby. To avoid problems with the poor global convergence of the method, we started with many different initial guess. We always found convergence to the same solution, except when $`N`$ was larger than $`30`$, where we found only convergence at high energy. We never got two different solutions (whithin our convergence criterion, see below) by starting at two different points.

Numerically, a convergence criterion is necessary. We stoped computations when one of this two conditions

$$\underset{i=1}{\overset{N}{}}\left|\delta \nu _i\right|<\mathrm{\hspace{0.25em}10}^{10}\mathrm{or}\underset{i=1}{\overset{N}{}}\left|F_i\right|<\mathrm{\hspace{0.25em}10}^{10}$$
(B6)   
was verified.

Each time the function $`F_i`$ was called, several integrals entering eq. (16) were performed numerically, using a Romberg algorithm . The integrands are smooth functions and it was possible to achieve a high precision with a relatively modest numerical effort.