On the clustering phase transition in self-gravitating N-body systems

The thermodynamic behaviour of self-gravitating $N$-body systems has been worked out by borrowing a standard method from Molecular Dynamics: the time averages of suitable quantities are numerically computed along the dynamical trajectories to yield thermodynamic observables. The link between dynamics and thermodynamics is made in the microcanonical ensemble of statistical mechanics. The dynamics of self-gravitating $N$-body systems has been computed using two different kinds of regularization of the newtonian interaction: the usual softening and a truncation of the Fourier expansion series of the two-body potential. $N$ particles of equal masses are constrained in a finite three dimensional volume. Through the computation of basic thermodynamic observables and of the equation of state in the $P - V$ plane, new evidence is given of the existence of a second order phase transition from a homogeneous phase to a clustered phase. This corresponds to a crossover from a polytrope of index $n=3$, i.e. $p=K V^{-4/3}$, to a perfect gas law $p=K V^{-1}$, as is shown by the isoenergetic curves on the $P - V$ plane. The dynamical-microcanonical averages are compared to their corresponding canonical ensemble averages, obtained through standard Monte Carlo computations. A major disagreement is found, because the canonical ensemble seems to have completely lost any information about the phase transition. The microcanonical ensemble appears as the only reliable statistical framework to tackle self-gravitating systems. Finally, our results -- obtained in a ``microscopic'' framework -- are compared with some existing theoretical predictions -- obtained in a ``macroscopic'' (thermodynamic) framework: qualitative and quantitative agreement is found, with an interesting exception.


INTRODUCTION
Many astrophysical problems, ranging from cosmology down to the formation of planets, demand a knowledge of the physical conditions for the existence of stable equilibria as well as for the appearance of instabilities in self-gravitating systems, either in the form of large gas clouds or of collections of concentrated objects (N -body systems). As a consequence, on this subject there is a long standing and important tradition of both theoretical and numerical investigations. In view of the aims of the present paper, there are three fundamental reference works. The first one is an old and interesting paper (Bonnor, 1956) where the equation of state (P − V ) was derived for a spherical mass of gas at uniform temperature in equilibrium under its own gravitation, and where the modified Boyle's perfect gas law was related with the onset of a gravitationally driven instability. The other two papers consider N gravitationally interacting particles in a spherical box which, under suitable physical conditions, can undergo a catastrophic state (the Gravothermal Catastrophe) characterized by the absence of a maximum of the entropy (Antonov 1962) or, in a suitable range of parameters, can exist in a clustered phase characterized by a negative heat capacity (Lynden-Bell & Wood 1968). Already in this latter paper, then confirmed by subsequent work (Aaronson & Hansen 1972), it was surmised that -in analogy with ordinary thermodynamics -self-gravitating systems undergo a phase transition. This hypothesis was also supported by the theoretical evidence of a condensation phenomenon in a system of self-gravitating hard spheres (Van Kampen 1964). Also numerical experiments on the N -body dynamics supported this analogy through the evidence of the core-halo formation (Aarseth 1963) or, in simulations of galaxy clustering (Aarseth 1979), through a phenomenology which is reminiscent of a phase transition occurring in atomic or molecular systems. More recent evidence of the existence of a phase transition in gravitational systems has been given for plane and spherical sheets (Reidl & Miller 1987;Miller & Youngkins 1998;Youngkins & Miller 2000) and for particles constrained on a ring (Sota, Iguchi & Morikawa 2000).
Though the "negative specific heat paradox" -arisen by Lynden-Bell & Wood in their explanation of Antonov's gravothermal catastrophe -seems to faint the analogy with laboratory systems, Hertel and Thirring (Hertel & Thirring 1971) showed that a negative specific heat -which is strictly forbidden in the canonical ensemble of statistical mechanics -can be legally found in the microcanonical ensemble as a consequence of the breakdown of ensemble equivalence; this ensemble inequivalence has been later confirmed for a simplified model of gravitationally interacting objects (Lynden-Bell & Lynden-Bell 1977) ⋆ .
All the results of the three above quoted papers were worked out in a macroscopic thermodynamic framework. Thus we can wonder if and how such a macroscopic phenomenology can be retrieved in the framework of a microscopic dynamical description. In other words, in the spirit of statistical mechanics, we can try to link the dynamics of the elementary constituents of a self-gravitating system with its large scale thermodynamics. (Throughout this paper we use "macroscopic" and "microscopic" in a relative sense, thus the microscopic level can be that of individual atoms or molecules of a gas cloud, as well as that of stars in a galaxy). This is actually one of the aims of the present paper, where we address the problem of the statistical mechanical treatment of self-gravitating systems in order to understand whether the known thermodynamics can be derived from purely gravitational dynamics. The other aim of the present paper, closely related to the previous one, is to tackle the clustering instability, i.e. the breakdown of the homogeneous phase into a core-halo phase, in the attempt to describe it in analogy with the other -better understood -laboratory phase transitions.
The paper is organized as follows. In Section 2 we briefly discuss why a-priori a statistical mechanical treatment of gravitational N -body systems runs into difficulties and, on the other hand, why these systems can be considered bona fide ergodic after suitable regularization; moreover, we present an unconventional regularization strategy of the exact Hamiltonian (transformed into a "Fourier-truncated" Hamiltonian), leading to models that do not afford any practical advantage in the numerical treatment of the Nbody dynamics but that are of great prospective interest for its analytic treatment in a statistical mechanical context and thus to apply concepts and methods that already proved powerful in the study of Hamiltonian dynamics. In Section 3, after a quick glance at the integration method and parameters of the equations of motion and at the Mon-teCarlo simulations of the canonical ensemble averages, we ⋆ The inequivalence of statistical ensembles is not limited to gravitational systems, but it occurs whenever the range of the potential is comparable with the size of a system (Cipriani & Pettini 2000) and has recently attracted a lot of interest (Gross 1997).
show how thermodynamics can be worked out through dynamics with the aid of microcanonical statistical ensemble and by borrowing from Molecular Dynamics useful formulae already in the literature. We report the caloric curves (temperature vs energy), the specific heat and an order parameter vs energy (suitably scaled with N ), and a comparison among the outcomes of the simulations obtained with the Fourier-truncated Hamiltonian and those obtained with the standardly softened Hamiltonian. A somewhat familiar phenomenology in the numerical study of phase transitions is found mainly through the order parameter. A strikingly strong ensemble inequivalence is found: the canonical Mon-teCarlo averages do not keep even the slightest trace of the change of the specific heat from positive to negative and of the phase transition. Having confirmed the a-priori expected ensemble inequivalence, we discuss the reasons to consider the microcanonical ensemble as the good representative statistical ensemble for self-gravitating N -body systems and we present the dynamically worked out P − V equation of state. In Section 4 we compare our results with their theoretical counterparts obtained within a thermodynamic macroscopic approach. Qualitative and quantitative agreement is found and everything fits into a coherent scenario, though with a remarkable difference between the clustering transition and the gravothermal catastrophe. Finally, in Section 5 some conclusions are drawn. In Appendix A some basic definitions and concepts of the Gibbsian ensemble formulation of statistical mechanics are briefly recalled. In Appendix B a technical detail about the mean-field approximation is sketched.

