text:
# Thermodynamics and collapse of self-gravitating Brownian particles in 𝐷 dimensions

## 1 Introduction

In an earlier paper , we studied a model of self-gravitating Brownian particles confined within a three-dimensional spherical box. We considered a high friction limit in which the equations of the problem reduce to a Smoluchowski-Poisson system with appropriate constraints ensuring the conservation of energy (in the microcanonical ensemble) or temperature (in the canonical ensemble) . The equilibrium states (maximum entropy states) correspond to isothermal configurations which are known to exist only above a critical energy or a critical temperature (see, e.g., ). When no hydrostatic equilibrium exists, we found that the system generates a finite time singularity (i.e., the central density becomes infinite in a finite time) and we derived self-similar solutions describing the collapse. This study was performed both in the microcanonical and canonical ensembles, with emphasize on the inequivalence of ensembles for such a nonextensive system. In the canonical ensemble, we showed that the scaling exponent for the density is $`\alpha =2`$ and we determined the invariant profile $`f(x)`$, satisfying $`f(x)x^\alpha `$ for $`x+\mathrm{}`$, analytically. In the microcanonical ensemble, the scaling exponent $`\alpha 2.21\mathrm{}`$ and the corresponding invariant profile $`f(x)`$ were determined numerically.

In this paper, we propose to extend our analysis to a space of arbitrary dimension $`D`$. The interest of this extension is twofold. First, we shall consider an infinite dimension limit $`D+\mathrm{}`$ in which the problem can be solved analytically. In particular, it is possible to determine the scaling exponent $`\alpha (D)`$ and the profile $`f(x,D)`$ in the microcanonical ensemble by a systematic expansion procedure in powers of $`D^1`$ (Sec. 3.4), while the canonical value is always $`\alpha =2`$ and the profile can be calculated exactly for any dimension (Sec. 3.3). We show that, already up to order $`O(D^2)`$, the results of the large $`D`$ expansion agree remarkably well with those found numerically for $`D=3`$. Moreover, we show that the nature of the problem changes at two particular dimensions $`D=2`$ and $`D=10`$. In Sec. 2, we compute the equilibrium phase diagram as a function of the dimension. For $`2<D<10`$, the $`TE`$ curve has a spiral shape like in $`3D`$. For $`D>10`$ and $`D<2`$, the $`TE`$ curve is monotonous. The dimension $`D=2`$ is critical and requires a particular attention that is given in Sec. 4. We show that for $`D=2`$ the system generates a Dirac peak (containing a finite fraction of mass) while for $`D>2`$, the central singularity contains no mass (at the collapse time). The case $`D=2`$ has interest in theoretical physics regarding 2D gravity  and string theory  (in connection with the Liouville field theory). It has also applications in the physics of random surfaces  and random potentials , 2D turbulence  and chemotaxis  (for bacterial populations). Finally, the dynamical equations considered in this paper and in Ref.  are receiving a growing interest from mathematicians who established rigorous results concerning the existence and unicity of solutions for arbitrary domain shape without specific symmetry. We refer to the papers of Rosier  and Biler et al. , and references therein, for the connection of our study with mathematical results.

## 2 Equilibrium structure of isothermal spheres in dimension $`D`$

### 2.1 The maximum entropy principle

Consider a system of particles with mass $`m`$ interacting via Newtonian gravity in a space of dimension $`D`$. The particles are enclosed within a box of radius $`R`$ so as to prevent evaporation and make a statistical approach rigorous. Let $`f(𝐫,𝐯,t)`$ denote the distribution function of the system, i.e. $`f(𝐫,𝐯,t)d^D𝐫d^D𝐯`$ gives the mass of particles whose position and velocity are in the cell $`(𝐫,𝐯;𝐫+d^D𝐫,𝐯+d^D𝐯)`$ at time $`t`$. The integral of $`f`$ over the velocity determines the spatial density

(1)
$$\rho =fd^D𝐯.$$

The total mass of the configuration is

(2)
$$M=\rho d^D𝐫.$$

In the meanfield approximation, the total energy of the system can be expressed as

(3)
$$E=\frac{1}{2}fv^2d^D𝐫d^D𝐯+\frac{1}{2}\rho \mathrm{\Phi }d^D𝐫=K+W,$$

where $`K`$ is the kinetic energy and $`W`$ the potential energy. The gravitational potential $`\mathrm{\Phi }`$ is related to the density by the Newton-Poisson equation

(4)
$$\mathrm{\Delta }\mathrm{\Phi }=S_DG\rho ,$$

where $`S_D`$ is the surface of a unit sphere in a $`D`$-dimensional space and $`G`$ is the constant of gravity. Finally, we introduce the Boltzmann entropy

(5)
$$S=f\mathrm{ln}fd^D𝐫d^D𝐯,$$

and the free energy (more precisely the Massieu function)

(6)
$$J=S\beta E,$$

where $`\beta =1/T`$ is the inverse temperature. If the system is isolated, the equilibrium state maximizes the entropy $`S`$ at fixed energy $`E`$ and mass $`M`$ (microcanonical description). Alternatively, if the system is in contact with a heat bath which maintains its temperature fixed, the equilibrium state maximizes the free energy $`J`$ at fixed mass $`M`$ and temperature $`T`$ (canonical description). It can be shown that for systems interacting via a long-range potential like gravity, this meanfield description is exact so that our “thermodynamical approach” is rigorous.

To solve this variational problem, we shall proceed in two steps. We first maximize $`S`$ (resp. $`J`$) at fixed $`M`$, $`E`$ (resp. $`T`$) and $`\rho (𝐫)`$. This yields the Maxwell distribution

(7)
$$f=\frac{1}{(2\pi T)^{D/2}}\rho (𝐫)e^{\frac{v^2}{2T}}.$$

It is now possible to express the energy and the entropy in terms of $`\rho (𝐫)`$ as

(8)
$$E=\frac{D}{2}MT+\frac{1}{2}\rho \mathrm{\Phi }d^D𝐫,$$

(9)
$$S=\frac{D}{2}M\mathrm{ln}T\rho \mathrm{ln}\rho d^D𝐫,$$

where we have omitted unimportant constant terms in the entropy (9). The entropy and the free energy are now functionals of $`\rho (𝐫)`$ and we consider their maximization at fixed energy or temperature. Introducing Lagrange multipliers to satisfy the constraints, the critical points of $`S`$ or $`J`$ are given by the Boltzmann distribution

(10)
$$\rho =Ae^{\beta \mathrm{\Phi }}.$$

Then, the equilibrium state is obtained by solving the Boltzmann-Poisson equation

(11)
$$\mathrm{\Delta }\mathrm{\Phi }=S_DGAe^{\beta \mathrm{\Phi }},$$

and relating the Lagrange multipliers to the appropriate constraints. Note that a similar variational problem occurs in the context of two-dimensional turbulence ($`D=2`$) to characterize large-scale vortices considered as maximum entropy structures .

It is easy to show that there is no global maximum of entropy at fixed mass and energy for $`D>2`$ (see, e.g., Ref.  for $`D=3`$). We can make the entropy diverge to $`+\mathrm{}`$ by approaching an arbitrarily small fraction of particles in the core ($`M_{core}M`$) so that the potential energy goes to $`\mathrm{}`$. Since the total energy is conserved, the temperature must rise to $`+\mathrm{}`$ and this leads to a divergence of the entropy to $`+\mathrm{}`$. Note that if we collapse all particles in the core, the entropy would diverge to $`\mathrm{}`$. Therefore, the formation of a Dirac peak is not thermodynamically favorable in the microcanonical ensemble. For $`D=2`$, there exists a global entropy maximum for all energies. On the other hand, there is no global maximum of free energy at fixed mass and temperature for $`D>2`$ (see, e.g., Ref.  for $`D=3`$) and if $`T<T_c=GM/4`$ for $`D=2`$ (see Appendix A). We can make the free energy $`J`$ diverge to $`+\mathrm{}`$ by collapsing all particles at $`r=0`$. Therefore, a canonical system is expected to form a Dirac peak. For $`D=2`$ and $`T>T_c`$, there exists a global maximum of free energy. For $`D<2`$, there exists a global maximum of entropy and free energy for all accessible values of energy and temperature. We refer to Refs.  for a rigorous proof of these results. When no global maxima of entropy or free energy exist, we can nevertheless look for local maxima since they correspond to metastable states which can be relevant for the considered time scales. Of course, the critical points of entropy at fixed $`E`$ and $`M`$ are the same as the critical point of free energy at fixed $`T`$ and $`M`$. Only the onset of instability (regarding the second order variations of $`S`$ or $`J`$ with appropriate constraints) will differ from one ensemble to the other. For $`D=3`$, this stability problem was considered by Antonov  and Padmanabhan  in the microcanonical ensemble and by Chavanis  in the canonical ensemble, by solving an eigenvalue equation. It was also studied by Lynden-Bell & Wood  and Katz  by using an extension of Poincaré theory of linear series of equilibria. We shall give the generalization of these results in Sec. 2.6 to the case of a system of arbitrary dimension $`D`$.

### 2.2 The $`D`$-dimensional Emden equation

To determine the structure of isothermal spheres, we introduce the function $`\psi =\beta (\mathrm{\Phi }\mathrm{\Phi }_0)`$, where $`\mathrm{\Phi }_0`$ is the gravitational potential at $`r=0`$. Then, the density field can be written

(12)
$$\rho =\rho _0e^\psi ,$$

where $`\rho _0`$ is the central density. Introducing the notation $`\xi =(S_DG\beta \rho _0)^{1/2}r`$ and restricting ourselves to spherically symmetrical configurations (which maximize the entropy for a non-rotating system), the Boltzmann-Poisson Eq. (11) reduces to the form

(13)
$$\frac{1}{\xi ^{D1}}\frac{d}{d\xi }\left(\xi ^{D1}\frac{d\psi }{d\xi }\right)=e^\psi ,$$

which is the $`D`$-dimensional generalization of the Emden equation . For $`D>2`$, Eq. (13) has a simple explicit solution, the singular sphere

(14)
$$e^{\psi _s}=\frac{2(D2)}{\xi ^2}.$$

The regular solution of Eq. (13) satisfying the boundary conditions

(15)
$$\psi =\psi ^{}=0\mathrm{at}\xi =0,$$

must be computed numerically. For $`\xi 0`$, we can expand the solution in Taylor series and we find that

(16)
$$\psi =\frac{1}{2D}\xi ^2\frac{1}{8D(D+2)}\xi ^4+\frac{1}{24}\frac{D+1}{D^2(D+2)(D+4)}\xi ^6+\mathrm{}$$

To obtain the asymptotic behavior of the solutions for $`\xi +\mathrm{}`$, we note that the transformation $`t=\mathrm{ln}\xi `$, $`\psi =2\mathrm{ln}\xi z`$ changes Eq. (13) in

