5 Results: I. Properties of equilibrium states in numerical and analytical models

Here, some N-body gravo-thermal experiments of systems with weak dissipation, i.e., of systems in quasi-equilibrium are presented. The results are compared with theoretical findings.

In order to cover a range of energy with the same experiment, we introduce a weak global dissipation scheme allowing us to describe a range of quasi-equilibrium states. Here weak means that $\tau_{\rm dis} \gg \tau_{\rm dyn}$, where $\tau_{\rm dis}$ is the dissipation time-scale. Then, the results can be compared with theoretical equilibrium states.

Follana & Laliena (2000) theoretically examined the thermodynamics of self-gravitating systems with softened potentials. They soften the Newtonian potential by keeping n terms of an expansion in spherical Bessel functions (hereafter such a regularized potential is called a Follana potential). This regularization allows the calculation of the thermodynamical quantities of a self-gravitating system. The form of their potential is similar to a Plummer potential with a corresponding softening length. Figure 2 shows a softened Follana potential with n=10, and a Plummer potential with $\epsilon = 0.05$.

Figure 2: The Follana (n=10) and the Plummer ( $\epsilon = 0.05$) potential, in units of Gm2/R, where m is the particle mass.

In their theoretical work Follana & Laliena found for a mild enough regularization (n<30) a phase transition below the critical energy $\varepsilon_{\rm c}\approx -0.335$ in a region with negative specific heat. The transition separates a high energy homogeneous phase from a low energy collapsed phase with a core-halo structure.

We want to reproduce these findings by applying a Plummer potential ($\xi =0$). Furthermore, the effect of a small-scale repulsive force ($\xi > 1/3$) is studied.

For these purposes, simulations with a weak global dissipation strength, $\alpha=0.025$, are carried out. The dissipation time is then $\tau_{\rm dis} = 80\;\tau_{\rm dyn} = 56.6\;\tau_{\rm ff}$, where the free fall time is $\tau_{\rm ff}=\tau_{\rm dyn}/\sqrt{2}$.

Before we discuss the results, let us briefly present some model properties. The initial state is a relaxed, unperturbed and confined N-body sphere with total energy $\varepsilon=1$.

Assigning a particle a volume of $4\pi\epsilon^3/3$, the volume filling factor is,

$$V_{\rm ff}=\frac{N \epsilon^3}{R^3},$$ (16)

where N is the particle number and R is the radius of the system. The volume filling factor is set as $V_{\rm ff}=1$, meaning that the force resolution is equal to the mass resolution. Then, the particle number is N=8000.

The applied weak dissipation strength maintains the system approximately thermalized and virialized (see Fig. 3). Here a system of N particles is called virialized when the moment of inertia I is not accelerated,

$$\frac{\ddot{I}}{2}=2T+\sum_{i=1}^N F_i r_i\approx0,$$ (17)

where T is the kinetic energy, r_i is the location of the i-th particle and F_i is the sum of all forces acting on the particle (gravitation, friction, confinement).

Figure 3: The potential U, kinetic T and total energy $\varepsilon$ as a function of time for a simulation with global dissipation and Plummer softening ($\xi =0$). The temporal evolution of the virial equation is indicated as well ( ${\rm Viriel} \equiv \ddot{I}/2$). The dashed vertical line depicts the time, when the system becomes fully self-gravitating. The dash-dotted and the dotted line mark the interval of negative specific heat.

The evolution of the inverse temperature $\beta=1/T$ as a function of the total energy $\varepsilon$ for two simulations with $\xi =0$ and $\xi =2/3$, respectively, as well as the semi-analytical curve calculated by Follana & Laliena are shown in the left panel of Fig. 4.

Figure 4: Comparison of theoretical predictions and simulated systems. The simulations are carried out with N=8000 and the softening length is, $\epsilon = 0.05$. Consequently, the volume filling factor is, $V_{\rm ff}=1$. Left: Inverse temperature $\beta$ versus energy $\varepsilon$ for models with softened potentials. The solid line indicates the theoretical result with the Follana potential, n=10. The other curves depict the evolution of the simulated systems. Crosses: repulsive potential ($\xi =2/3$). Dots: Plummer potential ($\xi =0$). The circle indicates the initial state of the simulations. The range of negative specific heat corresponds to the range where the slope of $\beta (\varepsilon )$ is positive. Middle: The dotted lines describe the evolution of the Lagrangian radii $\kappa$for the simulation with the Plummer softening potential ($\xi =0$). Each curve depicts the radius of a sphere containing a certain mass fraction. The different mass fractions are: $\Delta M/M=\{5\%, 10\%, 20\%, \ldots ,80\%, 90\%, 95\%\}$. The solid line depicts the theoretical $95\%$-Lagrangian radius for the Follana potential. The dashed vertical line depict the moment when the system becomes fully self-gravitating. The dash-dotted and the dotted line mark the interval of negative specific heat. Right: Idem for the simulation with the short range repulsive force.