SELF-GRAVITATING N-BODY SYSTEMS
In principle star clusters, galaxies, clusters of galaxies seem to naturally call for a statistical description of their dynamical behaviour. However, as the existence of negative specific heats reveals, there are some difficulties due to the very special nature of gravitational interaction. The gravitational interaction is always attractive, unscreened and of infinite range, therefore it is not stable, i.e. the potential energy U (r1, . . . , rN ) of N gravitationally interacting masses does not fulfil the condition U (r1, . . . , rN ) ≥ −N B, with B a positive constant, whence the so-called lack of saturation: in the N → ∞ limit the binding energy per particle and the free energy per particle diverge. This is also referred to as a breakdown of extensivity of these fundamental physical quantities, at variance with what is familiar for ordinary laboratory systems. Moreover, from a rigorous mathematical viewpoint, there is no equilibrium ground state: as N increases, at N/V = const, the energy per particle increases with N and cannot be defined in the thermodynamic limit † . The absence of an equilibrium state means that a gravitational system does not behave thermodynamically because standard thermodynamics does not apply to evolving systems. However, a thermodynamic description is still possible for self-gravitating N -body systems provided that they are † Though, from a physical point of view, the thermodynamic limit does not exist and has no meaning for self-gravitating systems.
not in a strongly unstable phase (Saslaw 1985), and in the case of slowly evolving system we can resort to statistical mechanics.
Moreover, statistical mechanics gives a correct thermodynamics if the thermodynamic potentials are extensive quantities. In a self-gravitating system only entropy in the microcanonical ensemble -at least on finite time scales and if the system is slowly evolving -is an extensive thermodynamic potential ‡ ; therefore, the results obtained in other ensembles (grancanonical and canonical) are a-priori expected to be in disagreement with the results obtained in the microcanonical ensemble.
In what follows, systems of N gravitationally interacting point masses will be considered that are described by the Hamiltonian function where ri ≡ (xi, yi, zi). For the sake of simplicity we shall put G = 1 and mi = 1, i = 1, . . . , N . We remark that this choice is convenient to keep the statistical and thermodynamic properties unaltered (e.g. the non-extensivity), at variance with what is implied by other possible choices (Heggie & Mathieu 1986).

Dynamics and statistics
Any physical phenomenon occurring in the system described by Hamiltonian (1) must have its origin in the dynamics. With the exception of those systems that can be treated with the methods of celestial mechanics (mainly perturbation theory), the dynamics of generic N -body systems can be worked out only through numerical simulations. The numeric approach, which is continually improved (Aarseth 1999;Meylan & Heggie 1997, and references quoted in these papers), provides the "experimental" counterpart of the theoretical investigation of these systems. In order to describe the physics of N -body systems through a few relevant macroscopic observables, dynamics must fulfil the requirements of ergodicity and mixing to justify the application of statistical mechanics. It is in general impossible for Hamiltonian systems of physical interest to ascertain whether ergodicity and mixing are rigorously verified. However, a generic non-integrable and chaotic Hamiltonian system with a large number N of degrees of freedom can be considered ergodic and mixing, at least in a physical sense. In fact, after a famous theorem due to Poincaré and Fermi (Poincaré 1892;Fermi 1923a,b), generic systems with three or more degrees of freedom are not integrable, i.e. there are no nontrivial invariants of motion besides total energy. Only global invariants (like total momentum and angular momentum), due to global symmetries of the Hamiltonian (1) (as space translations and rotations), can exist. Thus, once the initial condition of a gravitational N -body system is assigned, the representative point of the system ‡ From the fundamental relation 1/T = ∂S/∂E, being T proportional to the kinetic energy per particle and making use of the virial theorem, N T and E must have the same N dependence hence S is extensive. moves on a 6N − 10 dimensional hypersurface of its 6Ndimensional phase space. The lack of nontrivial integrals of motion (i.e. not related to global symmetries) entails that all the 6N − 10 dimensional hypersurface of phase space is accessible to a trajectory issuing from any initial condition. A possible source of non-ergodicity would be the coexistence of regular and chaotic motions: observed for the first time in the Hénon-Heiles model (Hénon & Heiles 1964), it is in general implied by the Kolmogorov-Arnold-Moser (KAM) theorem (Thirring 1978), which, however, has no practical relevance at large N § . Thus, self-gravitating N -body systems, after some suitable regularization to make finite the phase space volume, can be considered ergodic, so that time averages of physically relevant observables can be replaced with suitable statistical averages computed on a given 6N − 10 dimensional hypersurface of phase space.
The dynamical instability of the gravitational N -body systems (Miller 1964) implied one among the first numerical evidences of the existence of deterministic chaos in Hamiltonian dynamics. Since then, further numerical evidence of the existence of chaos in the self-gravitating systems has been provided by several authors (Kandrup 1990;Quinlan & Tremaine 1992;Goodman, Heggie & Hut 1993;Kandrup, Mahon & Smith 1994, and references therein;Cipriani & Pucacco 1994;Cerruti-Sola & Pettini 1995;Cipriani & Di Bari 1998;Miller 1999). On the other hand, chaotic dynamics in a many dimensional phase space implies a bona fide phase mixing (Casetti et al. 1999), which means that time averages of physical observables converge to their statistical counterparts in a finite time (whereas ergodicity implies an infinite time convergence).
Therefore, the use of microcanonical statistical mechanics is naturally motivated by the reasonably good ergodicity and mixing properties of the dynamics of regularized self-gravitating N -body systems (Cipriani & Pettini 2000). Taking into account the conservation of energy, of center of mass position and momentum, and of angular momentum, the microcanonical volume in phase space reads In principle, by means of ωN (E), we can compute statistical averages of any physical observable defined through a function A(p, r), and we can also compute the thermodynamics of a self-gravitating N -body system, the link being made by the entropy defined as S = kB log ωN (E); kB is the Boltzmann constant (see Appendix A).