(17)
$$\frac{d^2z}{dt^2}+(D2)\frac{dz}{dt}=e^z+2(D2).$$

For $`D>2`$, this corresponds to the damped oscillations of a fictitious particle in a potential $`V(z)=e^z2(D2)z`$, where $`z`$ plays the role of the position and $`t`$ the role of time. For $`t+\mathrm{}`$, the particle will come at rest at the bottom of the well at position $`z_0=\mathrm{ln}(2(D2))`$. Returning to original variables, we find that

(18)
$$e^\psi \frac{2(D2)}{\xi ^2}=e^{\psi _s},\mathrm{for}\xi +\mathrm{}.$$

Therefore, the regular solution of the Emden equation (13) behaves like the singular solution for $`\xi +\mathrm{}`$. To determine the next order correction, we set $`z=z_0+z^{}`$ with $`z^{}1`$. Keeping only terms that are linear in $`z^{}`$, Eq. (17) becomes

(19)
$$\frac{d^2z^{}}{dt^2}+(D2)\frac{dz^{}}{dt}+2(D2)z^{}=0.$$

The discriminant associated with this equation is $`\mathrm{\Delta }=(D2)(D10)`$. It exhibits two critical dimensions $`D=2`$ and $`D=10`$. For $`2<D<10`$, we have

(20)
$$e^\psi =\frac{2(D2)}{\xi ^2}\left\{1+\frac{A}{\xi ^{\frac{D2}{2}}}\mathrm{cos}\left(\frac{\sqrt{(D2)(10D)}}{2}\mathrm{ln}\xi +\delta \right)\right\},(\xi +\mathrm{}).$$

The density profile (20) intersects the singular solution (14) infinitely often at points that asymptotically increase geometrically in the ratio $`1:e^{2\pi /\sqrt{(D2)(10D)}}`$ (see, e.g., Fig. 1 of Ref.  for $`D=3`$). For $`D10`$, we have

(21)
$$e^\psi =\frac{2(D2)}{\xi ^2}\left\{1+\frac{1}{\xi ^{\frac{D2}{2}}}\left(A\xi ^{\frac{\sqrt{(D2)(10D)}}{2}}+B\xi ^{\frac{\sqrt{(D2)(10D)}}{2}}\right)\right\},(\xi +\mathrm{}).$$

For $`D=2`$, Eq. (17) can be solved explicitly and we get

(22)
$$e^\psi =\frac{1}{(1+\frac{1}{8}\xi ^2)^2}.$$

This result has been found by various authors in different fields (see, e.g., ). Note that $`e^\psi \xi ^4`$ at large distances instead of the usual $`\xi ^2`$ behavior obtained for $`D>2`$. This implies that the mass of an unbounded isothermal sphere is finite in $`D=2`$, although it is infinite for $`D>2`$.

For $`D<2`$, we can neglect the r.h.s. in Eq. (18) and we get

(23)
$$e^\psi e^{A_D\xi ^{2D}},(\xi +\mathrm{}),$$

where $`A_D`$ is a constant depending on the dimension $`D`$. For $`D=1`$, Eq. (13) can be solved exactly, yielding the result

(24)
$$e^\psi =\frac{1}{\mathrm{cosh}^2(\xi /\sqrt{2})},$$

establishing $`A_1=\sqrt{2}`$.

### 2.3 The Milne variables

As is well-known , isothermal spheres satisfy a homology theorem: if $`\psi (\xi )`$ is a solution of the Emden equation, then $`\psi (A\xi )2\mathrm{ln}A`$ is also a solution, with $`A`$ an arbitrary constant. This means that the profile of isothermal configurations is always the same (characterized intrinsically by the function $`\psi `$), provided that the central density and the typical radius are rescaled appropriately. Because of this homology theorem, the second order differential equation (13) can be reduced to a first order differential equation for the Milne variables

(25)
$$u=\frac{\xi e^\psi }{\psi ^{}}\mathrm{and}v=\xi \psi ^{}.$$

Taking the logarithmic derivative of $`u`$ and $`v`$ with respect to $`\xi `$ and using Eq. (13), we get

(26)
$$\frac{1}{u}\frac{du}{d\xi }=\frac{1}{\xi }(Dvu),$$

(27)
$$\frac{1}{v}\frac{dv}{d\xi }=\frac{1}{\xi }(2D+u).$$

Taking the ratio of the foregoing equations, we obtain

(28)
$$\frac{u}{v}\frac{dv}{du}=\frac{2D+u}{Duv}.$$

The solution curve in the $`(u,v)`$ plane is plotted in Fig. 1 for different values of $`D`$. The curve is parameterized by $`\xi `$. It starts from the point $`(u,v)=(D,0)`$ with a slope $`(dv/du)_0=(D+2)/D`$ corresponding to $`\xi =0`$. The points of horizontal tangent are determined by $`u=D2`$ and the points of vertical tangent by $`u+v=D`$. These two lines intersect at $`(u_s,v_s)=(D2,2)`$, which corresponds to the singular solution (14). For $`2<D<10`$, the solution curve spirals indefinitely around the point $`(u_s,v_s)`$. For $`D10`$, the curve reaches the point $`(u_s,v_s)`$ without spiraling. For $`D=2`$, we have the explicit solution $`v=2(2u)`$ so that $`(u,v)(0,4)`$ for $`\xi +\mathrm{}`$. For $`D<2`$, $`(u,v)(0,+\mathrm{})`$ for $`\xi +\mathrm{}`$ (see Fig. 2). More precisely,

(29)
$$\frac{ue^{\frac{v}{2D}}}{v^{\frac{D}{2D}}}\omega _D,(\xi +\mathrm{}),$$

where we have defined $`\omega _D=1/(A_D(2D)^{\frac{2}{2D}})`$. For $`D=1`$, $`\omega _1=1/\sqrt{2}`$.

### 2.4 The thermodynamical parameters

For bounded isothermal systems, the solution of Eq. (13) is terminated by the box at a normalized radius given by $`\alpha =(S_D\beta \rho _0)^{1/2}R`$. We shall now relate the parameter $`\alpha `$ to the temperature and energy. According to the Poisson equation (4), we have

(30)
$$GM=G\rho d^D𝐫=_0^RS_DG\rho r^{D1}𝑑r=_0^R\frac{d}{dr}\left(r^{D1}\frac{d\mathrm{\Phi }}{dr}\right)𝑑r=\left(r^{D1}\frac{d\mathrm{\Phi }}{dr}\right)_{r=R}.$$

Introducing the dimensionless variables defined previously (using $`r/R=\xi /\alpha `$), we get

(31)
$$\eta \frac{\beta GM}{R^{D2}}=\alpha \psi ^{}(\alpha ).$$

We note that, for $`D=2`$, the parameter $`\eta `$ is independent on $`R`$. This is a consequence of the logarithmic form of the Newtonian potential in two dimensions.

The computation of the energy is a little more intricate. First, using an extension of the potential tensor theory developed by Chandrasekhar (see, e.g., Ref. ), we can show that the potential energy in $`D`$-dimensions can be written

(32)
$$W=\frac{1}{D2}\rho 𝐫\mathrm{\Phi }d^D𝐫,$$

for $`D2`$. Now, the Boltzmann-Poisson equation (11) is equivalent to the condition of hydrostatic equilibrium

(33)
$$p=\rho \mathrm{\Phi },$$

with an equation of state $`p=\rho T`$. Substituting this relation in Eq. (32) and integrating by parts, we obtain

(34)
$$2K+(D2)W=DV_DR^Dp(R),$$

where $`V_D=S_D/D`$ is the volume of a hypersphere with unit radius. Eq. (34) is the form of the Virial theorem in $`D`$-dimension. The total energy $`E=K+W`$ can be written

(35)
$$E=\frac{D4}{D2}K+\frac{D}{D2}V_DR^Dp(R).$$

Expressing the pressure in terms of the Emden function, using $`p=\rho T`$ and Eq. (12), and using Eq. (31) to eliminate the temperature, we finally obtain

(36)
$$\mathrm{\Lambda }\frac{ER^{D2}}{GM^2}=\frac{D(4D)}{2(D2)}\frac{1}{\alpha \psi ^{}(\alpha )}\frac{1}{D2}\frac{e^{\psi (\alpha )}}{\psi ^{}(\alpha )^2}.$$

It turns out that the normalized temperature and the normalized energy can be expressed very simply in terms of the values of the Milne variables at the normalized box radius. Indeed, writing $`u_0=u(\alpha )`$ and $`v_0=v(\alpha )`$ and using Eqs. (31)(36), we get

(37)
$$\eta =v_0,$$

(38)
$$\mathrm{\Lambda }=\frac{1}{v_0}\left[\frac{D(4D)}{2(D2)}\frac{u_0}{D2}\right].$$

The curves $`\eta (\alpha )`$, $`\mathrm{\Lambda }(\alpha )`$ are plotted in Figs. 3-4. For $`2<D<10`$, they exhibit damped oscillations toward the values $`\eta _s=2`$ and $`\mathrm{\Lambda }_s=1/(D2)D/4`$, corresponding to the singular solution (14). For $`D10`$ the curves are monotonous. For $`D=2`$, we have explicitly

(39)
$$\eta =\frac{\alpha ^2}{2(1+\frac{1}{8}\alpha ^2)},\mathrm{\Lambda }=\frac{2}{\alpha ^2}\left(1+\frac{\alpha ^2}{8}\right)\left\{\frac{8}{\alpha ^2}\left(1+\frac{\alpha ^2}{8}\right)\mathrm{ln}\left(1+\frac{\alpha ^2}{8}\right)2\right\}.$$

The expression of the energy has been obtained directly from Eq. (8) with the boundary condition $`\mathrm{\Phi }(R)=0`$. The inverse temperature increases monotonically with $`\alpha `$ up to the value $`\eta _c=4`$. Using Eq. (22) and returning to the original variables, we can write the density profile in the form

(40)
$$\rho =\frac{4M}{\pi R^2(4\eta )(1+\frac{\eta }{4\eta }\frac{r^2}{R^2})^2}.$$

This density profile is represented in Fig. 5 for different temperatures. At the critical inverse temperature $`\eta _c=4`$, all the particles are concentrated at the center of the domain. The density profile approaches the Dirac distribution

(41)
$$\rho (𝐫)M\delta (𝐫)\mathrm{for}\eta \eta _c=4,$$

which has an infinite (negative) energy.

For $`D<2`$, the curves $`\eta (\alpha )`$ and $`\mathrm{\Lambda }(\alpha )`$ tend to $`+\mathrm{}`$ and $`0`$ as $`\alpha +\mathrm{}`$. For $`D=1`$, we have explicitly