For the simulation with the Plummer potential the interval of negative specific heat agrees with theoretical predictions. Also, in accordance with predictions, a phase transition takes place, separating a high energy homogeneous phase from a collapsed phase in a interval with negative specific heat. Yet, the simulated phase transitions occur at higher energies. To illustrate this the evolution of the Lagrangian radii $\kappa$ are shown in the middle and the right panel of Fig. 4.

In the high energy homogeneous phase the system is insensitive to the short-distance form of the potential. Thus all systems enter the interval of negative specific heat at the same energy. However, the collapsed phase is sensitive to the short-distance form. That is, the system with repulsive short-distance force re-enters the interval of positive specific heat at a higher energy than those without such forces (see left panel of Fig. 4).

Furthermore, the collapsed phase resulting from a simulation with a Plummer potential is hotter and denser than those resulting from a simulation with repulsive small scale forces. This can be seen in Fig. 4 as well.

For a smaller softening length, theory expects a phase transition at higher energies. Yet, already in simulations with a mild softening $(\epsilon=0.05)$ the collapsing transition occurs straight after the system has entered the interval of negative specific heat, which is shortly after the system has become self-gravitating. Thus, the collapsing transition cannot occur at significantly higher energies for a smaller softening length.

However, for systems with smaller softening length, the energy at which the system reenters the zone of positive specific heat after the collapse changes, i.e., it is shifted to smaller energies. This is shown in Fig. 5, where two simulations with $N=160\,000$ and $\epsilon = 0.05$ resp. $\epsilon =0.01$ are compared.

Figure 5: Same as Fig. 4 for two simulations with $N=160\,000$ and with a Plummer potential $\xi =0.0$. Dots: $\epsilon = 0.05$, $V_{\rm ff}=20.0$. Crosses: $\epsilon =0.01$, $V_{\rm ff}=0.16$.

Due to the higher particle number the dissipation time is in these simulations reduced to $\tau_{\rm dis}=4\tau_{\rm dyn}$. Thus, dynamical equilibrium is not as perfect as previously. That is, the acceleration of the moment of inertia I attains a maximum value of, $\ddot{I}/T=0.14$, where T is the kinetic energy. However, in general dynamical equilibrium is well approximated and we expect, due to the experience with simulations with varying particle number and dissipation time, that the effect of the temporal deviation from equilibrium does not affect qualitatively the results presented in Fig. 5.

The volume filling factor of the systems with $N=160\,000$ and $\epsilon = 0.05$ is $V_{\rm ff}=20$, meaning that the mass resolution is a factor 2.7 larger than the force resolution. Thus the system is less granular than in the previous simulations with N=8000 and $V_{\rm ff}=1$ (see Fig. 4). Yet, this does not affect the simulated quasi-equilibrium states and the deviation from theoretical predictions remains.

5.1 Discussion

Some predictions made by analytical models, such as the interval of negative specific heat, agree with findings resulting from N-body models. Yet, discrepancies also appear, such as the way the collapsing phase transition develops.

Figure 6: Evolution of the Lagrangian radii for a system with global dissipation scheme, dynamical friction and a local dissipation scheme, respectively. The Lagrangian radii depict, here and in the following figures, spheres containing the following mass fractions: $\Delta M/M=\{5\%, 10\%, 20\%, \ldots ,80\%, 90\%, 98\%\}$.

At first glance the deviating results are surprising. Yet the small scale physics is different in the two models, which may account for the discrepancies. Indeed, thermostatistics assumes a smooth density distribution and is thus not able to account for two-body relaxation effects in granular media. However, a certain degree of granularity is an inherent property of N-body systems. This may then lead to discrepancies, especially in the unstable interval of negative specific heat, where phase transitions may be sensitive to small scale physics.

While thermostatistics is too smooth to account for microscopic physics in granular self-gravitating media such as the interstellar gas, two-body relaxation is often too strong in N-body systems due to computational limitations. This is particularly so in high force resolution simulations where the force resolution is larger than the mass resolution. We refer to this in Sect. 6.1 where we discuss, among other things, the effect of the granularity on long-range correlations appearing in the interval of negative specific heat.

A further discrepancy between nature, analytical models and N-body models may be due to the entropy used in gravo-thermal statistics (Taruya & Sakagami 2001). Indeed, the entropy used in analytical models to find equilibrium states via the maximum entropy principle is the extensive Boltzmann-Gibbs entropy, that is in fact not applicable for non-extensive self-gravitating systems.

Generalized thermostatistics including non-extensivity are currently developed (Tsallis 1988; Sumiyoshi 2001; Latora et al. 2001; Leubner 2001). These formalisms suggest that non-extensivity changes not all but some of the classical thermodynamical results (Tsallis 1999; Boghosian 1999), which agrees with our findings.

Because currently it is a priori not known which thermostatistical properties change in non-extensive systems and consistent theoretical tools are not available, analytical results must be considered with caution.