Regularized N-body Hamiltonians
Let us now consider a system constrained in a finite volume. At variance with the customary choice of a spherical box § KAM theorem requires extremely tiny deviations from integrability to imply the existence of sizeable regions of regular motions in phase space, moreover these deviations from integrability cannot exceed a threshold value which drops to zero exponentially fast with N (Thirring 1978; Casetti et al. 1999) . (Bonnor 1956;Antonov 1962;Lynden-Bell & Wood 1968), let us consider a cubic box of side length equal to L. The reason for such a choice is to explicitly break the rotational invariance of the Hamiltonian (1) and, in so doing, the microcanonical volume ωN (E) simplifies to ¶ and thus we can borrow from the existing literature the analytic expressions, derived using ωN (E) of Eq.
(3), of some basic thermodynamic observables that are then used in our numerical computations. Working out anew the same analytic expressions using ωN (E) of Eq.
(2) is a non trivial task beyond the aims of our present investigation. However, this is not a severe restriction if we refer to almost non-rotating systems whose angular momentum, even if conserved, is negligible. The cubic box, allowing fluctuations of the total angular momentum, is thus equivalent to considering a whole ensemble of almost vanishing angular momentum systems. Even though the assumption of a confining box could seem unphysical, it is a simple way of idealizing different physical aspects which depend upon the chosen boundary conditions. In fact, the assumption of periodic boundary conditions is as if we took a fragment out of the bulk of a large system where -in the average -the number of particles remains constant and small local density and energy fluctuations of the subsystem take place. In this case when a mass exits the box in a given direction, another mass with the same energy enters the volume from the opposite side to keep constant the energy and the number of particles. The assumption of reflecting boundary conditions amounts to mimic the presence of a halo of diffuse matter whose gravitational potential field would act to confine the system. Both assumptions about the geometric constraints of the system, besides the explicit breaking of rotational symmetry, also guarantee the finiteness of the configuration space volume over which the integration in Eq. (3) is performed.
In order to guarantee also the boundedness in momentum space, so that the whole integral in Eq.(3) extends over a finite region of phase space, the two-body interaction potential must be regularized. We adopted two different kinds of regularization. The first one is the standard softening, adopted in numerical simulations, with the replacement (e.g. see Binney & Tremaine 1987) where η is a small softening parameter that bounds below the interaction potential and in so doing prevents the occurrence of arbitrarily large values of the momenta. This regularization is local in space. The second regularization is nonlocal in space. It makes use of the Fourier representation of the Green function G(r− r ′ ) for the Poisson equation (5) ¶ In fact, both periodic conditions and reflecting boundary conditions destroy the conservation of P and R.
In order to regularize the two-body potential in Hamiltonian (1), one can truncate the Fourier expansion in Hamiltonian (7) by summing l, m, n from 1 up to a finite number Nw.
The two regularizations are a-priori inequivalent: the former pertains events (close encounters) localised in real space, whereas the latter makes use of collective coordinates (the Fourier modes) which are not localised in real space. By truncating the Fourier expansion in Hamiltonian (7) we neglect all the dynamical details occurring at length scales smaller than the smallest spatial wavelength in the expansion. Loosely speaking, this is reminiscent of standard methods in statistical mechanics, mainly in the context of the renormalisation group theory, where relevant and irrelevant degrees of freedom at a given length scale are Fourier modes with wavelengths above and below some cutoff respectively. In order to ascertain to what extent a truncated model still retains some physically relevant feature of the exact N -body system (1), it is necessary to compare the outcomes of different truncations, i.e. different Nw, and to make a comparison between the results obtained by simulating the dynamics associated with Hamiltonian (1) and with a truncated version of the Hamiltonian (7).
There are different reasons for making numerically heavier the already heavy N -body problem. First of all, the spatial coarse-graining considerably lessens the dramatic effect that close encounters have on a reliable numerical computation of Lyapunov exponents (though they are not reported here), as well as on other observables of dynamical relevance; the frequency and the effect of close encounters in numerical simulations of N -body systems are in this respect too much enhanced in comparison with the physical reality of practically collisionless systems as galaxies, cluster of galaxies and, to a lesser extent, star clusters. Second, the approach to equilibrium, as well as many other nontrivial dynamical properties, are by far better revealed by the use of collective coordinates, as a long experience in a different context has widely witnessed (Casetti et al. 1999). Third, truncated Hamiltonians out of (7) have the interesting property of naturally allowing a mean-field-like decoupling of the degrees of freedom which is very interesting in view of analytical computations of both statistical mechanical properties and of chaotic properties of the dynamics in the framework of a geometric theory of Hamiltonian chaos (Casetti et al. 2000); moreover, such a mean-field-like representation of the Fourier truncated Hamiltonians leads to the definition of an order parameter that is useful to detect the occurrence of a phase transition. The regularization helps to smooth the local fluctuations on small spatial scales which do not significantly affect the macroscopic behaviour of the system, but cause a noisy variation of the averages on the time scales of interest. Finally, a wealth of models could be generated, even by retaining very few modes, which could reveal different aspects of the astonishingly rich dynamics of the exact N -body problem.
A generic truncated model Hamiltonian reads where we have chosen G = 1 and mi = 1, i = 1, . . . , N . From the associated Lagrangian function -being pxi =ẋi, pyi =ẏi, pzi =żi -the following equations of motion are derived i = 1, . . . , N , and where we have introduced S (i) l,m,n = S l,m,n − sin(k l xi) sin(kmyi) sin(knzi), with S l,m,n = N j=1 sin(k l xj) sin(kmyj) sin(knzj ), to put in evidence that the same quantities S l,m,n enter all the equations of motion, what obviously simplifies the numerical computations.

NUMERICAL RESULTS
The phase space trajectories of an Hamiltonian system are constrained on a constant energy surface in phase space; therefore, time averages computed along the numerical solutions of the equations of motion, of both the Fouriertruncated system and of the softened version of the Hamiltonian (1), are generically expected to converge to microcanonical ensemble averages (see Appendix A). Thus, in order to work out the thermodynamics of self-gravitating N -body systems, we borrow from Molecular Dynamics the Both regularizations do not qualitatively alter the chaoticity of the dynamics.
At variance, the numerical estimate of canonical ensemble averages requires to construct suitable random markovian walks in the full phase space. Along such random trajectories (no longer constrained on any energy surface), the averages of physical observables converge to canonical ensemble averages provided that the recipe to generate the random walk is that of the so-called Metropolis importance sampling of the canonical ensemble weight.

Numerical integration
The numerical integration of the equations of motion (9) and of the equations of motion derived from the Hamiltonian (1) with the replacement (4) has been performed by means of a symplectic algorithm (McLahlan & Atela 1992). Some runs to check the reliability of the results have been performed using also a bilateral scheme (Casetti 1995). Symplectic integrators update the coordinates and momenta of an Hamiltonian system through a canonical transformation of variables; for this reason, symplectic algorithms ensure a faithful representation of an Hamiltonian flow because, in addition to the conservation of phase space volumes and of the energy, they guarantee the conservation of all the Poincaré integral invariants of a system. Actually, there are interpolation theorems (Moser 1968;Benettin & Giorgilli 1994) stating that the numerical flows obtained through symplectic algorithms can be made as close as we may wish to the true flow of a given Hamiltonian. Though locally other integration schemes can give more precise results in presence of close encounters by using different time steps for individual particles (Aarseth 1985), the non-symplectic nature of these algorithms might a-priori somewhat alter the frequency with which different regions of phase space have to be visited by an ergodic dynamics, at least for very long runs. On the other hand, we are just interested in long runs, so that time averages of the chosen observables display a good stabilization, and, in order to safely replace microcanonical ensemble averages with time averages, dynamics has to properly sample the phase space.
The dynamics of the direct N -body system has been numerically computed using integration time steps ∆t ranging in the interval (5 · 10 −5 − 10 −3 ): the relative variations ∆E/E of the total energy were in the range (10 −10 , 10 −4 ) on integration runs of (10 6 , 10 7 ) time steps. The softening parameter η has been set to 0.01 and scaled as η = ηd, where For what concerns the Fourier-truncated model, computations have been done with Nw = 5, 7, 10, i.e. a total number of modes equal to 125, 343 and 1000 respectively. The initial conditions on the particle coordinates (xi, yi, zi) have been chosen at random with a uniform distribution in the interval (0, L). The initial momenta have been chosen at random with a gaussian distibution of zero mean and a suitable variance to approximately match the desired initial value of the total energy. An opportune velocity rescaling allowed to precisely fix the initial value of the total energy. With integration time steps ∆t ranging in the interval (0.01, 0.001), we got relative variations ∆E/E of the total energy in the range (10 −8 , 10 −6 ) with zero mean, i.e. without any drift, and on long integration runs of 10 6 time steps. For E < 0, the system has been then let evolve for about 3tD, where the dynamical time is defined as tD ∝ (N/|E/N | 3/2 ) (this makes the system virialize); for E > 0, we let the system evolve until it attains a stationary state of dynamical equilibrium between kinetic and potential energies.
By varying the side L of the box, the volume V has been varied. With both periodic and reflecting boundary conditions, the initial velocity of the center of mass has been set equal to zero. The initial total angular momentum can be made very small, but then it fluctuates because of the explicit breaking of rotational symmetry. The larger N , the smaller such fluctuations can be.