(42)
$$\eta =2\sqrt{2}\alpha \mathrm{tanh}(\alpha /\sqrt{2}),\mathrm{\Lambda }=\frac{3}{4\sqrt{2}}\frac{1}{\alpha \mathrm{tanh}(\alpha /\sqrt{2})}+\frac{1}{8}\frac{1}{\mathrm{sinh}^2(\alpha /\sqrt{2})}.$$

Note that according to Eq. (32), the energy is necessarily positive for $`D<2`$, so the region $`\mathrm{\Lambda }0`$ is forbidden. Returning to original variables, the density profile is given by

(43)
$$\rho =\frac{M}{4\sqrt{2}R}\frac{\alpha }{\mathrm{tanh}(\alpha /\sqrt{2})}\frac{1}{\mathrm{cosh}^2\left(\frac{\alpha r}{\sqrt{2}R}\right)},$$

where we recall that $`S_1=2`$. For $`\alpha +\mathrm{}`$, the profile tends to a Dirac peak $`M\delta (r)`$.

In Figs. 6-7, we have plotted the equilibrium phase diagram $`\mathrm{\Lambda }\eta `$, giving the temperature as a function of the energy, for different dimensions $`D`$. For $`2<D<10`$, the curve spirals around the limit point $`(\mathrm{\Lambda }_s,\eta _s)`$ corresponding to the singular solution. For $`D10`$, the curve is monotonous until the limit point. For $`D=2`$, the curve is explicitly given by

(44)
$$\mathrm{\Lambda }=\frac{1}{\eta }\left[\frac{4}{\eta }\mathrm{ln}\left(\frac{4}{4\eta }\right)2\right],$$

and is represented in Fig. 7, together with the case $`D=1`$.

We stress that the preceding results, obtained in the meanfield approximation, are exact in the thermodynamical limit $`N+\mathrm{}`$ such that $`\eta `$ and $`\mathrm{\Lambda }`$ are kept fixed. If the box radius is given, this implies that $`TN`$ and $`EN^2`$. Alternatively, if the temperature and the energy per particle are given, the thermodynamical limit is such that $`N+\mathrm{}`$ with $`N/R^{D2}`$ constant (for $`D>2`$).

### 2.5 The minimum temperature and minimum energy

For $`2<D<10`$, the curve $`\eta (\alpha )`$ presents an extremum at points $`\alpha _n`$ such that $`d\eta /d\alpha (\alpha _n)=0`$. Using Eqs. (37) and (27), we find that this condition is equivalent to

(45)
$$u_0=D2=u_s.$$

Since the curve $`u=u_s`$ passes through the center of the spiral in the $`(u,v)`$ plane, this determines an infinity of solutions (see Fig. 8), one at each extremum of $`v`$ (since $`\eta =v_0`$). Asymptotically, the $`\alpha _n`$ follow a geometric progression (see Ref.  for more details):

(46)
$$\alpha _ne^{2\pi n/\sqrt{(D2)(10D)}},(n+\mathrm{},\mathrm{integer}).$$

In Fig. 3, we see that an equilibrium state exists only for

(47)
$$\eta =\frac{\beta GM}{R^{D2}}\eta (\alpha _1),(2<D<10).$$

This determines a maximum mass (for given $`T`$ and $`R`$) or a minimum temperature (for given $`M`$ and $`R`$) beyond which no equilibrium state is possible. In that case, the system is expected to undergo an isothermal collapse (see Sect. 3.3). For $`D=2`$ and for $`D10`$, the $`\eta (\alpha )`$ curve is monotonous. An equilibrium state exists only provided that

(48)
$$\eta =\beta GM\eta _c=4,(D=2),$$

(49)
$$\eta =\frac{\beta GM}{R^{D2}}\eta _s=2,(D10).$$

We get comparable results for the energy. For $`2<D<10`$, the curve $`\mathrm{\Lambda }(\alpha )`$ presents an extremum at points $`\alpha _n^{}`$ such that $`d\mathrm{\Lambda }/d\alpha (\alpha _n^{})=0`$. Using Eqs. (38) and (26)(27), we find that this condition is equivalent to

(50)
$$4u_0^2+2u_0v_0+(D^28D+4)u_0+D(D2)(4D)=0.$$

We can check that the limit point $`(u_s,v_s)`$ is solution of this equation. Therefore, the intersection of the parabola $`(P)`$ defined by Eq. (50) with the spiral in the $`(u,v)`$ plane determines an infinity of points $`\alpha _n^{}`$ at which the energy is extremum (see Fig. 8). On Fig. 4, we see that an equilibrium state exists only for

(51)
$$\mathrm{\Lambda }=\frac{ER^{D2}}{GM^2}\mathrm{\Lambda }(\alpha _1^{}),(2<D<10).$$

This determines a minimum energy (for given $`M`$ and $`R`$) or a maximum radius (for given $`M`$ and $`E`$) beyond which no equilibrium state exists. In that case, the system is expected to collapse and overheat; this is called gravothermal catastrophe (see Sect. 3.4). For $`D10`$, the curve $`\mathrm{\Lambda }(\alpha )`$ is monotonous. An equilibrium state exist only for

(52)
$$\mathrm{\Lambda }=\frac{ER^{D2}}{GM^2}\mathrm{\Lambda }_s=\frac{1}{D2}\frac{D}{4},(D10).$$

For $`D=2`$, there exists an equilibrium state for each value of energy (see Fig. 7): there is no gravothermal catastrophe in the microcanonical ensemble in two-dimensions . For $`D<2`$, there exists an equilibrium state for all accessible values of energy and temperature.

### 2.6 The thermodynamical stability

We now study the thermodynamical stability of self-gravitating systems in various dimensions. We start by the canonical ensemble which is simpler in a first approach. A critical point of free energy at fixed mass and temperature is a local maximum if, and only if, the second order variations

(53)
$$\delta ^2J=\frac{(\delta \rho )^2}{2\rho }d^D𝐫\frac{1}{2T}\delta \rho \delta \mathrm{\Phi }d^D𝐫,$$

are negative for any variation $`\delta \rho `$ that conserves mass to first order, i.e.

(54)
$$\delta \rho d^D𝐫=0.$$

This is the condition of thermodynamical stability in the canonical ensemble. Introducing the function $`q(r)`$ by the relation

(55)
$$\delta \rho =\frac{1}{S_Dr^{D1}}\frac{dq}{dr},$$

and following a procedure similar to the one adopted in Ref. , we can put the second order variations of free energy in the quadratic form

(56)
$$\delta ^2J=\frac{1}{2}_0^R𝑑rq\left[\frac{G}{Tr^{D1}}+\frac{d}{dr}\left(\frac{1}{S_D\rho r^{D1}}\frac{d}{dr}\right)\right]q.$$

The second order variations of free energy can be positive (implying instability) only if the differential operator which occurs in the integral has positive eigenvalues. We need therefore to consider the eigenvalue problem

(57)
$$\left[\frac{d}{dr}\left(\frac{1}{S_D\rho r^{D1}}\frac{d}{dr}\right)+\frac{G}{Tr^{D1}}\right]q_\lambda (r)=\lambda q_\lambda (r).$$

with $`q_\lambda (0)=q_\lambda (R)=0`$. If all the eigenvalues $`\lambda `$ are negative, then the critical point is a maximum of free energy. If at least one eigenvalue is positive, the critical point is an unstable saddle point. The point of marginal stability, i.e. the value of $`\alpha `$ in the series of equilibria $`\eta (\alpha )`$ at which the solutions pass from local maxima of free energy to unstable saddle points, is determined by the condition that the largest eigenvalue is equal to zero ($`\lambda =0`$). We thus have to solve the differential equation

(58)
$$\frac{d}{dr}\left(\frac{1}{S_D\rho r^{D1}}\frac{dF}{dr}\right)+\frac{GF}{Tr^{D1}}=0,$$

with $`F(0)=F(R)=0`$. Introducing the dimensionless variables defined previously, we can rewrite this equation in the form

(59)
$$\frac{d}{d\xi }\left(\frac{e^\psi }{\xi ^{D1}}\frac{dF}{d\xi }\right)+\frac{F(\xi )}{\xi ^{D1}}=0,$$

with $`F(0)=F(\alpha )=0`$. If

(60)
$$\frac{d}{d\xi }\left(\frac{e^\psi }{\xi ^{D1}}\frac{d}{d\xi }\right)+\frac{1}{\xi ^{D1}}$$

denotes the differential operator which occurs in Eq. (59), we can check by using the Emden Eq. (13) that

(61)
$$(\xi ^{D1}\psi ^{})=\psi ^{},(\xi ^De^\psi )=(D2)\psi ^{}.$$

Therefore, the general solution of Eq. (59) satisfying the boundary conditions at $`\xi =0`$ is

(62)
$$F(\xi )=c_1(\xi ^De^\psi (D2)\xi ^{D1}\psi ^{}).$$

Using Eq. (62) and introducing the Milne variables (25), the condition $`F(\alpha )=0`$ can be written

(63)
$$u_0=D2.$$

This relation determines the points at which a new eigenvalue becomes positive (crossing the line $`\lambda =0`$). Comparing with Eq. (45), we see that a mode of stability is lost each time that $`\eta `$ is extremum in the series of equilibria, in agreement with the turning point criterion of Katz  in the canonical ensemble. In particular, the series becomes unstable at the point of minimum temperature (or maximum mass) $`\alpha _1`$. Secondary modes of instability appear at values $`\alpha _2`$, $`\alpha _3`$,… We obtain the same results by considering the dynamical stability of isothermal gaseous spheres with respect to the Navier-Stokes equations (see Ref.  for $`D=3`$). Therefore, dynamical and thermodynamical stability criteria coincide for isothermal gaseous spheres.

According to Eq. (55), the perturbation profile that triggers a mode of instability at the critical point $`\lambda =0`$ is given by

(64)
$$\frac{\delta \rho }{\rho _0}=\frac{1}{S_D\xi ^{D1}}\frac{dF}{d\xi },$$

where $`F(\xi )`$ is given by Eq. (62). Introducing the Milne variables (25), we get

(65)
$$\frac{\delta \rho }{\rho }=\frac{c_1}{S_D}(2v).$$

The density perturbation $`\delta \rho `$ becomes zero at point(s) $`\xi _i`$ such that $`v(\xi _i)=2`$. The number of zeros is therefore given by the number of intersections between the spiral in the $`(u,v)`$ plane and the line $`v=2`$ (see Fig. 8). For the $`n`$-th mode of instability we need to follow the spiral to the $`n`$-th extremum of $`v`$ (since $`\alpha _n`$ corresponds to an extremum of $`\eta `$, hence $`v`$). Therefore, the density perturbation $`\delta \rho `$ corresponding to the $`n`$-th mode of instability has $`n`$ zeros $`\xi _1,\xi _2,\mathrm{},\xi _n<\alpha _n`$. Asymptotically, the zeros follow a geometric progression with ratio $`e^{2\pi /\sqrt{(D2)(10D)}}`$ .

