Relaxation to thermal equilibrium in the selfgravitating sheet model
Abstract
We revisit the issue of relaxation to thermal equilibrium in the socalled “sheet model”, i.e., particles in one dimension interacting by attractive forces independent of their separation. We show that this relaxation may be very clearly detected and characterized by following the evolution of order parameters defined by appropriately normalized moments of the phase space distribution which probe its entanglement in space and velocity coordinates. For a class of quasistationary states which result from the violent relaxation of rectangular waterbag initial conditions, characterized by their virial ratio , we show that relaxation occurs on a time scale which (i) scales approximately linearly in the particle number , and (ii) shows also a strong dependence on , with quasistationary states from colder initial conditions relaxing much more rapidly. The temporal evolution of the order parameter may be well described by a stretched exponential function. We study finally the correlation of the relaxation times with the amplitude of fluctuations in the relaxing quasistationary states, as well as the relation between temporal and ensemble averages.
M. Joyce and T. Worrakitpoonpon]
M. Joyce and T. Worrakitpoonpon
1 Introduction
The socalled “sheet model” is an interesting toy model for the study of selfgravitating systems, or more generally of systems with longrange interactions. It is simply the one dimensional (1D) generalisation of Newtonian gravity, consisting of particles interacting by attractive forces independent of their separation (or, equivalently, infinite parallel planes embedded in three dimensions interacting via Newtonian gravity). Because the particle trajectories are exactly integrable between crossings, it has the nice feature that its numerical integration can be performed with an accuracy limited only by machine precision. It has been the subject of (mostly numerical) study in the literature for several decades (see, e.g., [1] for a review of the literature on the model) following earlier analytical studies [2, 3]. A fundamental question about this system — and more generally for any system with longrange interactions — is whether they relax to the statistical equilibrium calculated in the microcanonical or canonical ensemble. For this model the latter were first calculated exactly, for any particle number , by Rybicki [4]. The literature on this model — which we will discuss in greater detail in our conclusions section below — is marked by differing results (or, rather, interpretation of results) from different groups, and even some controversy. Work by two groups in the eighties (see, e.g. [5] for a summary) led to the conclusion that relaxation could not be observed, except perhaps in some special cases. Studies by two other groups over a decade ago [6, 7] found results indicating relaxation, and [6] gave a determination of the dependence of the characteristic time. However doubts about the interpretation of these latter results as establishing relaxation to equilibrium were raised by further analysis [8, 9]. In more recent work [10, 1] clear evidence for relaxation in a version of the model in which there are different particle masses has been found, but the dependence on has not been determined^{1}^{1}1Other variants of the model have also been studied in [11, 12, 13]. The mechanism of relaxation (if it indeed takes place) in these models remains, as in other longrange interacting systems, very poorly understood, and a basic subject of research in the statistical mechanics of longrange interacting systems (for recent reviews see e.g. [14, 15]). In this article we report an essentially numerical study of relaxation in the single mass sheet model. We introduce a simple but, as we will see, very useful tool for the characterisation of the longtime evolution and relaxation of the system. This tool allows us to resolve some outstanding issues about the relaxation in this system, and, in particular, to establish more definitively both that relaxation does indeed occur and the scaling with particle number of the time characterizing it. We consider a broader range of initial conditions, which allows us to establish also dependences of relaxation on these. We also study the fluctuations — both in time and over realizations of the initial conditions — about the average macroscopic evolution of the system, showing phenomenologically the correlation of their amplitude with the lifetimes of the intermediate “quasistationary” states.
We will discuss in greater detail in our conclusions the relation of our results to those in the previous literature, but it is useful at the outset to say a little more about the more general context of this study. In recent years there has been considerable interest in the statistical mechanics of long range interactions (see, e.g., [16, 17, 14]), stimulated by the need to understand the physics of various laboratory systems with interactions of this kind, as well as by the more classical case of selfgravitating matter relevant in astrophysics and cosmology. In this context one toy model in particular, and various variants of it, has been much studied: the Hamiltonian Mean Field (HMF) model (see e.g. [18, 19, 20] and references therein), which is a model of particles on a circle interacting by a cosine potential. Its study has shown that it shares many of the qualitative features well documented in the most studied of realistic longrange interacting systems — selfgravitating systems in astrophysics — and believed to be generic in such systems. Starting from generic initial conditions, the system evolves rapidly (by “violent relaxation”) to a virialized macroscopically stationary state. These states — commonly referred to in the more recent literature as “quasistationary states” (QSS) — are out of equilibrium states, which are described theoretically in the framework of Vlasov equation (more usually referred to as the “collisionless Boltzmann equations” in the astrophysical literature). On much longer time scales an evolution towards the true thermal equilibrium (i.e. that determined by the maximization of the Boltzmann entropy in a mean field approximation) is postulated. For realistic systems — such as Newtonian gravity in three dimensions — it is very difficult numerically to simulate the evolution on sufficiently long time scales to probe the relaxation. Studies in the literature (see e.g. [21, 22, 23, 24, 25]) provide some results but give still a very limited characterization and understanding of it. The HMF model has the particular feature that the potential energy of any particle may be expressed as a function of its (angular) position and the mean potential energy due to all particles — it is for this reason that it is “meanfield” — so that the calculation of the forces in a system with particles requires only of order operations (rather than in a typical longrange interacting system). Further the force is continuous at zero separation, so that the difficulties associated in the case of gravity with the regulation of the potential at small scales are avoided. This allows the regime of relaxation to be accessed numerically even for quite large particle numbers. The study of [18] found a scaling of the relaxation time in proportion to (but see also e.g. [19] which finds indications of longer lifetimes for other initial conditions).
It can clearly be of interest to study different toy models, to determine in particular features which are indeed generic. The “sheet model” is probably the oldest toy model of long range interactions — it was first explored in astrophysics as a toy model for selfgravitating systems in three dimensions — and is also, arguably, closer to reality than the HMF which is constrained on a circle. It has, further, as mentioned above the nice feature that it numerical integration can be performed up to machine precision. Despite this, the results concerning its dynamics and relaxation are less clearly determined than for the HMF, and the literature on the subject has, as we have discussed above, been marked by some controversy and results showing that the model has, apparently, some very peculiar behaviours — rapid relaxation to equilibrium for some classes of initial states [26, 27], persistent phase space structures impeding relaxation to QSS [28, 29], macroscopically chaotic behaviour in the long time evolution [30] — which indicate that it might not be a very useful toy model (in that its behaviours are perhaps nongeneric). In this article our main conclusion is that, 1) by using appropriate diagnostics of the macroscopic evolution and 2) by extending simulations to sufficiently large and/or averaging over sufficiently large numbers of realization, one finds behaviour in this toy model very similar in crucial respects to that in the HMF: to a very good first approximation a generic initial configuration relaxes to a longlived QSS, and then relaxes to its statistical equilibrium at sufficiently long time. This latter phase can be characterized apparently by a single timescale, with the evolution of the order parameter during relaxation well fit by a simple function (in our case a better fit is obtained using a simple stretched exponential, rather than a hyperbolic tangent in the HMF as in [18]). On the other hand the dependence of this time scale, linearly proportional to the number of particles , is different to that found in [18] for the HMF. This latter result, however, applies to spatially homogeneous states which, in the HMF, can occur due to the periodicity of the system. Relaxation which is slower than linear in is expected in this case, as shown using analysis of kinetic equations (see, e.g. contributions of P.H. Chavanis, and of F. Bouchet and J. Barré in [17]).
The article is organised as follows. In the next section we recall the basic definitions of the model, and relevant results on its statistical equilibrium. We then explain the choice of the macroscopic parameters (“order parameters”) we choose to monitor the evolution of the system. In the following section we first describe our numerical simulations and the initial conditions we study, and then give our results. In presenting them we give first results for single realizations, and then use temporal averages and finally ensemble averages to derive the scaling with of the relaxation time. This is followed by further study of the fluctuations about the average behaviours of the order parameters. Considering both temporal fluctuations and those in the ensemble, which we show to be very consistent with one another, we observe the correlation between their amplitude in the QSS and the observed relaxation time. In the conclusion sections we return to a more detailed discussion of the previous literature, presenting further results which allow one to understand the reasons for the divergence in conclusions in certain cases.
2 The sheet model
We first recall the model and fix our notation. We next summarize the results of [4] on statistical equilibrium, and then explain the rationale for our choice of “order parameter” in our study.
2.1 Definitions
We consider identical (equal mass) point particles in one dimension interacting by an attractive force independent of separation, i.e., the force on a particle at coordinate position exerted by a particle at the origin is given by
(1) 
where is the coupling. Equivalently it is the pair interaction derived from the pair potential
(2) 
which satisfies the 1D Poisson equation for a point source, (where is the Dirac delta function). Comparing with the three dimensional (3D) Poisson equation shows the equivalence with the case of an infinitely thin plane of infinite extent and surface mass density , which explains the widely used name “sheet model”. We will work in the one dimensional language, referring to “particles”. For a particle at coordinate position on the real axis the total force acting on it is thus
(3) 
where and are, respectively, the number of particles with coordinates greater than or less than (i.e. the force on a given particle is proportional to the difference in the number of particles on its right and its left).
To specify fully the dynamics we must prescribe what happens when two particles arrive at the same point. Since the force is bounded as the separation goes to zero, the natural physical prescription for the 1D model is that the particles simply cross (i.e. pass through one another). In one dimension, however, this is equivalent, up to a change in particle labels, to a hard elastic collision, as such a collision (of equal mass particles) simply results in an exchange of their velocities. Thus, up to particle labels, the sheet model for equal masses is equivalent to one in which particles experience always the same spatially constant force Eq. (3) and simply exchange velocities when they “collide”. As has been noted in some previous studies of models of this kind [31] it is convenient to exploit this equivalence in numerical simulation, as will be described below.
In contrast to Newtonian gravity in three dimensions, the pair potential (2) is positive and diverges at large separations, so that particles cannot escape from the system to infinity. It has therefore no particular interest to enclose the system in a finite box, and indeed such a confinement is not necessary in order to define the statistical equilibrium (in contrast to three dimensions). We will consider therefore always open boundary conditions. Likewise the fact that the potential has no divergence at short distances means that there is no equivalent of the socalled “gravothermal collapse” in three dimensions.
2.2 Thermal equilibrium
It has been shown by Rybicki [4] that the statistical equilibrium for this model can be derived exactly in the microcanonanical ensemble, for any . We will study here the limit (at fixed energy and mass, see [4] for full derivation), in which case the phase space distribution function (i.e. mass per unit phase space volume) becomes:
(4) 
where and are the characteristic scales of velocity and length, and is the total mass of the system. It is straightforward to verify that
(5)  
(6) 
where is the total energy of the system, which allows one to calculate explicitly as a function of and (and ) only.
As is typical of longrange systems, the statistical equilibrium is thus characterised by a space independent Maxwellian velocity distribution and an inhomogeneous spatial distribution. The same solution is recovered in the canonical ensemble. Thus, differently to many longrange systems (including gravity in three dimensions) the two ensembles are completely equivalent. This behaviour is associated also with the absence of microcanonical phase transitions which may arise in such systems.
This equilibrium solution in the continuum limit can be most easily derived by simply maximizing the Boltzmann entropy at fixed mass and energy, using the meanfield expression for the energy:
(7) 
This procedure gives simply
(8) 
where is the mean field potential
(9) 
and the associated mass density profile, which is therefore the solution to
(10) 
It is simple to verify that the expression Eq. (4) results, with appropriate identification of constants and choice of units.
2.3 Order parameters for relaxation
To monitor relaxation to equilibrium it is possible in principle to simply study the full distribution function as a function of time. In practice relaxation is extremely slow (in the characteristic time units of the system) and only accessible numerically for relatively small numbers (of order a thousand) particles, which makes the comparison of the full function subtle (because of finite fluctuations). In the previous literature various methods have been used — statistical tests on the velocity and spatial distributions (e.g. [26, 5]) and analyses based on the evolution of particle energies coarsegrained in time [6]). Here instead we study the evolution primarily using appropriately chosen macroscopic parameters, i.e., “order parameters” which take, in general, distinct values in and out of equilibrium. This is somewhat analogous to the approach used in the study of the HMF, where the magnetisation of the system plays the role of an order parameter used to characterize the evolution out of equilibrium (see e.g. [18]). Once the expected behaviour of the macroscopic parameter is identified, a more detailed analysis involving the distribution functions can be used to confirm that the system has indeed fully relaxed. We will see that, with the choice of parameter we make, it turns out that the single macroscopic parameter is sometimes a better indicator of relaxation that the full density or velocity distributions, and that indeed some of the controversial results in the previous literature may easily be sorted out using the tools used here.
An evident property of the distribution function Eq. (4) is that it is separable in its spatial and velocity coordinates. It is simple to show, as we will now verify, that it is in fact also the unique stationary solution of the Vlasov equation which is separable. Thus if the system is, during its evolution, very close at any time to a stationary solution of the Vlasov equation (which describes the collisionless limit) any parameter probing the degree of separability of the distribution function would be expected to be a useful indicative measure. This leads us to consider order parameters which are simply suitably normalized moments of the distribution function.
Let us first verify the result on separability. The Vlasov equation for the model is
(11) 
where is the mean field acceleration, i.e.,
(12)  
which can be conveniently expressed in terms of the mean field potential satisfying the Poisson equation with as source , i.e.,
(13) 
with an appropriate boundary condition from Eq. (12) [e.g. , where is the total mass]. Seeking a solution which is both stationary and separable we take
(14) 
and thus obtain, on substitution,
(15) 
Given that everywhere (except at the single point which divides the mass in two) we can write this as ^{2}^{2}2Note that this is not true in the HMF model, as the magnetization (which determines the acceleration) can indeed be zero everywhere when the stationary state is spatially uniform. This is a result of the periodic nature of the system. In this case there may thus exist QSS which are separable, uniform in space but with a non maxwellian velocity distribution. Such QSS are indeed observed and have been extensively studied in this model (see e.g. [18]).
(16) 
for any region of where . It follows that both sides are equal to a constant, say, and therefore
(17) 
The right hand side gives
(18) 
Differentiating with respect to and using the Poisson equation for the mean field Eq. (13), this gives exactly the same equation used above to determine the equilibrium solution for . The expression Eq. (4) is indeed therefore the only stationary separable solution of the Vlasov equation.
Given this observation we define the following order parameters:
(19) 
for nonzero and , where
(20) 
and is the value of the parameter for the th particle. By construction these quantities are zero in thermal equilibrium. While a finite number of such moments can of course be zero in a QSS with a nonseparable distribution function, we expect them generically to be nonzero in such states. We will use here both and . As detailed in the next section, we will consider both their temporal evolution in single realizations of our initial conditions, as well as averages of these temporal evolutions. These averages will be performed in two different ways: by averaging over a finite temporal window in a single realization, and by averaging over independent realizations of the initial conditions. Further we will consider the evolution as a function of time of the fluctuations of and with respect to these averages.
3 Numerical simulations
3.1 Algorithm
As remarked above in Sec. 2.1, it is convenient for the numerical integration of the model to exchange particles’ labels when they cross, which is equivalent to treating them as if they undergo an elastic collision in which they exchange their velocities when they meet. The force on each particle is then constant in space and time [and given by Eq. (3)], and the numerical algorithm must simply determine, at any time, the next crossing which occurs, and then exchange the velocities of the “colliding” particles. The optimal way to treat this kind of problem is, as has been pointed out and discussed in detail in [32], by using a socalled “heapbased” algorithm, which uses an object called a “heap” to store in an ordered way the next crossing times of the pairs (see [32] for details). This algorithm requires a number of operations of order to determine the next crossing (rather than of order for the evident direct algorithm in which one calculates and compares directly at each step the next crossing of each of the pairs). Given that the number of crossings per particle per unit time grows in proportion to , the simulation time thus grows in proportion to .
Because the particle trajectories are integrated exactly, the only limit on the accuracy of the numerical integration is thus the numerical precision. As is common practice we will use the total energy (which is conserved in the continuum model) as a control parameter. For the longest simulations we report the error in total energy of the order of .
3.2 Initial conditions
We will consider principally a simple class of spatially uniform initial conditions (IC), generated by randomly distributing the particles on a finite interval. As initial velocity distribution we will consider both the case that initial velocities are zero (“cold IC”) and the case that this distribution is also uniform in a finite interval. The latter are thus random samplings of a particular class of “waterbag” initial conditions in phase space (i.e. in which the phase space density in equal in the region in which it is nonzero), while the cold case can be considered as the limit in which the width of the velocity distribution goes to zero. In Fig. 1 the phase space distribution for a typical IC is shown.
These IC may be characterized solely by the particle number and a single parameter characterizing the waterbag. Rather than the width of the velocity distribution or phase space density, it is convenient to choose the parameter characterizing the waterbag to be the dimensionless initial virial ratio:
(21) 
where is the initial average total kinetic energy and is the initial average total potential energy. By “average” here we mean that the values of and are those calculated for the theoretical waterbag configuration. When we consider, as we do below, different realizations of the initial conditions at fixed and there are of course finite fluctuations about these values of and . The latter correspond then to their average values in the ensemble of initial conditions at fixed and . We note that this ensemble of initial conditions is clearly not a subset of the microcanonical ensemble because there are finite fluctuations also about the average total energy. As these fluctuations are small, however, we will assume below that the evolution of such an ensemble of initial conditions represents well that of an ensemble of such initial conditions with exactly the average energy^{3}^{3}3It is simple to show, given that the particle positions and velocities are both randomly sampled from a PDF which is uniform in a finite interval, that the relative fluctuations in and scale as for large . An exact calculation shows, for example, that at , the normalized variance of , which corresponds to that in the energy for the case , is . This means that the typical amplitude of the fluctuation in the energy for cold IC at is of order the difference between the mean energy of cold IC and IC with ..
We remark on a particularity of the cold IC which we will return to below. In the limit , the evolution from this IC becomes singular at a finite time: an element of mass initially at coordinate position feels a force , where is the initial mass density; all particles are in freefall under a force proportional to their distance, and therefore arrive at the origin at the same time, producing a density singularity. For a finite system, the corresponding behaviour is associated to the existence of a periodic oscillating mode when the particles are initially equally spaced (i.e. on a regular lattice). This “breathing mode” of such a cold system has been discussed in [33]. While there is no such mode in a three dimensional system of a finite number of particles, the limit of a cold spherical initial condition has the analogous singularity (see [34] for a detailed discussion).
3.3 Units and coordinates
For convenience we choose our coordinate system such that the centre of the mass of the system lies at and is at rest (i.e. after distributing the particles as described we add a spatial translation and constant velocity to all particles so that these conditions are satisfied).
We make the following choice of units: we set the particle mass and the coupling to unity, and take . This corresponds to a mass (and particle) density of unity, and a time unit equal to the dynamical time:
(22) 
which is the characteristic time for the system’s evolution under the mean field forces (the mean field forces, of order , moves a system particle of mass over the system size on this timescale). also coincides with the time of the singularity noted above in the smooth limit of the cold IC.
4 Results
The difficulty in this study of relaxation, as in such a study for any longrange system, is that one is interested in studying large systems — so that finite deviations from the meanfield behaviour are small — on a time scale which grows rapidly with (typically, one expects, in proportion to some power of ). Because of numerical limitations, particularly strong because of the computational cost of integrating a longrange interaction, it is in practice often difficult to arrive at definitive conclusions. In the case of gravity in three dimensions, notably, numerical studies exist (see references above) but they give only a very incomplete characterization and understanding of relaxation. As we have discussed in the introduction one of the attractive features of the HMF model is that, because of its mean field nature, the numerical cost of the force calculation is of order , allowing much larger particle numbers — [18] — to be simulated on the relevant long time scales than is feasible in other cases. The principal reason why the early literature on the sheet model was marked by controversy on the question of relaxation is simply, as we will discuss further below, that such relaxation could not be observed on the required time scales for systems sufficiently large so that the finite fluctuations were sufficiently small to allow the clear identification of the average behaviours. The study of [6], taking advantage of the greater numerical resources available already in the nineties, detected relaxation for from specific waterbag configurations and found a scaling of the relaxation time, over a small range in , linear in . This result was obtained, however, by doing a time average of their chosen diagnostic over a very broad time window (of order , only an order of magnitude less than the typical relaxation time for the cases explored), and its solidity has been placed in question in subsequent work [8, 9]. Exploiting the increase in numerical power since then, and aided greatly by the diagnostics we have defined in the previous section, we report here ^{4}^{4}4Evolution of particles to requires about minutes on a single processor; thus, given the scaling with of the computational cost per unit physical time, and a linear growth (see below) in the relaxation time itself, simulation times of order several weeks are required for the most rapidly relaxing case with . Our largest results are ensemble averages over systems with , obtained by running simultaneously on a large number of work stations. results showing relaxation for systems with . Further we obtain our results for the scaling of the relaxation time by doing ensemble averages (over realizations of the initial conditions) without time averages.
4.1 Temporal evolution of order parameters
Shown in Fig. 2 is the evolution of the virial ratio , and the order parameters and in a single realization of a system with , and . Note that the time axis (as it will be invariably here) is logarithmic. In Fig. 3 are plotted the same quantities for . While the fluctuations are very large, particularly in the first case, one can make out that there are, as expected, two stages in the macroscopic evolution probed by these parameters: a first stage () of “violent relaxation” during which all quantities (and notably the virial ratio) fluctuates strongly before settling down to behaviours which appear to fluctuate about a well defined average, and specifically about unity for the virial ratio. The averages of the parameters and are clearly nonzero on a much longer time scale than that characterizing the virialization. These nonzero averages, which appear to be approximately the same in each case for the two different , appear to remain roughly stable until at least about , after which both and start to evolve towards zero. The time scale at which the evolution sets in is clearly significantly shorter for . This behaviour should indicate, as we have discussed above, the relaxation to statistical equilibrium.
These behaviours can be seen more clearly by averaging in a temporal window, of width small compared to the characteristic times scales of this apparent evolution. Shown in Fig. 4 and 5 are the same quantities for the same simulations, but now each point represents the average over one hundred time slices, equally spaced in a window of width centred on the given time.
These behaviours are thus clearly indicative of the evolution expected, which is that believed to be typical of longrange interacting systems: violent relaxation brings one on a short time scale to a QSS, as a result of a mean field dynamics described by the coupled VlasovPoisson equations (and this independent of ). On longer, times scales, one relaxes to the meanfield equilibrium, given in this case by Eq. (4). That the decay to zero of and does indeed correspond to relaxation to the statistical equilibrium of Eq. (4) can be tested in further detail. Fig. 6 shows the velocity and space distributions for and particles, averaged again over a time window of width , at and . The continuous lines correspond to Eq. (4), clearly in very good agreement at the later time, and very different in the QSS phase. We have also checked (but do not show here) the agreement of the distribution of particle energies. These results indicate that and are very good diagnostics of the evolution towards equilibrium: indeed below we will see that they are typically more discriminating of relaxation than the full density and velocity distributions.
4.2 Dependence on initial virial ratio
Shown in Fig. 7 are the evolution of for and starting now from four different values of , as indicated () The results are averaged again in a time window of width . Note that results for , which extend up to , indicate that the evolution towards the statistical equilibrium is really a relaxation to a definitive equilibrium behaviour, i.e., which is stable and persists. This is further confirmed by Fig. 8 which shows the spatial and velocity distributions for the case and at (with the same time averaging window as used above in Fig. 6). Fig. 7 showv that there are, nevertheless, very significant fluctuations in and . These could indicate significant macroscopic, but stochastic, deviations from the equilibrium persisting over very significant times (see [30]). We will present evidence below that they are, as one would expect, finite effects, with an amplitude which decreases as increases.
We observe in Fig. 7 that, as expected, the QSS resulting from violent relaxation is clearly different depending on the initial condition, with very different values of . Further the relaxation toward equilibrium is evident in most cases, but at a time which depends not only on , but also on (or the intermediate QSS state). More specifically, the smaller is the shorter is the lifetime. Indeed for we just see the onset of the relaxation for the case , but do not see it at all for . For there is a difference of a factor of about one hundred in the time at which relaxation appears to becomes clearly visible in the cases and . In the respect we remark that earlier studies have not considered this kind of cold initial condition, in which relaxation occurs more rapidly.
4.3 Estimation of dependence using ensemble average
Let us focus now on the dependence of the relaxation. We wish to determine the scaling with of the characteristic time for relaxation, at a fixed value of the initial virial ratio. Given the very significant noise in the order parameters at the particle numbers we can simulate numerically up to times on which relaxation occurs to do so we must average out these fluctuations. This can be done using either a time average on a single realization (as above) or an average over realizations (or possibly some combination of both).
Shown in Fig. 9 is a plot of the evolution of in three different realizations for and and , up to . The quantities are again averaged in the same window as above. The variance, albeit clearly decreasing with , is in fact still so significant as to make an accurate determination of the scaling difficult. Averaging over larger time windows the curves become smoother, but such differences persist if we use a time window which is small compared to the time scale of the relaxation itself. In short the intrinsic finite fluctuations from realization to realization in the ( dependent) relaxation time are still so large for of this order as to limit significantly the determination of the average behaviour from a single realization.
We thus consider a simple ensemble average, over realizations of the initial conditions. While we could combine time averaging and such an ensemble average, we choose not to do so as this may complicate the interpretation of our result. More precisely, if we perform a time average, we would need to check carefully for any possible dependence of our results on the chosen averaging window. We will explore below in some detail the relation between time averages and ensemble averages over initial conditions.
Shown in Fig. 10 are plots of and averaged over the indicated number of realizations (and without any time average) for each of the indicated particle numbers, for . The error bars in this plot have been estimated by dividing randomly the realizations into two subsamples and recomputing the average in each of them (i.e. the error bar corresponds to the difference in the two averages).
Using these results we now determine the scaling with . Shown in Fig. 11 is a plot of of which , the characteristic time scale for relaxation, as a function of estimated from each of the curves for and . We have determined the value of in each case as that at which the order parameter reaches half its “plateau” value (i.e. in the QSS), i.e., we estimate the value of the parameter which corresponds to the approximate plateau and then determine the time at which half this value is attained. The error bars correspond to those estimated from those given in the previous figure. Shown also are linear behaviours, which in both cases provide a good fit to the results.