The N-scaling of the results
The number N of interacting bodies has been varied in the range 25 − 500, with some checks up to N = 2000. In the numerical simulations of extensive Hamiltonian systems, i.e. with short-range interactions so that energy and other physical observables are additive, the comparison among the results obtained by varying N is naturally made through densities: the values of the observables divided by N .
In what follows, we shall mainly vary the energy E by keeping the ratio ̺ = N/V constant. From a dimensional point of view, the potential energy . As we shall see in the following, this actually suggests the correct way of scaling the results obtained for different values of N .

MonteCarlo computations
For standard Hamiltonians H = i 1 2 p 2 i + U (q), the weight e −βH splits into a multidimensional gaussian, originated by the kinetic energy term, and into a configurational part e −βU (q) ; β is the inverse of the temperature. Only this latter part is dealt with by the MonteCarlo algorithm. The system under consideration has to be at equilibrium so that the transition probabilities W (1 → 2) and W (2 → 1) between any two configurations "1" and "2" satisfy the condition W (1 → 2) = W (2 → 1), the so-called detailed balance. A MonteCarlo move consists of a random update of the coordinates. If the system lowers its energy with the proposed configuration update, then the new configuration is accepted, otherwise the transition probability where N is the normalization factor, is compared with a random number ζ uniformly distributed in the interval [0, 1]; if W > ζ then the new configuration is accepted, otherwise the old configuration is counted once more. This is the essential of the Metropolis importance sampling algorithm (Binder 1979). The average acceptance rate is usually kept not far from 50 per cent by adjusting the mean variation of the coordinates at each update proposal. The MonteCarlo estimate of the canonical average of an observable A is then simply given by the sum being carried over the NMC accepted configurations {qi}.

Dynamical analysis of thermodynamical observables
In deriving the thermodynamics of self-gravitating N -body systems, one of our aims is to address the problem of the existence and of the characterization of the clustering phase transition. It is not out of place to remind here that the basic thermodynamical phenomenology of a phase transition is characterized by the sudden qualitative change of a macroscopic system when its temperature varies across a critical value Tc. This sudden change is qualitatively due to collective microscopic behaviours and is quantitatively reflected by the singular ⋆⋆ temperature dependence of the most relevant thermodynamical observables. Changes of state, like melting, condensation and so on are examples of first order phase transitions, the spontaneous magnetization of a paramagnetic material when temperature is lowered below the so-called Curie temperature is an example of a second order phase transition † † . The usual framework of both theoretical and numerical investigations of phase transitions is mainly that of Gibbs' canonical ensemble (see Appendix A).
Recently, the question of whether microscopic Hamiltonian dynamics displays some relevant change at a phase transition has been addressed by several works (Antoni & Ruffo 1995;Caiani et al. 1997;Dellago & Posch 1997, and earlier works therein cited; Gross 1997;Casetti et al. 2000;Cerruti-Sola, Clementi & Pettini 2000). The dynamical approach has its natural statistical mechanical counterpart in the microcanonical ensemble (see Appendix A). After a relaxation time which depends on the initial condition, the time averages of observables computed along numerical phase space trajectories provide good estimates of their microcanonical counterparts. In the following Sections we adopt this dynamical approach to investigate if also in self-gravitating N -body systems a phase transition is present and if the dynamics shows a corresponding qualitative change.

Caloric curve
The first goal of our numerical computations was to obtain the caloric curves of different self-gravitating systems resulting from different choices of the parameters. The caloric curve T (E) (temperature vs energy) represents a basic link between the microcanonical statistical ensemble and thermodynamics, thus between dynamics and thermodynamics. From the entropy definition S(E, V, N ) = kB log ωN (E) in the microcanonical ensemble, ωN being given by Eq. (3), temperature is derived as which, for systems described by standard Hamiltonians H = i 1 2 p 2 i + U (q), yields the expression ⋆⋆ True singularities are possible only in the N → ∞ limit. Here "singular" means that the larger N the sharper is the almost singular pattern of an observable as a function of temperature. † † According to the Ehrenfest's definitions, the first or the second derivative of the Helmoltz free energy with respect to temperature is discontinuous at a first order transition point or at a second order transition point, respectively.
where K −1 ω is the microcanonical average of the inverse of the kinetic energy K. The replacement of this average by means of the time average K −1 t provides the dynamical estimate of the temperature.
With the equivalent definition of microcanonical entropy as S(E, V, N ) = kB log ΩN (E) (see Appendix A), Eq.(10) yields the more familiar expression relating temperature with the mean kinetic energy per degree of freedom. In our numerical computations, the two temperatures have been always found almost coincident, well within the theoretically expected deviation of O(1/N ). Anyway, in the results reported here, T is computed according to Eq.(11).
The Fourier-truncated model, due to the finite number of modes Nw considered, underestimates the potential energy with respect to the direct model. However, since an additive constant in the Hamiltonian (8) does not affect the dynamics, we can shift a-posteriori the energy scale by computing the average energy difference of a set of random configurations whose potential energy is computed according to both Eqs. (1) and (8). Thus, in the following, in order to compare by superposition the results of the two models, the values of the energy densities obtained from the Fourier truncated model have been shifted towards the energy densities of the direct simulation.
In Fig. 1 the caloric curve is reported for the Fouriertruncated model (8) with Nw = 7, N = 50 interacting objects and L = 50, so that ̺ = 4 · 10 −4 . The results obtained with periodic and reflecting boundary conditions are synoptically displayed together with the outcomes of MonteCarlo simulations in the case of reflecting boundary conditions. In the case of MonteCarlo simulations, temperature is the input parameter and the average total energy is the outcome; anyway, since there is a one-to-one correspondence between temperature and the average total energy, the results can be reported without ambiguity also on the T vs E plane. The dynamical (microcanonical) caloric curves display the same qualitative features and are even quantitatively not very different one from another. The low energy branch has a negative slope, meaning that the specific heat is negative. At some value of the energy per degree of freedom, the slope of T (E) becomes positive and the curve bends towards an asymptotic straight line of slope 2/3, proper to a gas of independent particles. Fig. 2 shows the effect of different truncations in Eq. (8): qualitatively, the results are satisfactorily stable, quantitatively, the larger Nw, the closer is dT /dǫ (with ǫ = E/N ) to the newtonian slope −2/3.
Notice that Figs. 3 and 4 display also a slope −2/3 of the negative-cV branches, in agreement with the prediction of the virial theorem, in the regime in which the effect of the box is negligible (large negative energy values).  data collapsing ‡ ‡ is obtained by rescaling the energy with N 5/3 and temperature with N 2/3 (notice that temperature is already an energy divided by N ). The simple argument given in Section 3.1.1 appears correct.
There is a very good agreement between the results obtained with the Fourier-truncated model and those obtained with the direct N -body simulations. This is a very interesting result because of the relatively small number of Fourier modes considered, and in view of adopting the Fourier-truncated model for the analytic computation of statistical mechanical averages.
It is remarkable that the change of sign of the slope of the caloric curve implies the existence of a phase transition. In fact, T (E) can only change the sign of its slope either through a "V"-shaped pattern or through an "U"-shaped pattern: in the former case, the derivative ∂T /∂E is discontinuous at the cuspy point so that the specific heat makes a discontinuous jump (first-order phase transition), whereas in the latter case the same derivative vanishes at the minimum so that the specific heat diverges (second order phase transition).
The comparison among the dynamical (microcanonical) and the MonteCarlo (canonical) caloric curves shows a remarkable fact: the canonical caloric curve seems indistinguishable from that of a perfect gas, no reminiscence in the canonical results seems to be left of the feature shown by the microcanonical curves. As already mentioned throughout this paper, a negative specific heat is strictly forbidden in ‡ ‡ This is standard jargon in the numerical study of phase transitions, and it means that when different sets of results are compared after a suitable rescaling with the parameter that varies from a set to another, then all the results crowd on the same curve.   the canonical ensemble; however, in the case of a two dimensional model, it has been observed (Lynden-Bell & Lynden-Bell 1977) that the appearance of negative specific heat in the microcanonical ensemble corresponds to a signal of a phase transition in the canonical ensemble. Nothing similar is found here. This is a remarkable result: though a-priori the ensemble inequivalence was expected, such a radical loss of any signal of the transition in the canonical ensemble was unexpected.