In the microcanonical ensemble, the condition of thermodynamical stability requires that the equilibrium state is an entropy maximum at fixed mass and energy. This condition can be written

(66)
$$\delta ^2S=\frac{(\delta \rho )^2}{2\rho }d^D𝐫\frac{1}{2T}\delta \rho \delta \mathrm{\Phi }d^D𝐫\frac{1}{DMT^2}\left(\mathrm{\Phi }\delta \rho d^D𝐫\right)^2<0,$$

for any variation $`\delta \rho `$ that conserves mass to first order (the conservation of energy has already been taken into account in obtaining Eq. (66)). Now, following a procedure similar to that of Ref.  in $`D=3`$, the second variations of entropy can be put in a quadratic form

(67)
$$\delta ^2S=_0^R_0^R𝑑r𝑑r^{}q(r)K(r,r^{})q(r^{}),$$

with

(68)
$$K(r,r^{})=\frac{1}{DMT^2}\frac{d\mathrm{\Phi }}{dr}(r)\frac{d\mathrm{\Phi }}{dr}(r^{})+\frac{1}{2}\delta (rr^{})\left[\frac{G}{Tr^{D1}}+\frac{d}{dr}\left(\frac{1}{S_D\rho r^{D1}}\frac{d}{dr}\right)\right].$$

The problem of stability can therefore be reduced to the study of the eigenvalue equation

(69)
$$_0^R𝑑r^{}K(r,r^{})F_\lambda (r^{})=\lambda F_\lambda (r),$$

with $`F_\lambda (0)=F_\lambda (R)=0`$. The point of marginal stability ($`\lambda =0`$) will be determined by solving the differential equation

(70)
$$\frac{d}{dr}\left(\frac{1}{S_D\rho r^{D1}}\frac{dF}{dr}\right)+\frac{GF}{Tr^{D1}}=\frac{2V}{DMT^2}\frac{d\mathrm{\Phi }}{dr}(r),$$

with

(71)
$$V=_0^R\frac{d\mathrm{\Phi }}{dr}(r^{})F(r^{})𝑑r^{}.$$

Introducing the dimensionless variables defined previously, this is equivalent to

(72)
$$\frac{d}{d\xi }\left(\frac{e^\psi }{\xi ^{D1}}\frac{dF}{d\xi }\right)+\frac{F}{\xi ^{D1}}=\chi \frac{d\psi }{d\xi },$$

with

(73)
$$\chi =\frac{2}{D\alpha ^{D1}\psi ^{}(\alpha )}_0^\alpha \frac{d\psi }{d\xi }(\xi ^{})F(\xi ^{})𝑑\xi ^{},$$

and $`F(0)=F(\alpha )=0`$. Using the identities (61), we can check that the general solution of Eq. (72) satisfying the boundary conditions for $`\xi =0`$ and $`\xi =\alpha `$ is

(74)
$$F(\xi )=\frac{\chi }{D2u_0}\left(\xi ^De^\psi (D2)\xi ^{D1}\psi ^{}\right)+\chi \xi ^{D1}\psi ^{},$$

The point of marginal stability is then obtained by substituting the solution (74) in Eq. (73). Using the identities

(75)
$$_0^\alpha \psi ^{}\xi ^De^\psi 𝑑\xi =\alpha ^{D1}\psi ^{}(\alpha )(Du_0),$$

(76)
$$(D2)_0^\alpha \xi ^{D1}(\psi ^{})^2𝑑\xi =\alpha ^{D1}\psi ^{}(\alpha )(2D2u_0v_0),$$

which result from simple integrations by parts and from the properties of the Emden equation (13), it is found that the point of marginal stability is determined by the condition (50). Therefore, the series of equilibria becomes unstable at the point of minimum energy in agreement with the turning point criterion of Katz  in the microcanonical ensemble. The structure of the perturbation that triggers the instability can be determined with the graphical construction described in Ref. . It is found that the first mode of instability has a core-halo structure (i.e., two nodes) in the microcanonical ensemble, unlike the first mode of instability in the canonical ensemble .

The thermodynamical stability analysis presented in this section also shows that the equilibrium states for $`D2`$ and $`D10`$ are always stable since the series of equilibria do not present turning points of energy or temperature.

## 3 Dynamics of self-gravitating Brownian particles in dimension $`D`$

### 3.1 The Smoluchowski-Poisson system

We now consider the dynamics of a system of self-gravitating Brownian particles in a space of dimension $`D`$. Like in Ref. , we consider a high friction limit in order to simplify the problem. We thus study the Smoluchowski equation

(77)
$$\frac{\rho }{t}=\left[\frac{1}{\xi }(T\rho +\rho \mathrm{\Phi })\right],$$

coupled to the Newton-Poisson equation (4). In the microcanonical ensemble, the temperature $`T(t)`$ evolves in time so as to satisfy the energy constraint (8). In the canonical ensemble, the temperature is constant. The Smoluchowski equation can be obtained from a variational principle called the Maximum Entropy Production Principle . This variational approach is interesting as it makes a direct link between the dynamics and the thermodynamics. In the microcanonical description, the rate of entropy production can be put in the form

(78)
$$\dot{S}=\frac{1}{T\rho \xi }(T\rho +\rho \mathrm{\Phi })^2d^D𝐫0.$$

For a stationary solution, $`\dot{S}=0`$ and we obtain the Boltzmann distribution (10) which is a critical point of entropy. Considering a small perturbation around equilibrium, we can establish the identity :

(79)
$$\delta ^2\dot{S}=2\lambda \delta ^2S0,$$

where $`\lambda `$ is the growth rate of the perturbation defined such that $`\delta \rho e^{\lambda t}`$. This relation shows that a stationary solution of the Smoluchowski-Poisson system is dynamically stable if and only if it is a local entropy maximum. We get similar results in the canonical ensemble with $`J`$ in place of $`S`$. The relation (79) has been found for other kinetic equations satisfying a $`H`$-theorem (Chavanis, in preparation). In $`D=2`$, an equation morphologically similar to the Smoluchowski equation (77) can been derived from the $`N`$-body Liouville equation for a gas of point vortices .

### 3.2 Self-similar solutions of the Smoluchowski-Poisson system

From now on, we set $`M=R=G=\xi =1`$. The equations of the problem become

(80)
$$\frac{\rho }{t}=(T\rho +\rho \mathrm{\Phi }),$$

(81)
$$\mathrm{\Delta }\mathrm{\Phi }=S_D\rho ,$$

(82)
$$E=\frac{D}{2}T+\frac{1}{2}\rho \mathrm{\Phi }d^D𝐫,$$

with boundary conditions

(83)
$$\frac{\mathrm{\Phi }}{r}(0,t)=0,\mathrm{\Phi }(1)=\frac{1}{2D},T\frac{\rho }{r}(1)+\rho (1)=0,$$

for $`D>2`$. For $`D=2`$, we take $`\mathrm{\Phi }(1)=0`$ on the boundary. Integrating Eq. (81) once, we can rewrite the Smoluchowski-Poisson system in the form of a single integrodifferential equation

(84)
$$\frac{\rho }{t}=\frac{1}{r^{D1}}\frac{}{r}\left\{r^{D1}\left(T\frac{\rho }{r}+\frac{\rho }{r^{D1}}_0^r\rho (r^{})S_Dr^{{}_{}{}^{}D1}𝑑r^{}\right)\right\}.$$

The Smoluchowski-Poisson system is also equivalent to a single differential equation

(85)
$$\frac{M}{t}=T\left(\frac{^2M}{r^2}\frac{D1}{r}\frac{M}{r}\right)+\frac{1}{r^{(D1)}}M\frac{M}{r},$$

for the quantity

(86)
$$M(r,t)=_0^r\rho (r^{})S_Dr^{{}_{}{}^{}D1}𝑑r^{},$$

which represents the mass contained within the sphere of radius $`r`$. The appropriate boundary conditions are

(87)
$$M(0,t)=0,M(1,t)=1.$$

It is also convenient to introduce the function $`s(r,t)=M(r,t)/r^D`$ satisfying

(88)
$$\frac{s}{t}=T\left(\frac{^2s}{r^2}+\frac{D+1}{r}\frac{s}{r}\right)+\left(r\frac{s}{r}+Ds\right)s.$$

We look for self-similar solutions of the form

(89)
$$\rho (r,t)=\rho _0(t)f\left(\frac{r}{r_0(t)}\right),r_0=\left(\frac{T}{\rho _0}\right)^{1/2}.$$

In terms of the mass profile, we have

(90)
$$M(r,t)=M_0(t)g\left(\frac{r}{r_0(t)}\right),\mathrm{with}M_0(t)=\rho _0r_0^D,$$

and

(91)
$$g(x)=_0^xf(x^{})S_Dx^{{}_{}{}^{}D1}𝑑x^{}.$$

In terms of the function $`s`$, we have

(92)
$$s(r,t)=\rho _0(t)S\left(\frac{r}{r_0(t)}\right),\mathrm{with}S(x)=\frac{g(x)}{x^D}.$$

Substituting the ansatz (92) into Eq. (88), we find that

(93)
$$\frac{d\rho _0}{dt}S(x)\frac{\rho _0}{r_0}\frac{dr_0}{dt}xS^{}(x)=\rho _0^2\left(S^{\prime \prime }(x)+\frac{D+1}{x}S^{}(x)+xS(x)S^{}(x)+DS(x)^2\right),$$

where we have set $`x=r/r_0`$. The variables of position and time separate provided that there exists $`\alpha `$ such that

(94)
$$\rho _0r_0^\alpha =\kappa ,$$

where $`\kappa `$ is a constant. In that case, Eq. (93) reduces to

(95)
$$\frac{d\rho _0}{dt}\left(S(x)+\frac{1}{\alpha }xS^{}(x)\right)=\rho _0^2\left(S^{\prime \prime }(x)+\frac{D+1}{x}S^{}(x)+xS(x)S^{}(x)+DS(x)^2\right).$$

Assuming that such a scaling exists implies that $`(1/\rho _0^2)(d\rho _0/dt)`$ is a constant that we arbitrarily set equal to $`\alpha `$ (note that this convention is different from the one adopted in Ref. ). This leads to

(96)
$$\rho _0(t)=\frac{1}{\alpha }(t_{coll}t)^1,$$

so that the central density becomes infinite in a finite time $`t_{coll}`$. The scaling equation now reads

(97)
$$\alpha S+xS^{}=S^{\prime \prime }+\frac{D+1}{x}S^{}+S(xS^{}+DS).$$

For $`x+\mathrm{}`$, we have asymptotically

(98)
$$S(x)x^\alpha ,g(x)x^{D\alpha },f(x)x^\alpha .$$