100  3.4e+04  8.4e+04  2.7e+05   

200  6.5e+04  1.7e+05  5.3e+05   
400  1.3e+05  6.8e+05     
800  2.5 e+05       
Shown in Fig. 12 are the ensemble averaged evolution of for the three other initial conditions, for the same four values of . The determinations of the relaxation times, for each case where this is possible by the same method as used above, are shown in Table 1. As there are so few points we have not performed the same fitting procedure (with estimated error bars) as above, but it is clear that the given values are consistent with a scaling of the relaxation time linear in . In the case , however, we cannot deduce any reliable estimate of the scaling with , as we can just see the onset of the relaxation for but not in the other cases.
This last curve and the data in Table 1 allow us to see more quantitatively the dependence of the relaxation time on the initial value of also. At fixed we see that, between and the estimated relaxation time increases by a factor of about eight. These considerable differences translate into a very different appearance to the curves: in the case of the “QSS plateau” is much less visible as there is only a very small separation between the time scales for the establishment of the QSS () and the onset of relaxation.
The exact definition taken here for the relaxation time is somewhat arbitrary — we could equally consider the time at which deviates by from its plateau value, or, say, reaches of this value. Because the relaxation is very slow — to show the evolution of we must plot it as a function of the logarithm of time — such definitions would give enormously different results for the estimated time (differing by two to three orders of magnitude). Equally we see from Fig. 11 that if we use rather than , employing the same criterion we obtain times differing by an order of magnitude. That this factor indeed changes only the overall normalization of the times, and not their scaling with , is evident from the fact that, as can be seen by eye, the curves in the decay phase can be superimposed on one another well by a translation parallel to the time axis.
It is interesting to see if a simple functional behaviour can be fit to the decay of the order parameters. Shown in Fig. 13 are best fits to two simple functions for the case of initial conditions and , for which we have the best statistics. We have restricted to the range to cut out the initial (violent) relaxation phase. One employs a hyperbolic tangent given by
(23) 
in which, therefore, corresponds to the time estimated above. The best fit values of the parameters are , and . The other is a stretched exponential form:
(24) 
and gives the bestfit values and . The second function is clearly a significantly better fit. We note that the former function has been shown in [18] to give a good fit to the temporal evolution of the magnetisation during relaxation in the HMF model. Stretched exponential relaxation, on the other hand, is observed in a range of physical systems, and notably in the relaxation of structural and spin glasses (see, e.g., [35]).
We draw attention to one important feature of these results which introduces a systematic uncertainty into them, which could only be reduced by doing significantly larger simulations: in principle the intermediate QSS is independent of the number of particles , i.e., we are estimating the dependence of the relaxation time of a state which is, up to fluctuations, independent; in practice it is clear in our data that there is some residual dependence in the QSS at the we are simulating — the “plateau” in the curves of the time evolution of our order parameters do not exactly coincide. As we have seen that there is clearly a significant dependence of the lifetime on , which we interpret to be one on the intermediate QSS rather than the initial condition itself, it is possible that the dependence we measure at fixed is due to, or partially due to, this residual dependence of the QSS. We believe, however, that such an effect, if present, is probably negligible: the differences in the QSS “plateau” at given for different are very small compared to the differences between the QSS over the range of , and further, for the larger , the QSS do appear to converge. This is even the case for , where the dependence in the “plateau” is most evident. In this case, as we mentioned above when we discussed our initial conditions, an intrinsic dependence of the QSS might be anticipated: as the evolution becomes singular at , and the evolution at finite is regulated by the fluctuations about a uniform distribution which are dependent^{5}^{5}5For the analogous 3D problem — evolution from cold uniform initial conditions — the precise dependence of the virialized QSS state has been determined numerically in [34].. That such intrinsic dependence is weak, if present at all, is also indicated by the absence of visible dependence of in Fig. 10.
4.4 Relaxation and fluctuations in the QSS
Analytically the relaxation towards equilibrium of systems with long range interactions may be described by kinetic equations, derived for example from the BBGKY hierarchy. In practice these equations are intractable, and despite many attempts to develop appropriate approximation schemes which might make them tractable, there are really no solid results which allow us rigorously to model analytically the detailed phenomenology of relaxation observed in numerical simulations, and determine for example the observed dependence of the relaxation time.
Inspection of our results for the temporal evolution of the parameters and lead to one simple observation: the relaxation time appears to be correlated with the amplitude of the fluctuations about the relevant QSS, i.e., the smaller the fluctuations in the QSS, the longer is its lifetime. While this is somewhat trivial when we consider a given initial condition (i.e. ) at fixed — in postulating that there is a QSS we mean that the fluctuations about it are dependent (and decaying with ) — it is not evident that this should be so for the different at fixed . Theoretically such a correlation might not be surprising — in kinetic theory approaches the leading corrections to the collisionless (Vlasov) limit are, in perturbative approaches, sourced by fluctuations about the QSS (see, e.g., contributions of P.H. Chavanis, and of F. Bouchet and J. Barré in [17].).
Such a trend can be seen a little in Fig. 7, although in this case it is greatly obscured by the time averaging (i.e. it is much clearer if one plots a single realization in each case, which we have not done here). It is shown clearly to be present by the results in Figs. 14 and 15. The first shows the standard deviation, , of as a function of time, estimated in the indicated number of realizations of initial conditions , for each of the different values of indicated. The error bars in the plot correspond to the spread in when it is estimated in two subensembles defined by randomly dividing the realizations into two. As remarked above the fact that decreases with — and thus, given that the lifetimes of the states have been observed to increase with , that there is a correlation of the lifetime with their amplitude — is not surprising: it simply means that the fluctuations about the QSS are, predominantly, due to finite effects which will vanish as . We note that the amplitude of in the approximate “plateau” region — corresponding to the QSS – scales as , i.e., as they would if the fluctuations of is the sum of uncorrelated contributions from the particles.
Shown in Fig. 15 is the same quantity but now for the different values of , at two different fixed ( and ). In both cases we see clearly (except perhaps for the lowest curve in the lower panel, which is noisier due to the much smaller number of realizations) that the amplitude of fluctuations decreases as increases, i.e., that the amplitude is (inversely) correlated with the lifetimes we have observed for these states.
We note that these figures giving the behaviour of the variance of our macroscopic parameters also contain alot of other useful information beyond the correlation we have just observed. Indeed these curves themselves show very clearly the different timescales in the dynamics: the first period of “violent relaxation” is clearly identifiable by a very large variance, which decays on a time scale of order several tens of dynamical times; this is followed by an approximately stable value depending on the QSS (Fig. 15), which then evolves on a much longer time scale, dependent on and increasing with , towards a value which is independent of the initial state (i.e. thermal equilibrium).
These results also allow us to conclude more about the meaning of our quantitative results for the scaling of the relaxation time, which have been calculated using the ensemble average: the systematic decrease with of the variance of (and, we have verified, of ) implies that such an ensemble average, for sufficiently large , can indeed be interpreted consistently to give the macroscopic behaviour of a single realization in the ensemble. Thus, as we have been implicitly assuming, we can indeed take our determined relaxation times to represent those of single realizations, at sufficiently large .
It is interesting to go a little further and consider what the relation is, at finite (but large) , between the fluctuations in the ensemble average and the temporal fluctuations in a single realization. Indeed if, as we have postulated above, there is a real correlation between the amplitude of the fluctuations measured in the ensemble average and the lifetime of the corresponding QSS, it must be that these fluctuations measured in the ensemble bear some close relation to the temporal fluctuations in the same parameters in a single realization. That this is the case, to a good first approximation, can be seen from Figs. 16 and 17, which show exactly the same quantities as in the previous two figures, but that the standard deviations are calculated in one hundred time slices equally spread over a time window of width centred on the indicated point. We see the same quantitative behaviours as in the previous plots, and even, in particular for , quite good qualitative agreement.
To allow further more detailed comparison, in Fig. 18 and Fig. 19 are shown, for (left panels) and (right panels), and in both cases, the histogram of the values of , at the indicated times measured, in Fig. 18, in one hundred simulations from realizations of the same initial conditions, and, in Fig. 18, in one hundred snapshots in a window of width in a single simulation from one realization of the same initial conditions. Comparing the fours panels in the two figures one by one, we see that, although the fluctuations in each case are clearly not sampled from an identical distribution, the agreement is strikingly good: not only, as expected from what we have already seen above, do the averages and variances agree well in each case, but the general shape of the histograms, which is quite different in each QSS, resemble one another strongly. The results are also clearly in line with the conclusion drawn above for what concerns the relaxation to thermal equilibrium: at we see that the cold initial conditions have relaxed to a distribution centered on the value in thermal equilibrium, while for the case the system is still in a QSS but has evolved very slightly towards equilibrium.
It would be interesting to develop this study with greater statistics, varying the width of the time window to see how good agreement can be obtained, but we will not do so here. Such a study is related to the fundamental (and so far unanswered) question as to whether the properties of QSS may be determined by averaging over an appropriate (nonequilibrium) ensemble, determined by the initial conditions. The theory of violent relaxation formulated by LyndenBell, for example, postulates an answer to this question [36]: the appropriate ensemble is that of all configurations corresponding to a phase space distribution function permitted by the (collisionless) Vlasov dynamics. If this theory were correct (which is not the case for this system [37]) we should perform our ensemble average over such configurations rather than the more restricted one we have considered.
5 Conclusions and discussion
5.1 Summary
Our primary aim in this paper has been to establish and characterize more fully than in the previous literature the relaxation to thermodynamic equilibrium of one of the simplest toy models for longrange interacting systems: equal mass selfgravitating particles in one dimension (or infinite sheets in three dimensions). Compared to the much studied HMF model, notably, the basic properties of this model have remained somewhat unclear, and indeed marked by some controversy in the literature. The novelty of our work compared to previous studies is not just that we do more and larger simulations from a broader range of initial conditions, but that we have identified a tool which is very useful to characterize the evolution of the system: the measurement of appropriately normalized moments of the distribution function which characterize the “entanglement” of the one particle distribution function in configuration and velocity space. This is particularly appropriate as a measure simply because the thermal equilibrium has the property that such entanglement is absent while, we have shown, in any other stationary solution of the VlasovPoisson equations such entanglement is present. We note that this result, which we showed to be valid for any interaction in one dimension (but, as noted, excluding periodic systems like the HMF), can be generalized easily to three dimensions if we restrict to stationary solutions which have radial symmetry in space and velocity. This suggests that these “order parameters” may also be useful indicators of relaxation in much more general, and perhaps, as we discuss below, even useful tools for understanding the mechanisms of such relaxation.
With the aid of these macroscopic measures, we have shown in our numerical study, of a range of simple “waterbag” and cold initial conditions, that the system manifests the behaviour thought to be generic in longrange systems: there are essentially two phases in the evolution with two completely different timescales. An initially nonstationary state evolves first, on timescales characterized by the “dynamical time” (roughly the crossing time of a particle in the meanfield due to the others), to a QSS, an out of equilibrium state, which then evolves on a much longer time scale, dependent on the number of particles, to thermal equilibrium. In other words it is reasonable to suppose that the system is ergodic (and mixing) on these very long time scales, but not so on the shorter time scales. Further we can identify clearly that the QSS resulting from different initial conditions (i.e. different values of ) are very different macroscopically, characterized by very different phase space entanglement.
Focussing on the the dependence of the relaxation, averaging over a very large number of realizations to average out the fluctuations, we have concluded that the characteristic time scale for relaxation behaves, to a very good approximation, as
(25) 
where is a numerical factor which depends on the initial condition, which we have denoted in this way as we expect that this dependence is not strictly on the initial condition but on the QSS which results from it. We have seen that this prefactor increases as does, by about a factor of ten between and , and approximately a further factor of ten for . We have noted that the overall normalization of is rather arbitrary, as it depends greatly on the exact criterion used to define the relaxation timescale. Given that the evolution towards zero of , which is what we have used to determine this time scale, is in fact well fit by the simple functional behaviour as a function of the time on a logarithmic scale, the normalisation of can differ by two orders of magnitude by a trivial change in its definition. More specifically we have seen, that in the case where we have accumulated the greatest statistics allowing us to constrain the temporal evolution, a very good fit to our order parameter is obtained to a stretched exponential form.
Although the relaxation of this system, and in general longrange interacting systems, is not well understood, we can say that this finding of a linear scaling — besides the fact that it is, as we will discuss below, in line with less complete previous analyses — is not a surprising result: such a scaling can be anticipated both on the basis of simple naive estimates of the effects of collisionality, as well from general considerations based on kinetic theory.
A simple “derivation” of this scaling, along the lines of those applied originally by Chandrasekhar (see [38] or [39]) to obtain such an estimate for 3D selfgravitating systems, may be given as follows. Relaxation is in principle due to the discretisation, in particles, of a continuous mass distribution. Let us suppose that this latter field density varies spatially on a scale, . The typical fractional change in the velocity of a test particle due to its interaction with one (localized) particle, rather than the continuous mass distribution, can be estimated as . As it crosses the system (in a time ) such a particle will interact with all particles. Assuming the contribution from each particle adds incoherently, one estimates
(26) 
for the normalized variance of the velocity in . Scaling with at fixed mass and energy (and fixed ) requires , and therefore . It follows that the relaxation time scales linearly with in units of the . A slightly different argument leading to the same result may be found in [6], and a more developed analysis in [40]. In the framework of approaches based on kinetic theory, a linear scaling is obtained as collisional terms arise as the leading corrections in a formal expansion in which gives the collisionless (Vlasov) limit at leading order (see, e.g. [41, 42, 43]).
This scaling linear in is to be contrasted with the case of the scaling observed for the lifetime in QSS in the HMF, proportional to . While the exponent found is not understood, the fact that it is larger than unity is consistent with these considerations as this result applies for spatially homogeneous QSS (which are possible in the HMF because of its periodicity) for which it has been shown that the leading collisional term vanishes (see, e.g., contributions of P.H. Chavanis, and of F. Bouchet and J. Barré in [17]).
We note that our study suggests also that the “order parameters” we have defined and studied may be relevant quantities for understanding relaxation in this and other longrange systems. Indeed in all cases we have observed that, at sufficiently large , these parameters start from a nonzero value in the initial QSS and evolve monotonically towards zero, i.e., the relaxation of the QSS can apparently be described as a process of progressive “disentanglement” of the one particle phase space density. In this respect the very different, much less efficient, relaxation observed in the HMF might be interpreted as a result of the absence of such entanglement in spatially uniform QSS. Further, in the case where we have enough statistics to provide a precise fit to the evolution of the parameter , we found that it is well fit by a simple stretched exponential form. It would be interesting to see in further study whether this fit is more than an approximate numerical fit for the case we have studied, and, if so, whether the exponent characterizing it is the same or not. As we have remarked such a functional form has been observed in other slowly relaxing (e.g. glassy) systems and theoretical tools derived in this context to understand relaxation may be relevant. In [35], for example, this behaviour is linked to the existence of a fractal structure in a bounded accessible region of phase space.
5.2 Comparison with previous literature
Let us now turn finally to compare our findings in greater detail with those in the previous literature.
An early numerical study by Hohl and Broaddus [44] which concluded a relaxation time proportional to was found to be incorrect by two groups, who studied the problem in greater detail (and with greater numbers of particles). However, these groups found conflicting results: Miller et al. found no evidence for relaxation at all to thermal equiibrium in their simulations [45], while Luwel et al. [26] found relaxation on a time scales even shorter than . Further study (see [5, 27], which also contain a detailed synthesis of the previous literature) by Miller et al. concluded that the discrepancy was related to the very specific initial condition studied by the other group. Studying this case in detail they found that it indeed appears to thermalize very rapidly, but some further, but not completely conclusive analysis of the evolution at longer times, suggested that this thermalization was not complete.
To determine whether these cases are consistent with our findings — and see whether our analysis using the parameters and can throw light on these previous findings — we have resimulated the relevant initial conditions. These are “counterstreamed” waterbag initial conditions, an example of which is shown in Fig. 20. We have simulated a range of such initial conditions, in particular the cases (one of which is that shown in the figure) considered by [26] and [5]. Shown in Fig. 21 are the density profile and velocity distribution at obtained starting from a realization of initial conditions like those shown in Fig. 20, but with . We see that the profiles indeed agree very well with the equilibrium ones. In Fig. 22 is shown the evolution of as a function of time for the indicated values of averaged in each case over the number of realizations indicated. We observe that, although small and fluctuating, its value is clearly on average nonzero, indicating that the state, despite the good agreement of the density and velocity profiles, is not in fact an equilibrium. Just as in the cases we studied we see clearly the relaxation towards equilibrium at longer times, and indeed that the characteristic time increases on . Although we haven’t done the more extensive study required to determine precisely this dependence, the results are quite consistent with Eq. (25) with a value of of order that found for the case .
This case illustrates the usefulness of the parameters and as discriminants of relaxation: indeed we have just seen that the single measure of is sufficient to discard the interpretation of Luwel et al. [26] of their results. This is simply because they are physically very appropriate indicators, for the reasons we have explained in introducing them: the property they probe — of entanglement of the phase space distribution — is one which must evolve significantly during relaxation, because the phase distribution must become separable. While and being zero does not imply thermalization, of course, we have not found a single QSS, despite exploring a broad range of initial conditions (considerably more extended that those reported here) in which they are both zero (within the uncertainty of fluctuations), i.e., the only states we have found in which they are both zero are states which we have concluded, using a range of other measures, are indistinguishable from the equilbrium state of Rybicki. It is not difficult, on the other hand, to find initial conditions which lead to a QSS in which or . Indeed for the waterbag initial conditions we have studied both and actually change sign as varies over the range we have considered, and one can thus find by trial and error the appropriate which make them zero individually.
Another evident quantity to measure, which we have in fact considered systematically but have not reported in detail, is the distribution of the individual particle energy . This is in fact generally a better discriminant of relaxation than either and , i.e., we have found that in quite alot of cases and are not easy to distinguish from the equilibrium profiles, but that allows one to see more clearly that one is indeed not in the equilibrium state. An example is the counterstreamed case just considered above. In Figs. 23 are shown, for , the evolution of the ensemble averaged at a few different times. We have not plotted the equiibrium curve, as it is indistinguishable from the measured curve at the final time shown. One can see clearly see that, despite the good agreement of and shown in Fig. 21, the system is not in equilibrium at the early times: has a clearly visible excess of particles at high energies compared to that at the much later times at which the evolution of indicated relaxation (and indeed approaches very accurately its predicted equilibrium form). While such a measure of , averaged over a large ensemble of realizations, can, in all the cases for which we’ve studied it, clearly discriminate relaxation, the use of just (and possibly ) is an extremely efficient shortcut to “diagnose” relaxation.
Subsequent to [5], in the nineties, Tsuchiya et al. reported an analysis of larger and more importantly, longer, simulations in order to clarify the issue. A first paper [6] they reported the evolution of a rectangular waterbag initial condition corresponding to our case , and reported a detection of relaxation to thermal equilibrium. These authors made a distinction between two time scales of relaxation: one of “microscopic relaxation”, the other for “macroscopic relaxation”. These are identified, and both found to be proportional to , by considering the evolution of the mean standard deviation of the particle energies averaged over a time window from their equipartition value. The former is estimated from the slope at short time of this function, and the latter from the position of “peaks” which are observed to occur at much longer times. While the latter is interpreted in terms of of macroscopic relaxation in the sense we have used here, the former is interpreted as a time scale on which particles sample the energy distribution but on which there is no macroscopic evolution. The justification for these interpretations are not clear, and no direct comparison with the equilibrium distribution derived by Rybicki, Eq. (4), has been reported which might show their correctness. Indeed both a subsequent article by the same authors [8] and a study by Yawn and Miller [9] place in doubt the correctness of the interpretation in terms of relaxation.
Nevertheless, in light of the results we have given here, it would be reasonable to infer that the results given by Tsuchiya et al. in [6] are indeed correct, at least for what concerns the dependence of the relaxation. Further comparison could of course clarify the relation of the behaviour of their measured quantities and the macroscopic relaxation as we have probed it here (and should be much easier for the shorter lived, smaller , initial conditions rather than the case studied by these authors). We do not believe, however, that there is any clear basis for either an operational or physical distinction between “microscopic” and “macroscopic” relaxation as described by these authors: as we have discussed there is an arbitrariness in the definition of the relaxation time because of the very slow nature of this relaxation. As we have noted, we could easily, for example, have obtained here estimates of the relaxation time differing by several of orders in magnitude in their prefactor, just like the two different time scales determined by Tsuchiya et al. [6], by using slightly different definitions, or choosing to use a different order parameter. This point can be illustrated by considering the evolution of for one of the cases we have considered: shown in Figs. 24 is this quantity for the case and , averaged over realizations. While we have associated (see Table 1 above) the time scale to the relaxation in our analysis, one can discern by inspection of these figures significant evolution (in particular of the initially clear “corehalo” structure) in already by , i.e., there is evolution of the energy distribution on the time scale of “microscopic” relaxation (of order ) identified in [6]. While it is possible that there are different time scales associated to different physical processes as argued in [6], it seems a more plausible interpretation to us to suppose that there is single physical relaxation process leading, albeit very slowly, to macroscopic relaxation of the system, and to characterize this relaxation by a function and the scaling of its parameters with . In this respect it is interesting to note that the specific stretched exponential form we fitted to the temporal behaviour has the known property [46] that it is can be written as a weighted integral over simple exponentials (i.e. it can be interpreted as arising from the superposition of an infinite number of relaxation processes each with a single characteristic time).
In [8] Tsuchiya et al. have also described chaotic “itinerant” behaviour of these systems, starting from the same () initial conditions i.e., in which the system shows apparently stochastic macroscopic behaviour. In our analysis this would correspond to such behaviour for the parameters or . While we have seen that there are indeed very significant fluctuations in these parameters, which correspond to very significant differences in the “macroscopic” evolution of these systems, we have studied carefully their dependence on and found them to decay monotonically. The results of [8] were obtained for , a range in which we still see fluctuations of which are order unity. Only when we reach of order several hundred do we see these fluctuations diminish significantly so that the macroscopic trajectory of the system becomes quite localized. We thus believe that as increases these effects will becomes negligible, even on the time scales on which relaxation occurs, and an effectively deterministic macroscopic evolution will occur.
It is interesting to compare our results also to those of Yawn and Miller [1], who have analyzed in detail relaxation in a version of the sheet model in which there are sheets of different masses. In this case the relaxation towards thermal equilbrium may be clearly distinguished by testing for equipartition of the kinetic energy, and the associated spatial segregation of the mass populations. In simulations starting from waterbag type initial conditions with a virial ratio of two, for a range of different mass ratios and up to particles, clear evidence was found in [1] for such relaxation using appropriately defined indicators. Like the order parameters we have employed here, these show characteristic behaviours corresponding to the principal phases of the dynamical evolution (violent relaxation, QSS phase, relaxation towards thermal equilibrium). Although we cannot compare our results directly, we note that the time scales observed for relaxation of systems with particles are quite consistent with those we have observed for the equal mass system with initial virial ratio . Yawn and Miller [1] also measure temporal correlation properties and find weak but persisting correlations characterized by a powerlaw decay (in time), which they interpret as evidence for the incompleteness of relaxation. In the present study we have found, in contrast, that our principal observables decay in time with a functional form which allows the identification of characteristic time scales. Further all deviations of these observables from their equilibrium values decrease clearly as increases, and thus we have interpreted the associated “incompleteness” of relaxation simply in terms of finite effects. It would be interesting certainly to perform a more direct comparison of the results in the two models, and in particular to extend the study of Yawn and Miller to allow a determination of the dependence of the parameters they study. We note also that Yawn and Miller argue that the powerlaw decay suggests the existence of a fractal structure in phasespace, which, as we have been mentioned above, is also proposed as an explanation in [35] for the appearance of relaxation characterized by a stretched exponential behaviour.
Let us finally mention some other issues of importance concerning aspects of the dynamics of this system which have been treated elsewhere but which we have not discussed here. As we have discussed, we interpret our results in line with those of many previous studies of this and other longrange systems: the evolution from an arbitrary out of equilibrium initial condition is characterized a first phase of relaxation to a QSS, interpreted as a finite particle sampling of a stationary solution of the Vlasov equation, on a time scale independent of , followed by a slow relaxation to thermal equilibrium on an dependent time scale. Studies of the single mass sheet model for other specific initial conditions suggest that this simple scheme may be too limiting, for this model (and possibly, for all such models). On the one hand Reidl and Miller have reported numerical results [47] for specific “two cluster” initial conditions which show a dependence on in the time scale for relaxation to a QSS. On the other hand, as mentioned in the introduction, Rouet et al. [28, 29] have shown, using both particle simulations and simulations of the Vlasov equation, for yet other initial conditions that “holes” which rotate in phase space may be present after violent relaxation and persist on very long time scales. Although it is not evident that there is necessarily a relation between either finding and the mechanism of relaxation to thermal equilibrium, a study incorporating such initial conditions would certainly be more complete that that reported here. Extension of the study reported here to larger still would likewise be desirable, despite the extremely rapidly growing numerical cost of such simulations with .
5.3 Acknowledgements
The simulations were carried out in large part at the Centre de Calcul of the Institut de Physique Nucléiare et Physique des Particules. We are particularly grateful to Laurent Le Guillou for advice and help. We thank also Duccio Fanelli for providing us with his own code which allowed us to perform checks of ours. We thank P. Astier, J. Barré, A. Gabrielli, B. Marcos, P. Viot, F. Sicard, F. Sylos Labini for useful conversations, comments or suggestions.
References
 [1] K. R. Yawn and B. N. Miller. Incomplete relaxation in a twomass onedimensional selfgravitating system. Phys. Rev. E, 68(5):056120, Nov 2003.
 [2] G.L. Camm. Selfgravitating star systems. Mon. Not. R. Astron. Soc., 110:305, 1950.
 [3] J. H. Oort. The force exerted by the stellar system in the direction perpendicular to the galactic plane and some related problems. Bull. Astron. Inst. Neth., 6:249–+, August 1932.
 [4] G. Rybicki. Exact statistical mechanics of a one dimensional selfgravitating system. Astrophys. Sp. Sci., 14:56, 1971.
 [5] C. Reidl and B.N. Miller. Gravity in onedimension: selective relaxation. Astrophys. J., 318:248, 1987.
 [6] T. Tsuchiya, N. Gouda, and T. Konishi. Relaxation processes in a onedimensional selfgravitating manybody system. Phys. Rev., E53:2210, 1996.
 [7] Lj. Milanović, H.A. Posch, and W. Thirring. Statistical mechanics and computer simulations of systems with attractive power law potentials. Phys. Rev., E57:2763, 1998.
 [8] T. Tsuchiya, N. Gouda, and T. Konishi. Chaotic Itinerancy and Thermalization in a OneDimensional SelfGravitating System. Astrophy. Space Sci., 257:319–341, 1997.
 [9] K. R. Yawn and B. N. Miller. Ergodic properties and equilibrium of onedimensional selfgravitating systems. Phys. Rev. E, 56(3):2429–2436, Sep 1997.
 [10] K. Yawn and B. N. Miller. Equipartition and mass segregation in a one dimensional gravitating system. Phys. Rev. Lett., 79:3561, 1997.
 [11] B. N. Miller and P. Youngkins. Phase Transition in a Model Gravitating System. Physical Review Letters, 81:4794–4797, November 1998.
 [12] P. Valageas. Thermodynamics and dynamics of a 1d gravitational system. Astron. Astrophys., 450:450, 2006.
 [13] P. Valageas. Relaxation of a 1d gravitational system. Phys. Rev.., E74:016606, 2006.
 [14] A. Campa, T. Dauxois, and S. Ruffo. Statistical mechanics and dynamics of solvable models with longrange interactions. Phys. Reports, 480:57–159, September 2009.
 [15] F. Bouchet, S. Gupta, and D. Mukamel. Thermodynamics and dynamics of systems with longrange interactions. ArXiv: 1001.1479, January 2010.
 [16] T. Dauxois, S. Ruffo, E. Arimondo, and M. Wilkens. Dynamics and Thermodynamics of Systems with Long Range Interactions. Springer, Berlin, 2002.
 [17] A. Campa, A. Giansanti, G. Morigi, and F. Sylos Labini. Dynamics and Thermodynamics of Systems with Long Range Interactions: Theory and experiments. AIP Conference Proceedings, 2008.
 [18] Y. Y Yamaguchi, J. Barré, F. Bouchet, T. Dauxois, and S. Ruffo. Stability criteria of the vlasov equation and quasistationary states of the hmf model. Physica A, 337:36–66, 2004.
 [19] A. Campa, A. Giansanti, and G. Morelli. Longtime behaviour of quasistationary states of the Hamiltonian meanfield model. Phys. Rev. E, 76(4):041117–+, October 2007.
 [20] P. H Chavanis and A. Campa. Inhomogeneous Tsallis distributions in the HMF model. ArXiv: 1001.2109, January 2010.
 [21] Ch. Theis and R. Spurzem. On the evolution of shape in nbody simulations. Astron. Astrophys., 341:361–370, 1999.
 [22] J. Diemand, B. Moore, J. Stadel, and S. Kazantzidis. Twobody relaxation in cold dark matter simulations. Mon. Not. R. Astron. Soc., 348:977–986, March 2004.
 [23] J. Binney and A. Knebe. TwoBody Relaxation in Cosmological Simulations. Mon. Not. Roy. Astron. Soc., 325:845, 2001.
 [24] A. El Zant. Twobody relaxation in simulated cosmological haloes. Mon. Not. Roy. Astron. Soc., 370:1247, 2006.
 [25] Y. Levin, R. Pakter, and F. Rizzato. Collisionless relaxation in gravitational systems: from violent relaxation to gravothermal collapse. Phys. Rev., E78:021130, 2008.
 [26] M. Luwel, G. Severne, and P.J. Rousseeuw. Numerical study of the relaxation of one dimensional gravitational systems. Astrophys. Sp. Sci., 100:261, 1984.
 [27] C. Reidl and B. N. Miller. Gravity in onedimension: a correction for ensemble averaging. Astrophys. J., 371:371, 1991.
 [28] J.L. Rouet and M.R. Feix. Persistence of collective fluctuations in nbody metaequilibrium gravitating and plasma systems. Phys. Rev., E59:73, 1999.
 [29] P. Mineau and M .R. Feix and J. L. Rouet. Numerical simulations of violent relaxation and formation of phase space holes in gravitational systems Astron. Astrophys., 228:344349, 1990
 [30] T. Tsuchiya and N. Gouda. Relaxation and lyapunov time scale in a onedimensional gravitating sheet system. Phys. Rev., E61:948, 2000.
 [31] E. Aurell, D. Fanelli, and P. MuratoreGinanneschi. On the dynamics of a selfgravitating medium with random and nonrandom initial conditions. Physica D, 148:272–288, 2001.
 [32] A. Noullez, E. Aurell, and D. Fanelli. A heapbased algorithm for the study of one dimensional particle systems. J. Comp. Phys., 186:697–703, 2003.
 [33] C. Reidl and B.N. Miller. Gravity in one dimension: The critical population Phys. Rev., E48:42504256, 1993.
 [34] M. Joyce, B. Marcos, and F. Sylos Labini. Energy ejection in the collapse of a cold spherical selfgravitating cloud. Mon. Not. R. Astron. Soc., 397:775, 2009.
 [35] R.M.C. Del Almeida, N. Lemke, P. Jund, R. Julien, I. A. Campbell, and D. Bertrand. Dynamics of complex systems above the glass temperature. Jour. NonCryst. Sol., 287:201, 2001.
 [36] D. LyndenBell. Statistical mechanics of violent relaxation in stellar systems. Mon. Not. R. Astr. Soc., 167:101–121, 1967.
 [37] Y.Y. Yamaguchi. Onedimensional selfgravitating sheet model and lyndenbell statistics. Phys. Rev. E, 78:1114, 2008.
 [38] S. Chandrasekhar. Stochastic Problems in Physics and Astronomy. Reviews of Modern Physics, 15:1–89, January 1943.
 [39] J. Binney and S. Tremaine. Galactic Dynamics. Princeton University Press, 1994.
 [40] B. N. Miller. Source of relaxation in a one dimensional gravitating system. Phys. Rev, E53:R4279, 1996.
 [41] R. Balescu. Equilibrium and nonequilibrium statistical mechanics. Wiley, New York, 1975.
 [42] P.H. Chavanis. Hamiltonian and Brownian systems with longrange interactions: II. Kinetic equations and stability analysis. Physica A, 361:81–123, February 2006.
 [43] P.H. Chavanis. Kinetic equations for systems with longrange interactions: a unified description. J. Stat. Mech. 5:19, 2010
 [44] F. Hohl and T. Broaddus. Thermalization effects in a onedimensional selfgravitating system. Phys. Lett., A25:713, 1967.
 [45] H. L. Wright, B. N. Miller, and W. E. Stein. The relaxation time of a onedimensional selfgravitating system. Astrophys. Space Sci., 84:421–429, June 1982.
 [46] E. Montroll and D. Bendler. On Levy (or stable) distributions and the WilliamsWatts model of dielectric relaxation. Jour. Stat. Phys., 34:129, 1984.
 [47] C. Reidl and B.N. Miller. Population dependence of early relaxation Phys. Rev., E51:884888, 1995.