Specific heat
From the above given definition of temperature the microcanonical specific heat is computed according to the formula Hence the knowledge of the caloric curve implies the knowledge also of the specific heat, but in practice the points on the T − E plane are inevitably affected by a "noise" which can make too rough the computation of CV as the numerical derivative of the caloric curve. The interest of an independent numerical derivation of the specific heat is therefore twofold: on one side, it constitutes a necessary and useful check of the previous results, on the other side it should give more information about the order of the phase transition. The already mentioned Laplace transform technique (Pearson et al. 1985) yields the following expression which has been used in our numerical simulations; K is the total kinetic energy. The specific heat vs energy qualitatively displays a pattern in full agreement with what can be guessed from the caloric curves. At small energy values a branch of negative specific heat is found. At an energy value slightly below the minimum of the caloric curve, the specific heat shows two marked spikes, stressing the occurrence of a sudden jump between negative and positive values. The high energy asymptotic value is 3/2, in quantitative agreement with the high energy slope of 2/3 of the caloric curve. Again the data collapsing is obtained by scaling the total energy with N 5/3 . In Fig. 5 the specific heat obtained with periodic and reflecting boundary conditions respectively is reported vs energy density for the Fourier-truncated model. Also through this observable it is confirmed that there is a weak sensitivity of the outcomes upon the boundary conditions, at least as far as thermodynamical observables are concerned.
For the Fourier-truncated model, Fig. 6 and Fig. 7 show -at different values of the density -how the pattern of cV changes when N is increased. The corresponding results obtained for the direct N -body simulations are reported in Fig. 8 and Fig. 9. A common feature of these results is the ambiguous behaviour of the height of the peaks of cV , close to the transition point, when N is increased, at variance with other systems with short-range interactions (Caiani et al. 1998;Cerruti-Sola et al. 2000). The doubtful absence of a systematic increase with N of the peak of cV could be due  to three main reasons: i) a very fine tuning of the energy value could be necessary and the right energy values could have been missed; ii) the transition could be too mild, for example with only a logarithmic divergence of cV ; iii) the pattern of T (E) could be "V"-shaped, thus ∂ 2 S/∂E 2 should make a finite jump at the transition point, with the physical consequence that the system should undergo a first order phase transition.  The high energy branches of cV , above the transition energy, obtained with the Fourier-truncated model and with the direct N -body simulation are in very good agreement, whereas the quantitative agreement is less good in the branch of negative values. Notice that in the Fouriertruncated model there is a lower bound to the energy values which depends on the number of modes retained in the trun- cated Fourier expansion, the larger the number of modes, the lower is the minimum attainable value of negative energies § § . Figures 5 through 8 show that the results of the computations of cV according to Eq.(13) are in very good agreement with the corresponding patterns of T (E). In both the direct newtonian simulations and in the Fourier-truncated model, a kink in the negative slope region of T (E) is found and is confirmed by the cV (E/N 5/3 ) pattern close to the transition point. We attribute the appearance of this feature to the competing effects of the fixed size of the box and the increasing (with energy) gravitational radius of the system.