### 3.3 Canonical ensemble

In the canonical ensemble in which the temperature $`T`$ is a constant, we have <sup>1</sup><sup>1</sup>1The case $`T=0`$ is treated in the Appendix B

(99)
$$\alpha =2,\kappa =T.$$

In that case, the scaling equation (97) can be solved analytically. Following a procedure similar to the one developed in Ref. , we find that

(100)
$$S(x)=\frac{4}{D2+x^2}.$$

Then, Eqs. (92) and (91) yield

(101)
$$g(x)=\frac{4x^D}{D2+x^2},f(x)=\frac{4(D2)}{S_D}\frac{D+x^2}{\left(D2+x^2\right)^2}.$$

According to Eqs. (89) and (96), the central density evolves with time like

(102)
$$\rho (0,t)=\rho _0(t)f(0)=\frac{2D}{(D2)S_D}(t_{coll}t)^1.$$

According to Eqs. (89) and (90), the core radius and the core mass evolve like

(103)
$$r_0(t)=\sqrt{2T}(t_{coll}t)^{1/2},M_0(t)=\frac{1}{2}(2T)^{D/2}(t_{coll}t)^{\frac{D}{2}1}.$$

Note that for $`D>2`$, the core mass goes to zero at the collapse time. At $`t=t_{coll}`$, we get the singular profile

(104)
$$\rho (r,t=t_{coll})=\frac{4T(D2)}{S_Dr^2},M(r,t=t_{coll})=4Tr^{D2}.$$

### 3.4 Microcanonical ensemble

In the microcanonical ensemble, the exponent $`\alpha `$ is not determined by simple dimensional analysis. In Ref. , we found numerically that the scaling equation (97) has physical solutions only for $`\alpha \alpha _{\mathrm{max}}`$, with $`\alpha _{max}2.21`$ for $`D=3`$. We also argued that the system will select the exponent $`\alpha _{max}`$, since it leads to a maximum increase of entropy. In this section, we show that in the limit of large dimension $`D`$, we can explicitly understand the occurrence of such a $`\alpha _{\mathrm{max}}`$. In addition, we will present the derivation of perturbative expansions for $`\alpha _{\mathrm{max}}`$ and the scaling function $`S`$, in powers of $`D^1`$.

Eq. (97) can be formally integrated as a first order differential equation (writing $`S^{\prime \prime }=S^{}\times [S^{\prime \prime }/S^{}]`$), leading to an expression of $`S(x)`$ as a function of $`x`$, $`S(x)`$ itself, and $`S^{\prime \prime }(x)/S^{}(x)`$:

(105)
$$\left|\frac{\alpha }{DS(x)}1\right|=\left|\frac{\alpha }{DS(0)}1\right|\mathrm{exp}\left[\alpha _0^x\frac{ydy}{y^2(1S(y))y\frac{S^{\prime \prime }(y)}{S^{}(y)}(D+1)}\right].$$

We now define $`x_0`$, such that $`S(x_0)=\alpha /D`$. Since $`S`$ should be analytic, the foregoing relation implies for $`xx_0`$,

(106)
$$_0^x\frac{\alpha y}{F(y)}𝑑y\mathrm{ln}|xx_0|,$$

where $`F(y)`$ is the function which occurs in the denominator of the integral in Eq. (105). From Eq. (106), we must have $`F(y)=\alpha x_0(yx_0)`$ for $`yx_0`$, which implies $`F(x_0)=0`$ and $`F^{}(x_0)=\alpha x_0`$. These conditions can be rewritten

(107)
$$x_0^2\left(1\frac{\alpha }{D}\right)x_0\frac{S^{\prime \prime }(x_0)}{S^{}(x_0)}(D+1)=0,$$

(108)
$$(\alpha 2)x_0=\frac{d}{dx}\left[x^2S(x)+x\frac{S^{\prime \prime }(x)}{S^{}(x)}\right]_{x_0}.$$

This preparatory work now allows the introduction of a systematic expansion in large dimension $`D`$ for the scaling function $`S`$, the scaling exponent $`\alpha `$, and $`x_0`$. In this limit, let us neglect the contribution of the terms which are not of order $`D`$ in the right-hand side of Eq. (97). This actually amounts to take $`F(y)y^2D`$ in Eq. (105). Within this approximation, we find

(109)
$$\left|\frac{\alpha }{DS(x)}1\right|=\left|\frac{\alpha }{DS(0)}1\right|\times \left|\frac{x^2}{D}1\right|^{\alpha /2},$$

which is an analytic function only if $`\alpha =2`$. This leads to $`x_0=\sqrt{D}`$, and to the more explicit form for $`S`$,

(110)
$$S(x)=\frac{S(0)}{1+\left(\frac{DS(0)}{2}1\right)\frac{x^2}{D}}.$$

$`S(0)`$ remains undetermined, and will be fixed by the next order approximation. Indeed, we can iteratively solve the full scaling equation Eq. (97), by reinserting the zeroth order solution into Eq. (105), and eventually continue this process with the new improved scaling function, and so forth… Thus, expressing the conditions of Eq. (107) and Eq. (108), and defining $`z=\frac{DS(0)}{2}`$ (which will be of order $`O(1)`$), we obtain

(111)
$$x_0^2=D+\frac{4}{z}+O(D^1),\mathrm{or}x_0=\sqrt{D}\left(1+\frac{2}{zD}+O(D^2)\right),$$

and

(112)
$$\alpha 2=\frac{4}{D}\left[\frac{1}{z}\frac{2}{z^2}\right]+O(D^2).$$

Eq. (112) provides a relation between the possible values for $`\alpha `$ and the associated value of $`S(0)=\frac{2z}{D}`$. Note that the function of $`z`$ in the right-hand side of Eq. (112) has a well defined maximum. Hence, up to order $`O(D^1)`$, we find that $`\alpha \alpha _{\mathrm{max}}`$, with

(113)
$$\alpha _{\mathrm{max}}=2+\frac{1}{2}D^1+O(D^2),$$

which is associated to the value $`z=4+O(D^1)`$ or $`S(0)=\frac{8}{D}+O(D^2)`$. As $`\alpha `$ is necessarily greater than 2 (as the temperature cannot vanish), a solution exists for any $`\alpha [2,\alpha _{\mathrm{max}}]`$. As already mentioned, $`\alpha _{\mathrm{max}}`$ is dynamically selected as it leads to the maximum divergence of the entropy and the temperature (see Eq. (123) below).

Inserting Eq. (110) into Eq. (105), we find the next order approximation for $`S`$

(114)
$$\left|\frac{\alpha }{DS(x)}1\right|=\left|\frac{\alpha }{DS(0)}1\right|\times \left|\frac{x^2}{x_0^2}1\right|^{\frac{\alpha }{2}(1\varphi )}\times \left[\frac{x^2}{x_1^2}+1\right]^{\frac{\alpha \varphi }{2}},$$

where $`x_0`$ is given by Eq. (111), and $`x_1`$ and $`\varphi `$ are defined by

(115)
$$x_1^2=\frac{D}{z1}+\frac{2(z2)}{z(z1)}+O(D^1),\varphi =\frac{2}{D}\left[\frac{1}{z}\frac{2}{z^2}\right]+O(D^2).$$

Again, the analyticity condition imposes that $`\frac{\alpha }{2}(1\varphi )=1`$, which exactly leads to Eq. (112), and to the following explicit form for $`S`$:

(116)
$$S(x)=\frac{\alpha }{D}\left[1+\left(1\frac{\alpha }{2z}\right)\left(\frac{x^2}{x_0^2}1\right)\left(\frac{x^2}{x_1^2}+1\right)^{\frac{\alpha }{2}1}\right]^1.$$

This improved scaling function can be inserted again into the conditions expressed by Eq. (107) and Eq. (108), leading to the next order term in the expansion of $`\alpha `$. After elementary but cumbersome calculations, we end up with

(117)
$$\alpha 2=\frac{4}{D}\left[\frac{1}{z}\frac{2}{z^2}\right]+\frac{8}{D^2}\left[\frac{5}{z}\frac{26}{z^2}+\frac{31}{z^3}\frac{6}{z^4}\left(\frac{1}{z}\frac{7}{z^2}+\frac{14}{z^3}\frac{8}{z^4}\right)\mathrm{ln}z\right]+O(D^3).$$

This function has again a well defined maximum for

(118)
$$z=\frac{D}{2}S(0)=4+\left(\frac{41}{2}6\mathrm{ln}2\right)D^1+O(D^2),$$

associated to the value

(119)
$$\alpha _{\mathrm{max}}=2+\frac{1}{2}D^1+\frac{11}{16}D^2+O(D^3).$$

This expansion gives $`\alpha _{\mathrm{max}}=2.24\mathrm{}`$ in $`D=3`$, in fair agreement with the exact value $`\alpha _{\mathrm{max}}=2.2097\mathrm{}`$ obtained numerically in . In addition, the exponent $`\alpha =2`$ is associated to $`z=2+4D^1+O(D^2)`$. In principle, these expansions can be systematically pursued to the prize of increasingly complicated calculations.

Finally note that Eqs. (89) and (96) lead to the following exact asymptotic for the central density $`\rho (0,t)`$:

(120)
$$\rho (0,t)K_D(\alpha )(t_{coll}t)^1,K_D(\alpha )=\frac{2z(\alpha )}{\alpha S_D},$$

where we have used $`f(0)=DS(0)/S_D`$ and the definition of $`z`$. The function $`z(\alpha )`$ is determined implicitly by Eq. (117), up to order $`O(D^2)`$. For the special cases $`\alpha =2`$ and $`\alpha =\alpha _{\mathrm{max}}`$, we respectively find

(121) $`K_D(2)`$ $`=`$ $`2S_D^1\left(1+2D^1+O(D^2)\right),`$
(122) $`K_D(\alpha _{\mathrm{max}})`$ $`=`$ $`4S_D^1\left(1+\left({\displaystyle \frac{39}{8}}{\displaystyle \frac{3}{2}}\mathrm{ln}2\right)D^1+O(D^2)\right),`$
which shows that $`K_D(\alpha _{\mathrm{max}})`$ is substantially greater than $`K_D(2)`$ (twice bigger in the infinite $`D`$ limit, the ratio being even bigger for finite $`D`$, as $`\frac{39}{8}\frac{3}{2}\mathrm{ln}23.835\mathrm{}>2`$). This substantial difference was noted in Ref. , in the case $`D=3`$. Finally, as expected in the microcanonical ensemble, the temperature diverges during the collapse as $`T(t)(t_{coll}t)^{\alpha _T}`$ with $`\alpha _T=12/\alpha `$; see Eqs. (89)(94)(96). The strongest divergence is obtained for $`\alpha =\alpha _{\mathrm{max}}`$. According to Eq. (119), we have

(123)
$$\alpha _T=2\frac{2}{\alpha _{\mathrm{max}}}=\frac{1}{4}D^1+\frac{9}{32}D^2+O(D^3).$$

If we plug $`D=3`$ in Eq. (123), we find the estimate $`\alpha _T0.11\mathrm{}`$ in fair agreement with the exponent measured numerically in , $`\alpha _T0.1`$.

## 4 The two-dimensional case

### 4.1 The critical temperature

In two dimensions, the dynamical equation (85) for the mass profile reads

(124)
$$\frac{M}{t}=4Tu\frac{^2M}{u^2}+2M\frac{M}{u},$$

after the change of variable $`u=r^2`$ has been effected. Looking for a stationary solution, and using $`uM^{\prime \prime }=(uM^{})^{}M^{}`$, Eq. (124) is readily integrated leading to

(125)
$$M(u)=\frac{4T}{4T1}\frac{u}{1+\frac{u}{4T1}},$$

Note that $`M(1)=1`$, ensuring that the whole mass is included in this solution. Using $`\rho =M^{}/\pi `$, we find that the density profile is given by

(126)
$$\rho (r)=\frac{4\rho _0}{\pi }\frac{1}{(1+(r/r_0)^2)^2},$$

with

(127)
$$r_0=\sqrt{4T1}\mathrm{and}\rho _0r_0^2=T.$$

This solution exists provided that $`T>T_c=1/4`$, which defines the collapse temperature. We have thus recovered the result (40) by a slightly different method. Note that the value of $`T_c`$ and the dependance of $`r_0`$ and $`\rho _0`$ with the temperature coincide with the exact results obtained within conformal field theory . In the following, $`T`$ is set constant as we have already seen that the gravothermal does not exist in the microcanonical ensemble in two dimensions.

### 4.2 Scaling collapse for $`T=T_c`$

We now address the dynamics at the critical temperature $`T=T_c=1/4`$. We note that contrary to what happens in other dimensions, the central density diverges at $`T_c`$. Thus, in analogy with critical phenomena, we anticipate a scaling solution for $`M(u,t)`$, of the form

(128)
$$M(u,t)\frac{(a(t)+1)u}{1+a(t)u},$$

which preserves the scaling form obtained above $`T_c`$, and which satisfies the boundary condition $`M(1,t)=1`$. The corresponding density profile is

(129)
$$\rho (r,t)=\frac{a(t)+1}{\pi }\frac{1}{(1+a(t)r^2)^2}.$$

The central density

(130)
$$\rho (0,t)=\frac{a(t)+1}{\pi },$$

is expected to diverge for $`t+\mathrm{}`$, so that $`a(t)`$ is also expected to diverge.

Inserting the ansatz Eq. (128) into Eq. (124) shows that the left-hand term is indeed negligible compared to both terms of the right-hand side, to leading order in $`a`$. So far, this prevents us from determining a dynamical equation for $`a`$. In order to achieve that, we must solve Eq. (124) to the next order in $`a^1`$. We thus look for a solution of the form

(131)
$$M(u,t)=\frac{a(t)u}{1+a(t)u}+a(t)^1h(u,t),$$

where $`h(u,t)`$ is expected to be of order $`O(1)`$, and should satisfy $`h(0,t)=0`$ and $`h(1,t)=1`$ (the total integrated mass should be 0 and 1, respectively, for $`u=0`$ and $`u=1`$), and $`\frac{h}{u}(0,t)=0`$, which ensures that Eq. (130) is exactly obeyed, defining $`a(t)`$ without any ambiguity. The contribution of $`\frac{M}{t}`$ in the left-hand side of Eq. (124) is dominated by the time derivative of Eq. (128) :

(132)
$$\frac{}{t}\left[\frac{(1+a(t))u}{1+a(t)u}\right]=\frac{u(1u)}{(1+au)^2}\frac{da}{dt},$$

which will be checked self-consistently hereafter. In addition, non linear terms in $`h`$ in the right-hand side are also negligible. Therefore, $`h`$ satisfies

(133)
$$\frac{au(1u)}{(1+au)^2}\frac{da}{dt}=u\frac{^2h}{u^2}+2\frac{}{u}\left(\frac{au}{1+au}h\right).$$

Actually, for a given time, this equation becomes an ordinary differential equation involving only one variable $`u`$, as $`a`$ and $`da/dt`$ appear as parameters. Eq. (133) can be integrated leading to a first order equation in $`h`$, which can be easily solved. Defining $`v=au`$, we finally get

(134) $`h(u,t)=a^1(1+{\displaystyle \frac{2}{a}}){\displaystyle \frac{da}{dt}}(1+v)^2\times `$
$`\left[(v^21)\mathrm{ln}(1+v)+v(12v)+2v{\displaystyle _0^v}{\displaystyle \frac{\mathrm{ln}(1+z)}{z}}𝑑z{\displaystyle \frac{2v^2+v^3}{2(a+2)}}\right],`$
which depends on time only through the variable $`a`$ and $`da/dt`$. Now, $`da/dt`$ is determined by imposing the boundary condition $`h(1,t)=1`$, which leads to

(135)
$$\frac{da}{dt}=\frac{a}{\mathrm{ln}a5/2}\left[1+O(\mathrm{ln}a^2)\right].$$

One can solve iteratively Eq. (124), by adding the time derivative of the above solution to the left-hand side, in order to compute an improved $`h`$. To leading order, the solution of Eq. (134) is preserved. However, new terms are generated which are important for $`va`$ ($`u1`$), and which generate terms of order $`O(a/\mathrm{ln}^3a)`$ is the expansion for $`da/dt`$. This explains the form of the error term in Eq. (135).

Integrating Eq. (135) for large time, we get the exact asymptotic expansion for large time

(136)
$$a(t)=\mathrm{exp}\left(\frac{5}{2}+\sqrt{2t}\right)\left[1+O(t^{1/2}\mathrm{ln}t)\right].$$

For $`t+\mathrm{}`$, the central density diverges like $`a(t)`$ and the core radius goes to zero like $`a(t)^{1/2}`$. In addition, the scaling solution (129) at $`T=T_c`$ goes to a Dirac peak containing the whole mass (see Eq. (41)), as the decay exponent of the scaling function is 4, which is strictly greater than 2.

### 4.3 Collapse for $`T<T_c`$

For $`D=2`$, the scaling equation associated to Eq. (85) does not display any physical solution when numerically solved. In this section, we thus present a special treatment adapted to this case. The principal difference with other dimensions is the divergence of the central density at $`T_c`$, and the occurrence of a scaling solution at this temperature.

Strictly below $`T_c`$, we expect a finite time collapse. Close to the center, the solution takes the form

(137)
$$M(u,t)4T\frac{a(t)u}{1+a(t)u},$$

where again the left-hand side of Eq. (124) is negligible compared to each term on the right-hand side. We thus look for a solution of the type

(138)
$$M(u,t)=4T\frac{a(t)u}{1+a(t)u}+h(u,t),$$

where $`h`$ is now of order $`O(1)`$ as it contains a finite fraction of the total mass, since the first term contains a mass of order $`4T<1`$. Inserting this ansatz in the dynamical equation Eq. (124), we obtain

(139)
$$\frac{1}{4T}\frac{h}{t}+\frac{da}{dt}\frac{u}{(1+au)^2}=u\frac{^2h}{u^2}+2\frac{}{u}\left(\frac{au}{1+au}h\right)+2h\frac{h}{u}.$$

One can look for a scaling solution of the type

(140)
$$h(u,t)=a^{\gamma 1}H(au),\mathrm{with}H(v)cv^{1\gamma },\mathrm{when}v+\mathrm{},$$

so that the mass included in this scaling profile $`h(1,t)=c=O(1)`$. With this definition, the density profile decays for large distance as $`\rho r^\alpha `$, with $`\alpha =2\gamma `$. Inserting this ansatz in Eq. (139), we obtain

(141)
$$\left[\frac{1}{4T}\left(vH^{}+(\gamma 1)H\right)+a^{1\gamma }\frac{v}{(1+v)^2}\right]\frac{da}{dt}a^2=vH^{\prime \prime }+2\left(\frac{v}{1+v}H\right)^{}+2a^{\gamma 1}HH^{},$$

where derivatives are with respect to the variable $`v`$. We are free to choose $`a(t)=\pi \rho (0,t)/(4T)`$, so that $`H^{}(0)=H(0)=0`$. For small $`v`$, Eq. (141) leads to

(142)
$$\frac{da}{dt}=H^{\prime \prime }(0)a^{\gamma +1}.$$

Eq. (141) has a global scaling solution only for $`\gamma =1`$. However, we know that in this case the scaling equation obtained by setting $`\gamma =1`$ does not display any physical solution. Thus, we conclude that there is no scaling solution obtained by imposing that all terms in Eq. (141) scale the same way. However, as we will see in the section devoted to numerical simulations, the direct simulation of Eq. (124) seems to display a scaling solution with $`\gamma 0.60.7`$ for numerically accessible densities. Strictly speaking, this is totally excluded by the above equation, except if one allows $`\gamma `$ to very slowly depend on the density or $`a`$. For a given $`a`$, we thus want to solve Eq. (141), where the boundary conditions will ultimately select the effective value of $`\gamma `$, and that of $`da/dt`$. More precisely, once we impose $`H^{}(0)=H(0)=0`$, and the condition of Eq. (142), we end up with a shooting problem for $`H^{\prime \prime }(0)`$ and $`\gamma `$. For large $`a`$, and $`va`$, it is clear that the non linear term of the right-hand side of Eq. (141) becomes irrelevant, and we drop it from now on.

In order to understand the origin of this shooting problem, and to obtain an accurate estimate of $`\gamma `$, let us solve Eq. (141) in the limit of very large $`a`$, in the range $`1va`$. In this regime, Eq. (141) simplifies to the following equation

(143)
$$\left[\frac{1}{4T}\left(vH^{}+(\gamma 1)H\right)+a^{1\gamma }v^1\right]\omega =vH^{\prime \prime }+2H^{},$$

where

(144)
$$\omega =\frac{da}{dt}a^2=H^{\prime \prime }(0)a^{\gamma 1}.$$

Let us now multiply this equation by $`v^{\gamma 2}`$ and integrate the resulting equation. After elementary manipulations, we obtain

(145)
$$H^{}+\left[\frac{3\gamma }{v}\frac{\omega }{4T}\right]H=\frac{\omega c}{4T}v^{1\gamma }\frac{\omega a^{1\gamma }}{2\gamma }v^1+(2\gamma )(3\gamma )v^{1\gamma }_v^+\mathrm{}w^{\gamma 3}H(w)𝑑w,$$

where $`cO(1)`$, which has been defined in Eq. (140), appears here as an integration constant. Then, one can integrate this differential equation which leads to the following self-consistent relation for $`H`$:

(146)
$$H(v)=v^{\gamma 3}\mathrm{exp}\left(\frac{\omega v}{4T}\right)_v^+\mathrm{}w^{3\gamma }\mathrm{exp}\left(\frac{\omega w}{4T}\right)F(w)𝑑w,$$

where $`F`$ is defined as the opposite of the right-hand side of Eq. (145):

(147)
$$F(v)=\frac{\omega c}{4T}v^{1\gamma }+\frac{\omega a^{1\gamma }}{2\gamma }v^1(2\gamma )(3\gamma )v^{1\gamma }_v^+\mathrm{}w^{\gamma 3}H(w)𝑑w.$$

Eq. (143) implies that $`H(v)\mathrm{ln}v`$, when $`v0`$ (of course, this apparent divergence does not occur in the full dynamical equation Eq. (141)). Considering the prefactor $`v^{\gamma 3}`$ in Eq. (146), this behavior can be obtained if and only if

(148)
$$_0^+\mathrm{}w^{3\gamma }\mathrm{exp}\left(\frac{\omega w}{4T}\right)F(w)𝑑w=0.$$

As $`\omega `$ is expected to go to zero for large $`a`$ as $`\gamma <1`$ (see Eq. (144)), the dominant contribution of the integral of the third term in the definition of $`F`$ comes from the large $`w`$ region, for which $`H`$ can be replaced by its asymptotic form. Hence, defining $`\mathrm{\Gamma }(x)=_0^+\mathrm{}w^x\mathrm{exp}(w)𝑑w`$ and $`\epsilon =1\gamma `$, and using Eq. (144), the condition expressed in Eq. (148) can be rewritten

(149)
$$c=\frac{\mathrm{\Gamma }(1+\epsilon )}{\epsilon (1+\epsilon )^2\mathrm{\Gamma }(1+2\epsilon )}H^{\prime \prime }(0)\left(\frac{H^{\prime \prime }(0)}{4T}\right)^\epsilon a^{\epsilon ^2}.$$

As $`c`$ is of order $`O(1)`$, we find that $`\epsilon 0`$ as $`a+\mathrm{}`$. More precisely, in this limit, $`\epsilon `$ is the solution of the following implicit equation

(150)
$$\epsilon =\sqrt{\frac{\mathrm{ln}(K/\epsilon )}{\mathrm{ln}a}},$$

where $`K=H^{\prime \prime }(0)/c+O(\epsilon )`$. Finally, we obtain

(151)
$$\epsilon =1\gamma =\sqrt{\frac{\mathrm{ln}\mathrm{ln}a}{2\mathrm{ln}a}}\left(1+O([\mathrm{ln}\mathrm{ln}a]^1)\right).$$

In conclusion, although the solution is not strictly speaking a true scaling solution, the explicit dependence of $`\gamma `$ on $`a`$ is so weak that an apparent scaling should be seen with an effective $`\gamma `$ almost constant for a wide range of values of $`a`$. Hence, the total density profile is the sum of the scaling profile obtained at $`T_c`$ with a $`T/T_c`$ weight (behaving as a Dirac peak of weight $`T/T_c`$, at $`t=t_{coll}`$) and of a pseudo-scaling solution associated to an effective scaling exponent slowly converging to $`\alpha =2`$.

Let us illustrate quantitatively the time dependence of $`\alpha =2\gamma `$. For example, taking arbitrarily $`K=1`$ (the dependence on $`K`$ is weak and vanishes for large $`a`$), Eq. (150) and Eq. (151) respectively lead to $`\gamma (a=10^3)=0.624\mathrm{}`$ and $`\gamma (a=10^3)=0.626\mathrm{}`$, and to $`\gamma (a=10^5)=0.684\mathrm{}`$ and $`\gamma (a=10^5)=0.674\mathrm{}`$ (note that the error between the asymptotic expansion of Eq. (151) and the implicit expression first grows before slowly decaying for $`a10^{12}`$ !). Finally, for the maximum value of $`a`$ accessible numerically of order $`a10^4`$, we expect to observe an apparent scaling solution with $`\gamma 0.65`$, or $`\alpha =2\gamma 1.3`$.

### 4.4 Numerical simulations

In this subsection, we present numerical simulations concerning the two-dimensional case. Indeed, the three-dimensional case has been extensively studied in . It has been shown that the scaling function as well as the corrections to scaling (which have been calculated for the canonical ensemble in ) are perfectly described by the theory. In addition, as the system behaves qualitatively the same for any dimension $`D>2`$ in both thermodynamical ensembles, we have decided to focus on the numerical study of the $`D=2`$ case, which displays some very rich behaviors, as exemplified in the present section.

We consider the system in the canonical ensemble, as the gravitational collapse does not occur in the microcanonical ensemble. In Fig. 9, we show the scaling function at $`T_c`$, as given by Eq. (129), finding a perfect agreement. In Fig. 10, we also display $`a(da/dt)^1`$ as a function of $`\mathrm{ln}a`$, and find an asymptotic behavior fully compatible with that given by Eq. (135).

Below $`T_c`$, and in the accessible range of $`a`$ (up to $`a10^5`$), we find an apparent scaling regime with $`\alpha =2\gamma 1.3`$, as predicted in Sec. 4.3. This is illustrated in Fig. 11, for $`T=T_c/2=1/8`$. Note that the effective $`\gamma `$ or $`\alpha `$ can also be extracted from the time evolution of $`a(t)`$ or the central density (see Eq. (142)). In Fig. 12, we show that this way of measuring $`\gamma `$ is fully compatible with the value of the effective exponents $`\alpha =2\gamma 1.3`$

## 5 The one dimensional case

When an equilibrium state exists, there is little hope to be able to solve the full Smoluchowski-Poisson system analytically in order to study the relaxation towards equilibrium. We shall consider a simpler problem in which a test particle evolves in a medium of field particles at statistical equilibrium (thermal bath approximation). The particles are assumed to create a stationary potential $`\mathrm{\Phi }_{eq}(𝐫)`$ which induces a drift of the test particle along the gradient of $`\mathrm{\Phi }_{eq}`$. In addition, the test particle is assumed to experience a diffusion process. If $`\rho `$ denotes the density probability of finding the test particle in $`𝐫`$ at time $`t`$, we expect the evolution of $`\rho `$ to be determined by a Smoluchowski equation of the form

(152)
$$\frac{\rho }{t}=(T\rho +\rho \mathrm{\Phi }_{eq}),$$

where $`\mathrm{\Phi }_{eq}(𝐫)`$ is solution of the Boltzmann-Poisson equation (11). This means that we replace the true potential by its equilibrium value but still allow the density $`\rho `$ to vary with time. As we shall see, it is possible to solve the Smoluchowski equation (152) analytically in $`D=1`$ by using an analogy with a problem of quantum mechanics. An equation of the form (152) has been proposed in Refs.  to model the motion of a test vortex in a bath of field vortices at statistical equilibrium. In that context, Eq. (152) can be formally derived from the $`N`$-body Liouville equation of the system by using projection operator technics.

It is well-known that a Fokker-Planck equation like Eq. (152) can be formally transformed into a Schrödinger equation with imaginary time. Indeed, performing the change of variable

(153)
$$\rho =\psi e^{\frac{1}{2T}\mathrm{\Phi }_{eq}},$$

we find that the evolution of $`\psi `$ is determined by an equation of the form

(154)
$$\frac{\psi }{t}=T\mathrm{\Delta }\psi +\left(\frac{1}{2}\mathrm{\Delta }\mathrm{\Phi }_{eq}\frac{1}{4T}(\mathrm{\Phi }_{eq})^2\right)\psi .$$

This can be written as a Schrödinger-type equation

(155)
$$\frac{\psi }{t}=T\mathrm{\Delta }\psi V(𝐫)\psi ,$$

with a potential $`V(𝐫)=\frac{1}{2}\mathrm{\Delta }\mathrm{\Phi }_{eq}+\frac{1}{4T}(\mathrm{\Phi }_{eq})^2`$. So far, this transformation is general. If we now consider the one-dimensional case, the Boltzmann-Poisson equation (152) can be solved analytically and the potential $`V(r)`$ determined explicitly. Introducing the notations $`\xi =\alpha r/\sqrt{2}R`$ and $`\tau =\alpha ^2Tt/2R^2`$ and using Eq. (24) we can rewrite Eq. (154) in the form

(156)
$$\frac{\psi }{\tau }=\mathrm{\Delta }\psi +\left(\frac{2}{\mathrm{cosh}^2\xi }1\right)\psi .$$

A separation of the variables can be effected by the substitution

(157)
$$\psi (\xi ,t)=\varphi (\xi )e^{\lambda t}(\lambda 0),$$

where $`\varphi `$ is solution of the ordinary differential equation

(158)
$$\frac{d^2\varphi }{d\xi ^2}+2\left(E+\frac{1}{\mathrm{cosh}^2\xi }\right)\varphi =0,$$

where we have set $`\lambda 1=2E`$. The solutions of this Schrödinger equation are described in detail in Ref. . The spectrum of positive energies is continuous. The spectrum of negative energies is discrete and reduces to $`E_0=1/2`$ (fundamental state). The first excited state in the continuum is $`E_1=0`$. We can check that the corresponding eigenfunctions are $`\varphi _0=1/\mathrm{cosh}\xi `$ and $`\varphi _1=\mathrm{tanh}\xi `$. In order to obtain the qualitative behavior of the time dependent solution of Eq. (152), we neglect the contribution from the continuum states with $`E>0`$, only keeping the $`E=1/2`$ and $`E=0`$ eigenstates.

Within this approximation and for sufficiently large times, we obtain

(159)
$$\psi (\xi ,\tau )=\frac{A}{\mathrm{cosh}\xi }+B\mathrm{tanh}\xi e^\tau ,(\tau +\mathrm{}).$$

where $`A`$ and $`B`$ are constant. Returning to original variables, we get

(160)
$$\rho (r,t)=\rho _{eq}(r)\left\{1+C\mathrm{sinh}\left(\frac{\alpha r}{\sqrt{2}R}\right)e^{\frac{\alpha ^2T}{2R^2}t}+\mathrm{}\right\},$$

where $`\rho _{eq}`$ is given by Eq. (43) and $`C=B/A`$ is a constant. We find that the relaxation time is given by $`t_{relax}=2R^2/\alpha ^2T`$.

## 6 Conclusion