Order parameter
Both in the theoretical and numerical studies of phase transitions, the choice of a good order parameter for a given system is a main point. An order parameter is typically a collective variable that bifurcates at the transition temperature Tc. A classic example is the spontaneous magnetization of a ferromagnetic material: it is zero above Tc, then it suddenly bifurcates away from zero at Tc and increases as temperature is lowered below Tc. True singularities exist only in the N → ∞ limit; at finite N they are smoothed but it is still possible to detect the existence of a phase transition by changing N and observing if the smoothed signals tend to sharpen or not.
In many cases the definition of an order parameter is § § Also for a softened newtonian system there is a lower bound on the total energy. With the values of the softening chosen here, the lower bound is much smaller than that of the Fourier-truncated model.
not trivial. A-priori this is true also for self-gravitating Nbody systems. However, let us notice an interesting property of the Hamiltonian (8) and of the equations of motion (9): at large N the coefficients S (i) l,m,n can be approximated by the S l,m,n with an error O(1/N ), and the replacement of the S l,m,n by the averages S l,m,n , computed through a consistency equation (see Appendix B), decouples all the degrees of freedom. The Hamiltonian (8) is approximated by Nw l,m,n=1 S l,m,n (l 2 + m 2 + n 2 ) sin(k l xi) sin(kmyi) sin(knzi) .
which is now in the form of a sum over independent degrees of freedom, because the S l,m,n are collective variables or, equivalently, order parameters, that take the same values for all the coordinates xi, yi, zi. The replacement of the coefficients S l,m,n with the averages S l,m,n in Hamiltonian (14) is a typical mean-field approximation in statistical mechanics (the S l,m,n are like Weiss molecular field of the early times of statistical mechanics of magnetic materials). The recasting of Hamiltonian (1) into the approximate form (14), which can be as precise as we may want if N , Nw and L are sufficiently large, is of great prospective interest for theoretical -analytic or semi-analytic -computations of statistical mechanical kind (see Appendix B). Out of the huge family of order parameters S l,m,n , it is natural to define the following two global order parameters in analogy with other more standard contexts (Caiani et al. 1998 l + m + n (l 2 + m 2 + n 2 ) |S l,m,n | .
The coefficients (l +m+n)/(l 2 +m 2 +n 2 ) are the sum of the corresponding coefficients in the equations of motion (9) for the three variables x, y, z. They decrease as the norm of the wavevector k = (l, m, n) increases, compensating the growing number of the S l,m,n terms as the norm of k = (l, m, n) increases. We easily realize that M1 cannot be very sensitive to a change of the particle distribution in the box because the dishomogeneities measured at different scales tend to cancel each other. Thus, to get rid of this problem, M2 is defined through the sum of the absolute values of the S l,m,n . Actually, the numerical computations, performed at different energies and N , confirmed that M1 is always very close to its minimum value (which depends on Nw). At variance, the order parameter M2 is very effective to detect the clustering transition and to make strict the analogy between this transition and a thermodynamic phase transition. In order to measure the deviation from the homogeneous spatial distribution of particles, it is convenient to slightly modify the order parameter as follows: we denote with M0 the value of M2 that corresponds to a uniform density of particles in the box, then we modify the order parameter M2 to µ2 = ν(M2 − M0), where ν is a normalization constant. Thus µ2 varies in the interval [0, 1], being zero for a uniform distribution of particles and one for a maximally clustered configuration. Such a normalization makes easy the comparison of the results obtained at different N . In the case of the Fourier truncated model, the maximally clustered configuration is obtained when all the particles are concentrated in a cell whose side is equal to the shortest wavelength. However, for simplicity, we empirically normalized µ2 by using the largest value of M2 measured at the lowest accessible energy for the given truncation of the Fourier expansion.
The values M0 of M2, corresponding to a uniform density of particles in the box, have been numerically computed by averaging M2 on a few hundreds of uniformly distributed random configurations of 50, 100, 200 and 400 particles. For N = 50, 100, 200, 400 and Nw = 7 we found M0 = 3.49, 2.60, 1.97, 1.54 respectively.
The averages |S l,m,n | are numerically computed as time averages to give M2, and the order parameter µ2 is worked out as a function of the energy and of the number of particles by using the above reported values of M0.
The same order parameter µ2 has been computed for the direct N -body simulations. Also in this case we considered Nw = 7, though in principle this choice is arbitrary because Nw is independent of the dynamics. Since a phase transition is a collective effect, driven by the long wavelength modes, there is no need for a large value of Nw.
The final results are reported versus E/N 5/3 in Fig. 11 for the Fourier-truncated model, and in Fig. 12 for the direct N -body simulations.
The data reported in these figures display typical patterns that are found in presence of a phase transition in other contexts. In fact, the order parameter µ2 displays a classic bifurcation pattern at the transition point, smoothed by the finite and not very large number of particles but with a clear tendency to become sharper at larger N as is shown by Fig.  11. The simulation of the direct N -body system allows to work out the bifurcation pattern of the order parameter in a larger energy interval, because there is no lower bound to the energy ¶ ¶ . The inset of Fig. 12 displays the pattern of µ2 on a large energy interval. The vertical dotted lines in both figures indicate the energies of the peaks of specific heat. In general, the inflection point of the order parameter and the location of the peak of the specific heat do not coincide but tend to get closer at larger N , thus making more precise the phase transition point. The bifurcation pattern of µ2 gives a reliable indication about the order of the phase transition: it strongly indicates a second order transition, thus eliminating the ambiguity of the information given by cV . In fact, in presence of a first order transition, the order parameter µ2 should make a finite jump at the transition point which is clearly not our case. We have formulated in Section 3.3.2 two possible explanations (i) and ii)) -alternative to the first order phase transition -of the behaviour of cV which can reconcile the outcomes on cV with those on µ2.
Notice that M2 is also a measure of the average modulus | Fi| of the force acting on each mass and -independently of N -this is never negligible for a bounded and clustered system, whereas in the homogeneous phase | Fi| depends on the relative density fluctuations and thus it is expected to decrease with increasing N . The physical consequence is that, in the homogeneous phase, the larger N the smaller is the relative average weight of the potential energy with ¶ ¶ As already noticed, the smoothing parameter entails a lower bound to the energy also in this case, but it occurs at a very large negative value.
respect to the kinetic energy so that the system behaves like a collection of almost independent particles, i.e. not far from a perfect gas.

Equation of state
The equation of state of a system contains the basic information of how pressure changes with volume at constant temperature. This information is contained in a family of isothermal curves in the P − V plane. If a phase transition is there, then the isotherms at T < Tc display a substantial difference from those at T > Tc, as an example, the isotherms of the van der Waals equation, describing the liquid-vapour transition, display a kink at T < Tc (Huang 1963).
Since dynamics has its natural statistical counterpart in the microcanonical ensemble, where energy rather than temperature is the independent variable, the isotherms are replaced by isoenergetic curves in the P − V plane.
From thermodynamics we have ∂S ∂V E = P T and using the above given definition of temperature one gets where ΩN is defined in the Appendix A. Hence whence, with the aid of the already mentioned Laplace transform technique, one obtains (Pearson et al. 1985) which is the microcanonical equation of state yielding pressure vs volume at constant energy; K is the total kinetic energy and U is the potential energy. Temperature is given by Eq.(10) as a function of the energy. Again, the averages in Eq. (17) are computed as time averages along the numerical phase space trajectories. By replacing L = V 1/3 in the potential part of the Hamiltonian (8), and noting that after the replacements xi = Lχi, yi = Lηi and zi = Lζi (where χi, ηi, ζi are adimensional variables) the arguments of the sines do not depend on L, the following simple result is found so that the microcanonical equation of state finally reads Its explicit analytic computation seems hard but still feasible with the aid of the mean-field decoupling of the degrees of freedom in Eq. (14), however such an attempt is beyond our present aims and again we use dynamics in order to estimate the microcanonical averages through time averages. Moreover, by extracting out of U its prefactor V −1/3 [see Eq. (8)], we see that Eq. (19) is of the form where U(E) would be the outcome of the analytic computation of the microcanonical averages, i.e. U(E) = 1 3 V 1/3 U/K K −1 . Equation (20) agrees with the already proposed modification of Boyle's perfect gas law for selfgravitating systems (see the next Section) which is here qualitatively rederived in the microcanonical ensemble and quantitatively computed numerically. The results are reported in Fig.13, where we have differently marked the points that correspond to positive and negative specific heat respectively. The one-to-one correspondence polytrope-negative specific heat (perfect gas law-positive specific heat) is well evident. The transition from a polytropic V −4/3 law to a V −1 perfect gas law is well evident. This reflects the competition between the two terms in the r.h.s. of Eqs. (19) and (20). Figure 14 shows P V /N 2 vs V , obtained at N = 100, 200 and at different values of the energyẼ which have not been shifted towards their corresponding newtonian values as it was done in the preceding Sections (the reason should be clarified by the concluding remarks of the present Section). In order to compare the results obtained at different N ,Ẽ is scaled with N 2 as well as P V . The results of Fig. 14 directly compare to Eq.(20) and make clearer the cross-over between V −4/3 and V −1 . Such a cross-over seems rather smooth, though a discontinuity in the derivative dP/dV cannot be excluded, and this seems to confirm that the gravitational clustering phenomenon is a peculiar phase transition with respect to all the known laboratory phase transitions.
The non-trivial physics behind the cross-over is phenomenologically displayed by Figs. 15 through 19, where a few snapshots of the spatial distribution of the N interacting masses are projected onto the x − y plane. It turns out that the V −4/3 branch of the equation of state corresponds to the clustered phase, whereas the same picture obtained in the V −1 branch displays an homogeneous distribution of particles, typical of a gas phase. Figures 15 and 16 refer to the Fourier-truncated model and are obtained keeping energy constant and varying the volume of the box, whereas Figs. 17 through 19 refer to the direct N -body system and have been obtained by keeping the volume constant and varying the energy.
A remarkable property of this equation of state is the absence of a critical point of flat tangency, i.e. dP/dV = 0, as it occurs in the liquid-gas transition described by the van der Waals equation. The kink in the van der Waals isotherms below the critical one (corresponding to the phase coexistence and to the existence of the latent heat) is here absent, thus the "gravitational condensation" appears rather different from the gas-liquid condensation.
A final comment about the energy values of the P − V curves is in order. At any finite Nw, the Hamiltonian (8) implies a truncation error in the potential energy with respect to the Hamiltonian (1). At fixed volume this implies just a shift in energy scale, making easy the comparison between the Fourier-truncated model and the direct model. But when volume is varied, also the truncation error ∆U varies. It is numerically confirmed that the expected Vdependence of the correction to the potential energy scale is proportional to V −1/3 (for example, at N = 100, Nw = 7, V = 10 2 and V = 2.5 10 5 , we found ∆U/N 5/3 = −0.85 and ∆U/N 5/3 = −0.065 respectively). Thus, the truncation can only shift the value of V at which the crossover occurs, without qualitatively changing the phenomenology. The larger Nw the smaller the shift will be.

COMPARISON WITH THEORETICAL PREDICTIONS
It was suggested by Terletsky (1952) that for a large mass M of gas in a volume V and at temperature T , composed of N particles under a boundary pressure p, the equation of state of the perfect gas (Boyle's law) should be corrected to where α accounts for the shape of the mass and G is the gravitational coupling constant. Shortly after, it was shown that such a correction of Boyle's law results in gravitational instability (Bonnor 1956). Two theoretical milestones followed with the papers of Antonov (1962) and Lynden-Bell & Wood (1968), where a gravitating system of point particles in a spherical box and the thermodynamics of a self-gravitating isothermal gas sphere were respectively considered. Antonov's very interesting result is about the existence of the gravothermal catastrophe when the density contrast ̺c/̺e between centre and edge of the box exceeds the value 709. Lynden-Bell and Wood (1968) related Antonov's result with the existence of negative heat capacity and found the remarkable result that stable isothermal spheres with negative specific heat exist in the density contrast range between 32.2 and 709; the specific heat diverges at ̺c/̺e = 32.2 and becomes positive below this threshold value; above ̺c/̺e = 709 there is a runaway in the catastrophic phase. Moreover, the minimum attainable  temperature by an equilibrium isothermal sphere of radius re and mass M is and is achieved for a density contrast ̺c/̺e = 32.2 where cV becomes negative.  Our present study naturally represents the microscopic counterpart of the thermodynamic framework where all the above mentioned results were obtained.
In the preceding section we have seen that the dynamical-microcanonical equation of state of a selfgravitating N -body system yields a numerical result in perfect agreement with that of Eq.(21), at least as far as the functional form of the equation of state is concerned, and  even though we obtained isoenergetic P − V curves instead of isothermal ones as is implicit in Eq.(21). The transition point, separating Boyle's law from the polytropic one, actually corresponds to the loss of stability of the homogeneous gas phase which goes into a clustered state, as is also qualitatively shown by Fig.s 15 and 16. Moreover, we have found that the transition actually occurs at the minimum attainable temperature in fairly good agreement (taking into We have G = m = k B = 1, re = L and M = N . Trivial algebra gives T min /N 2/3 = (N/V ) 1/3 /2.52, whence T min /N 2/3 ≃ Figure 20. Density contrast vs energy per particle for the Fourier-truncated model with N = 100, Nw = 7 and reflecting boundary conditions. The vertical dashed line separates c V < 0 (left) from c V > 0 (right). The horizontal dashed line indicates the density contrast centre to edge ̺c/̺e = 32. account the difference of the geometry of the box) with the prediction of Eq.(22), and that below the energy that corresponds to the minimum attainable temperature the specific heat is negative, as shown in Section 3.3.2. The increase of the density N/V implies a shift of the transition energy density towards more negative values -as shown by Figs. 3 and 4 -and this is coherent with the Lynden-Bell & Wood prediction of a critical density contrast above which cV < 0 and clustering occurs. In fact, increasing the density amounts to decrease the density contrast (if E/N is kept constant) and increasing E/N amounts to decrease the density contrast (at constant ̺), thus, in order to keep constant a given ratio ̺c/̺e (for example the threshold value) it is necessary to decrease E/N .
Finally, in Fig.20, we show that the density contrast between centre and edge when cV becomes negative is actually in strikingly good agreement with the prediction of Lynden-Bell & Wood. Even if these authors warned that "both these instabilities * * * depend on the presence of the heat-bath; thermally isolated systems do not suffer this type of instability", we found that the clustering instability shows up also in the absence of a heat-bath, whereas this is not the case of Antonov's gravothermal instability.
There are several possible physical explanations for the absence of a dynamical counterpart of the gravothermal catastrophe. A first possibility is that its growth rate is so small that a huge time is required to observe it; this is what has been found for a hybrid model -mean-field with 0.4 for ̺ = 1, whereas we found 0.6 -as shown in Figure 4 -and T min /N 2/3 ≃ 0.03 for ̺ = 4 10 −4 , whereas we found 0.04, as shown in Figure 3. * * * Occurring at ̺c/̺e = 32.2 and at the Antonov threshold ̺c/̺e = 709. local perturbations -where the time scale to develop the gravothermal catastrophe is larger than the Hubble time (Sygnet et al. 1984). Another dynamical mechanism that can inhibit the collapse is the formation of binary systems (Heggie & Aarseth 1992).

CONCLUSIONS
The present paper has mainly dealt with the dynamic ("microscopic") origin of the thermodynamic behaviour of selfgravitating systems and -in close connection with this topic -with the problem of their appropriate statistical mechanical treatment.
The dynamics of self-gravitating N -body systems, described by two different kinds of regularization of the newtonian interaction, has been numerically investigated. Dynamics provides the time averages of certain quantities whichwith the aid of suitable formulae developed in the field of Molecular Dynamics -yield thermodynamics. The link between dynamics and thermodynamics is made in the microcanonical ensemble of statistical mechanics.
The regularizations of the newtonian interaction are the standard softening and a truncation of the Fourier expansion series of the two-body potential. N particles of equal masses are constrained in a finite three dimensional volume.
The introduction of a Fourier-truncated model makes possible a mean-field-like decoupling of the degrees of freedom, which is of great prospective theoretical interest for a statistical mechanical treatment of self-gravitating systems. Moreover, just through such a decoupling, an order parameter can be naturally introduced to signal the occurrence of a phase transition. Actually, through the computation of the caloric curve, of the specific heat and of the order parameter, clear evidence is found for a dynamical signature of a phase transition of clustering type which appears as a second order phase transition. Thus, the gravitational condensation seems to take place trough a transition milder than that of a gas-liquid condensation, which is a first order transition with a finite jump in the entropy. This is also coherent with our numerical computation of the microcanonical equation of state: the counterpart of the clustering transition on the P − V plane is a cross-over between a polytrope of index n = 3, i.e. P = KV −4/3 and a perfect gas law P = K ′ V −1 , without any kink in the P vs V curves as it would be expected for a first order transition.
The very good agreement between direct simulations of the N -body systems and Fourier-truncated models is also interesting because it reveals that the relevant informations for thermodynamics are mainly conveyed by the large spatial scales rather than by the small ones. This is in agreement with the fact that the statistical mechanical peculiarities of N -body self-gravitating systems are due to the long-range nature of the newtonian interaction and are not influenced by its short-scale singularity.
The scale invariance (Heggie & Mathieu 1986) of the gravitational N -body systems is broken by constraining the particles in a box. The physical reason is that the box introduces a new length scale besides the gravitational one. The clustering phase transition is a consequence of the existence of these two independent length scales. Therefore, the physical meaning of what we found is that whenever an extra length scale is present besides the gravitational one, an N -body system can undergo a clustering transition. The use of the box is perhaps the simplest way of breaking scale invariance, but a halo of diffuse matter or any other external potential could in principle work.
We have found that newtonian dynamics naturally yields a number of classical results obtained within a phenomenological (thermodynamic) framework. A remarkable difference exists between the theoretical prediction of two kinds of transitions, the clustering and the gravothermal catastrophe, and our numerical results: only clustering is recovered in the dynamical approach. Actually, the possibility of a dynamical suppression of the gravothermal catastrophe (motivated either by a huge instability time or by the formation of binaries) has been discussed in the more recent literature.
Finally, we have compared the microcanonical (dynamical) averages to their canonical ensemble counterparts, obtained through standard Monte Carlo computations. A remarkable result found here concerns the ensemble inequivalence. This is a-priori expected because of the appearance of negative specific heat -forbidden in the canonical ensemble -in the dynamically worked out thermodynamics. However, the inequivalence is so strong that no signal of the clustering transition seems to survive in the canonical ensemble. In the case of other models some trace is left (Lynden-Bell & Lynden-Bell 1977), but perhaps the transition here is too soft. A reliable statistical mechanical approach to (regularized) self-gravitating systems seems possible only in the microcanonical ensemble. but for abstract systems, there are only trivial examples of physical systems where it is possible to observe a violation of ergodicity (with the exception of amorphous materials, like glasses or spin-glasses), i.e. where the predictions of statistical mechanics are in disagreement with the outcomes of experimental measurements at equilibrium.
The real problem is the dynamical realization of equilibrium. In statistical mechanics it is ab initio assumed that a system is at equilibrium, but the equilibrium concept is very subtle, and depends in a crucial way on the time scales over which a phenomenon is observed. In Feynman's words, at equilibrium all the "fast" things have already happened, while the "slow" ones not yet. This is the reason why from a physical point of view it is mixing, rather than ergodicity, the most important feature of a dynamical system. A mixing dynamics is also ergodic; the converse is not true. Under a mixing dynamics, a generic probability distribution in phase space ̺(q, p, t) converges exponentially fast to the stationary equilibrium distribution. Therefore mixing is responsible for the relaxation to equilibrium, and for the convergence of time averages to statistical equilibrium values on finite time scales. This sheds a new light on the role of dynamics in explaining thermodynamical phenomena, because there is a tight relationship between mixing and that kind of dynamical instability which is called deterministic chaos.
The reason to believe that chaos is at the origin of phase space mixing is that in all the systems where the mixing property can be rigorously ascertained, mixing is found to be a consequence of chaos.

A2 From statistics to thermodynamics
Statistical mechanics also provides a link between the microscopic (dynamic) and macroscopic (thermodynamic) descriptions of large collections of objects. The main frameworks are (Huang 1963) the microcanonical ensemble (when energy and number of particles are given), the canonical ensemble (when temperature and number of particles are given) and grand-canonical ensemble (when total energy and number of particles are allowed to fluctuate). The equivalence of these ensembles, in the large N limit (thermodynamic limit), is a fundamental point, so that the choice of the statistical ensemble is only a matter of practical convenience.
In the microcanonical ensemble the link with thermodynamics is made through the following definition of entropy (V is the physical volume, N is the number of particles) S(E, V, N ) = kB log ωN (E). (A8) Another definition of the microcanonical ensemble volume in phase space, alternative to ωN in Eq. (A5), is where Θ(·) is the Heaviside step function, and the entropy is now given bỹ S(E, V, N ) = kB log ΩN (E).
The two definitions of entropy differ by O(1/N ) terms and therefore coincide in the thermodynamic limit. With both definitions of entropy, temperature is obtained as Strictly speaking, the above definitions of entropy are arbitrary and justified a-posteriori mainly by showing that they are consistent with the laws of thermodynamics (Huang 1963).
In the canonical ensemble the link with thermodynamics is based on the following definition of the Helmoltz free energy where ZN is the partition function defined in Eq. (A7); from the function F all the other thermodynamic functions are obtained through standard Maxwell's relations. Also the definition of F (T, V, N ) is a-priori arbitrary and is given validity a-posteriori. At the macroscopic level, a many-body system has a good thermodynamic behaviour only if the microscopic interaction potential fulfils two basic properties: stability and temperedness (Ruelle 1969). Loosely speaking, temperedness means that the positive part of the interaction energy between particles at large distances is vanishingly small. The negative part of the interaction is left arbitrary by the definition of temperedness, thus it has to be accompanied by the stability condition -which is the relevant condition for gravitational interactions - where B ≥ 0. From these conditions the existence of the thermodynamic limit and the equivalence of the statistical ensembles follow. The existence of the thermodynamic limit means that by letting N → ∞ and V → ∞ so that N/V remains finite, then also the energy density, entropy density and Helmoltz free energy density remain finite. In other words, energy, entropy, free energy are extensive quantities. The ensemble equivalence holds in the thermodynamic limit and means that the same thermodynamics can be described by means of microcanonical and canonical ensembles, for example.

APPENDIX B:
Let us briefly sketch here how the computation of the coefficients of the mean-field Hamiltonian (14) should proceed and, similarly, how this mean-field version of the gravitational N -body Hamiltonian (8)  This integral equation in the unknown S l,m,n is the consistency equation for the mean-field approximation. Now the following identity is very useful S l,m,n = where L stands for Laplace transform and L −1 for its inverse. In fact, the momenta are easily integrated (Pearson et al. 1985) leaving S l,m,n = d xi e −sU i ( x i , S l,m,n ) sin(k l xi) sin(kmyi) sin(knzi) , where C0 is a numerical constant, s is the variable of the Laplace transform, and thanks to the decoupling of the degrees of freedom in H mf , all the triple integrals in square brackets are equal and independent one from another because the potential energy U in Eq. (14) is the sum of independent contributions Ui. If we denote the generic triple integral with F (s, S l,m,n ), we are finally left with S l,m,n = L −1 1 ωN C0 1 s which, in general, can be worked out through some approximate method, the most popular being the saddle point method (Huang 1963).
Once the coefficients S l,m,n have been computed, the mean-field Hamiltonian H mf is completely specified and can be used to compute other microcanonical averages through the same guideline depicted above. This paper has been produced using the Royal Astronomical Society/Blackwell Science L A T E X style file.