In this paper, we have studied the Boltzmann-Poisson equation and the Smoluchowski-Poisson system in various dimensions of space. Our study shows in particular how the nature of the problem changes as we pass from $`D=3`$ to $`D=2`$. We showed that the dimension $`D=2`$ is critical in the sense that the results obtained for $`D>2`$ diverge if they are naively extrapolated to $`D=2`$. On a physical point of view, the two-dimensional problem differs from the $`D>2`$ case in two respects: in the 2D case, the central density of the equilibrium state is infinite at the critical temperature $`T_c`$ while it is finite at $`T_c`$ in higher dimensions. On the other hand, in $`D=2`$, the self-similar collapse results in a Dirac peak which contains a finite fraction of mass, while for $`D>2`$, the mass contained in the core tends to zero at the collapse time (but a Dirac peak is always formed in the canonical ensemble after $`t_{coll}`$ as discussed in Appendix B). We have also evidenced another characteristic dimension $`D=10`$ at which the nature of the problem changes. For $`D10`$ the classical spiral characterizing isothermal spheres in the physical three-dimensional space disappears. However, since the points in the spiral correspond to unstable states, that are therefore unphysical, this transition at $`D=10`$ is not so critical and indeed the nature of the self-similar collapse does not show any transition at that dimension. It is interesting to note that the dependance of the phase diagram in the $`(E,T)`$ and $`(u,v)`$ planes with the dimension of space $`D`$ shows some resemblance with the dependance of the phase diagram of confined polytropic spheres with the index $`n`$ of the polytrope . An extension of our study would be to relax the high friction limit and consider solutions of the Kramers-Poisson system and other relaxation equations described in Ref. . These equations are expected to display qualitatively similar behaviors than those described here (i.e. gravitational collapse, finite time singularity, self-similar solutions,…) but their study appears to be of considerable difficulty since we now need to consider the evolution of the system in phase space. We hope to come to that problem in future publications.

## Appendix A Absence of global maximum of free energy in $`D=2`$ for $`T<T_c=GM/4`$

We give a proof similar to the one given in Ref.  for $`D=3`$. In two dimensions, we consider a homogeneous disk of mass $`M`$ and radius $`a`$ at temperature $`T`$. It is easy to show that the total energy (8) of this disk is

(161)
$$E=MT+\frac{1}{2}GM^2\left(\mathrm{ln}a\frac{1}{4}\right),$$

with the convention $`\mathrm{\Phi }GM\mathrm{ln}r`$ at large distances. According to Eqs. (6) and (9), its free energy reads

(162)
$$J=M\mathrm{ln}TM\mathrm{ln}\left(\frac{M}{\pi a^2}\right)M\frac{GM^2}{2T}\left(\mathrm{ln}a\frac{1}{4}\right).$$

For $`a0`$, the free energy behaves like

(163)
$$J2M\left(1\frac{GM}{4T}\right)\mathrm{ln}a.$$

Therefore, if $`T<T_c=GM/4`$ the free energy goes to $`+\mathrm{}`$ when we contract the system to a point. This is sufficient to prove the absence of a global maximum of free energy below $`T_c`$. This also shows the natural tendency (in a thermodynamical sense) of a canonical self-gravitating system to collapse to a Dirac peak for $`T<T_c`$.

## Appendix B The case of cold systems ($`T=0`$)

For $`T=0`$, Eq. (88) reduces to

(164)
$$\frac{s}{t}=\left(r\frac{s}{r}+ds\right)s.$$

Looking for a self-similar solution of the form (92) and imposing the conditions (94) and (96), we find that the scaling profile satisfies

(165)
$$xS^{}+\alpha S=\left(xS^{}+DS\right)S.$$

Of course, for $`T=0`$, the exponent $`\alpha `$ cannot be determined on dimensional grounds, as the definition $`r_0=\sqrt{T/\rho _0}`$ is not relevant anymore. As we will see, $`\alpha `$ will be determined by imposing that the scaling solution is analytic. Eq. (165) can be readily solved leading to the following implicit equation for $`S`$:

(166)
$$\left(\frac{\alpha }{D}S(x)\right)^{1\frac{\alpha }{D}}=Kx^\alpha S(x),$$

where $`K`$ is an integration constant. Now, from the definition of $`S`$, we expect a small $`x`$ expansion of the form $`S(x)=S(0)+\frac{1}{2}S^{\prime \prime }(0)x^2+O(x^4)`$, which first implies that

(167)
$$S(0)=\frac{\alpha }{D},$$

and that $`(x^2)^{1\alpha /D}x^\alpha `$, which finally leads to

(168)
$$\alpha =\frac{2D}{D+2}\mathrm{and}K=\frac{D+2}{2}\left(\frac{1}{2}|S^{\prime \prime }(0)|\right)^{\frac{D}{D+2}}.$$

In terms of the scaling function $`g(x)`$ associated to the mass profile, Eq. (166) can be rewritten

(169)
$$g(x)=\frac{2x^D}{D+2}\frac{|S^{\prime \prime }(0)|}{2}\left[\frac{D+2}{2}g(x)\right]^{\frac{D+2}{D}},$$

where $`S^{\prime \prime }(0)<0`$ is arbitrary. This leads to the exact large $`x`$ asymptotic behavior

(170)
$$g(x)\frac{2}{D+2}\left(\frac{4}{(D+2)|S^{\prime \prime }(0)|}\right)^{\frac{D}{D+2}}x^{\frac{D^2}{D+2}}.$$

Moreover, using Eq. (96) and (167), we get the exact universal asymptotic behavior of the central density

(171)
$$\rho (0,t)S_D^1(t_{coll}t)^1,$$

Finally, we note that the implicit equation (169) can be written as a parametric set of equations

(172)
$$g(y)=\frac{2}{D+2}y,x(y)=\left[y+\frac{D+2}{4}|S^{\prime \prime }(0)|y^{\frac{D+2}{D}}\right]^{\frac{1}{D}}.$$

These results can be obtained by a different, more physical, method. We have indicated in Ref.  that, for $`T=0`$, the particles have a deterministic motion with equation

(173)
$$\frac{d𝐫}{dt}=𝐮=\mathrm{\Phi }.$$

For a spherically symmetrical system, this can be rewritten

(174)
$$\frac{dr}{dt}=\frac{M(r,t)}{r^{D1}},$$

where $`M(r,t)`$ is the mass within $`r`$. If $`a`$ denotes the initial position of the particle located at $`r`$ at time $`t`$, we have

(175)
$$M(r,t)=M(a,0),$$

so Eq. (174) can be integrated explicitly in

(176)
$$r^D=a^DDM(a,0)t.$$

If $`M(a,0)`$ behaves like

(177)
$$M(a,0)=A(a^DBa^{D+2})+\mathrm{},$$

close to the origin (which is a regular expansion), then

(178)
$$M(r,t)=Aa^D(1Ba^2),\mathrm{with}r^D=(1DAt)a^D+DABa^{D+2}t.$$

Introducing the collapse time $`t_{coll}=1/DA`$ and considering the limit $`tt_{coll}`$, we obtain

(179)
$$M(r,t)=\frac{a^D}{Dt_{coll}},\mathrm{with}r^D=\frac{1}{t_{coll}}(t_{coll}t)a^D+Ba^{D+2}.$$

Introducing the scaling variables

(180)
$$x=\frac{r}{(t_{coll}t)^{\frac{D+2}{2D}}},y=\frac{1}{t_{coll}}\left[\frac{a}{(t_{coll}t)^{1/2}}\right]^D,$$

we can put the solution in a parametric form

(181)
$$M(r,t)=\frac{1}{D}(t_{coll}t)^{D/2}y,\mathrm{with}x=(y+Cy^{\frac{D+2}{D}})^{\frac{1}{D}},$$

where $`C`$ is a constant. At the collapse time $`t=t_{coll}`$,

(182)
$$M(r,t=t_{coll})=\frac{1}{DC^{\frac{D}{D+2}}}r^{\frac{D^2}{D+2}},\rho (r,t=t_{coll})=\frac{D}{(D+2)S_DC^{\frac{D}{D+2}}}r^{\frac{2D}{D+2}}.$$

These results are of course equivalent to those obtained previously.

We can now use this method to study the evolution of the system for $`t>t_{coll}`$ (post-collapse solution). For $`t=t_{coll}+\delta t`$, according to Eqs. (176) and (182), the mass contained inside the sphere of radius $`a_{coll}=C^{1/2}\delta t^{\frac{D+2}{2D}}`$ at $`t=t_{coll}`$ has collapsed at $`r=0`$, resulting in a Dirac peak of weight

(183)
$$M(0,t)=\frac{1}{DC^{D/2}}(tt_{coll})^{D/2}.$$

Note that in a bounded domain the final collapse to a central Dirac peak containing the whole mass occurs in a finite time after $`t_{coll}`$. For $`r>0`$ (associated to $`a>a_{coll}`$), one has

(184)
$$M(r,t)=M(0,t)+\frac{1}{DC^{\frac{D}{D+2}}}\left(a^{\frac{D^2}{D+2}}a_{coll}^{\frac{D^2}{D+2}}\right),$$

(185)
$$r^D=a^D\left(1\left(\frac{a_{coll}}{a}\right)^{\frac{2D}{D+2}}\right).$$

Introducing the scaling variables

(186)
$$x=\frac{r}{a_{coll}},y=\left(\frac{a}{a_{coll}}\right)^{\frac{D^2}{D+2}}1,$$

we obtain

(187)
$$M(r,t)=M(0,t)(1+y),\mathrm{with}x=(1+y)^{\frac{D^2}{D+2}}\left(1(1+y)^{\frac{2}{D}}\right)^{\frac{1}{D}}.$$

Subtracting the Dirac peak at $`r=0`$, and considering $`x1`$, for which $`y\frac{D}{2}x^D`$, we find that the leading contribution to the mass profile for small $`r`$ is

(188)
$$M(r,t)_>\frac{r^D}{2\delta t}.$$

Hence the density profile does not diverge at $`r=0^+`$ for $`t>t_{coll}`$. Instead, the density approaches the constant value

(189)
$$\rho (0+,t)=\frac{D}{2S_D\delta t},$$

which decreases with time. The density profile is thus depleted on a scale $`ra_{coll}\delta t^{\frac{D+2}{2D}}`$, which increases with time. For $`ra_{coll}`$, the density profile remains essentially unaffected.

In principle, the same phenomenon arises for any $`0<T<T_c`$: the density profile obtained at $`t_{coll}`$ ultimately collapses into a central Dirac peak at a time $`t_{}>t_{coll}`$. This solves the apparent paradox that the solution at $`t=t_{coll}`$ has a vanishing central mass and a finite free energy . In fact, if we allow singular profiles to develop, the evolution continues for $`t>t_{coll}`$ and the Dirac peak with infinite free energy (predicted by statistical mechanics ) is formed during the post collapse regime of our Brownian model <sup>2</sup><sup>2</sup>2As discussed in Sec. 2.1, the results should be different in the microcanonical ensemble. We shall reserve for a future communication the full description of the post-collapse regime.. In practice, degeneracy effects (of quantum or dynamical origin) lead to a finite small core of finite density, controlled by the maximum allowed degeneracy .