From chaos to cosmology: insights gained from 1D gravity

The gravitational force controls the evolution of the Universe on several scales. It is responsible for the formation of galaxies from the primordial matter distribution and the formation of planets from solar nebulae. Because the gravitational force is singular and has infinite range, making predictions based on fully three-dimensional models may be challenging. One-dimensional (1D) Newtonian gravity models were proposed as toy models for understanding the dynamics of gravitational systems. They can be integrated exactly and were used for computer simulations starting in the 1960s, providing the first demonstration of violent relaxation and the rapid development of long-lived quasi-stationary states (QSS). The present review provides the bases of the physics of 1D gravitational systems. It is divided into two main parts, the first concerning the approach to equilibrium and the second applications to cosmology. Each part is self-contained and can be read independently of the other. In the first part, we provide an introduction to the equilibrium thermodynamics of the one-dimensional gravitational sheet (OGS) system in the Vlasov limit. Both fixed and periodic boundary conditions are considered. The relaxation to equilibrium of the OGS is studied through numerical simulations which establish the role played by QSS and violent relaxation. We also survey existing work on the Lyapunov exponents of the OGS and on the chaotic dynamics of 1D systems with few particles, focusing on the 1D three-body problem. The second part summarizes work on dynamical structure formation in cosmology using 1D systems. By transforming to comoving coordinates, which follow the global expansion of the Universe, the 1D approach provides a useful laboratory for studying structure formation in various cosmological scenarios, from Einstein-de Sitter and ΛCDM to more recent, alternative cosmological models. A key result is the appearance of scale-free behavior with fractal dimension, which can be reliably studied in 1D for large systems over many epochs. Finally, an appendix gives some details on the numerical simulation methods used in these studies.

equilibrium of the OGS is studied through numerical simulations which establish the role played by QSS and violent relaxation. We also survey existing work on the Lyapunov exponents of the OGS and on the chaotic dynamics of 1D systems with few particles, focusing on the 1D three-body problem. The second part summarizes work on dynamical structure formation in cosmology using 1D systems. By transforming to comoving coordinates, which follow the global expansion of the Universe, the 1D approach provides a useful laboratory for studying structure formation in various cosmological scenarios, from Einstein-de Sitter and ΛCDM to more recent, alternative cosmological models.
A key result is the appearance of scale-free behavior with fractal dimension, which can be reliably studied in 1D for large systems over many epochs. Finally, an appendix gives some details on the numerical simulation methods used in these studies.
Keywords: one-dimensional, gravity, cosmology, dynamical systems, fractals, Vlasov, computer simulations (Some figures may appear in colour only in the online journal)

Introduction
Self-gravitating systems of particles appear in many problems of astrophysics and cosmology. Systems as varied as stars, star clusters and galaxies are essentially N-particle systems with masses m i interacting by the two-body attractive Newtonian potential energy In the language of classical mechanics, the N-body problem is expressed as the solution of a system of coupled ordinary differential equations with G Newton's gravitational constant, and initial condition {⃗ r i (0),⃗ v i (0)} where the velocity is⃗ v i (t) =⃗ r i ′ (t). Here, the prime denotes differentiation with respect to the time t. The existence and uniqueness of the solutions {⃗ r i (t)} to these equations is ensured for −t 2 < t < t 1 where t 1,2 are associated with times when r(t) → 0, with r(t) = min i,j |⃗ r i (t) −⃗ r j (t)|. See chapter 2.1 in [1] for a detailed discussion. Except where otherwise noted, in the following we assume that all masses are equal.
The solutions of the Newtonian N-body problem for N > 2 have a complex behavior. In the absence of a total angular momentum, total collapse can occur, where all particles come together at the same time. This result is known as the Sundman theorem: total collapse cannot occur unless the total angular momentum is zero. See chapter 2.4 in [1] for a discussion. Chaotic behavior is also present, for example in the Solar System [2] and the galaxy [3].
Ejection to infinity in finite time can also occur for N ⩾ 5. See [4] for a detailed explanation and a historical perspective. Collapse and ejection to infinity are due to the properties of the gravitational potential in three dimensions V(r ij ) = −Gm 2 |⃗ r i −⃗ r j | −1 . This is singular as r ij → 0 and it vanishes as r ij → ∞, with r ij = |⃗ r i −⃗ r j |. The singular behavior at r ij → 0 can be remedied by adding a hard core to the interaction, which prevents the particles from approaching each other too closely, or a soft core, which modifies the potential 1/r ij at small distances.
A simpler approach is to study gravitational systems in 1D or 2D. The gravitational potential is determined by the solution of the Poisson equation in Euclidean space with the requisite dimensions. Thus, for these cases, the gravitational interaction potential has the form 1D : U(r) = 4πGm r (2) 2D : U(r) = 2πGm log r.
The dynamics of highly symmetric, higher dimensional, gravitational systems, such as a system of concentric spheres, could in principle be considered as one-dimensional since each 'shell' is characterized by its radial coordinate (see for example [5]). However, the potential is not a solution of the 1D Poisson equation and those systems will not be considered in this review. One-dimensional treatments are natural for systems with planar symmetry, which are uniform in two directions orthogonal to an axis. For example, the density profile of a spiral galaxy (seen edge-wise) can be studied as a 1D model because to a good approximation the galaxy can be seen as being made up of parallel 'sheets' of stars, see figure 1.  real values, which corresponds to an unbounded gravitational system. The Hamiltonian of the system is This model was originally proposed by Oort [6] and Camm [7] to describe galactic density distributions along a line perpendicular to the galactic plane. In a parallel development, a modification replacing the attractive potentials with Coulomb potentials has been considered as the one-dimensional plasma model by Lenard [8] and Eldrige and Feix [9]. Historically the main use of the model has been to study the complex dynamics of relaxation to equilibrium of gravitational systems. 1D gravitational systems are easily simulated in dynamical simulations, and the approach to equilibrium can be studied precisely. These studies reveal a complex pattern of approach to equilibrium, as discussed in section 2.3.
The OGS model can be simulated dynamically exactly, since the equations of motion have a simple form between particle crossings. The main concern of the simulation is reduced to that of finding the next crossing, which can be done in closed form. Consequently, individual particle trajectories can be followed up to the computer precision and the simulations can deal with a large number of particles. Moreover, in the case of neutral systems (see sections 2.5 and 3), simulations can deal with a large slice L of the system compared to the unit of length provided by the Jeans length [3,10]. 1D simulations are able to follow systems with both a large number of particles per Jeans length and a large number of Jeans lengths. Thus, in Fourier space, functions can be computed over a large number of decades. This feature has led researchers to use it as a first attempt to understand complex behavior in self-gravitating systems. As a consequence N-Body simulations have been the principal tool for investigating the evolution of the 1D system. The simulation algorithms are summarized in the section B.1. In addition to N-Body simulations, evolution has also been studied in the Vlasov, or continuum, limit. Besides direct integration of the Vlasov equation in phase space, recently the approach to equilibrium was also investigated via a normal mode decomposition [11][12][13].
The statistical mechanics of self-gravitating systems requires special treatment, because the gravitational interaction falls off with the inter-particle distance too slowly to be integrable.
Recall that in three dimensions we have U(r) ∼ r −1 . Potentials of this type U(r) ∼ r −α with α < d with d the space dimension, are known as long-range potentials.
The thermodynamics of systems with long-range interactions have distinctive properties which are very different from those of the systems with short-range interactions: • Non-extensive thermodynamics [14] • Negative specific heat [15,16] • Slow relaxation to thermodynamical equilibrium, and the existence of quasi-stationary states (QSS) in which the system can be found for a long time.
The 3D gravitational potential has two features which lead to peculiar properties of gravitational systems. It is singular as r → 0 and it vanishes at infinity. The singularity leads to the Antonov singularity-the collapse of a gravitational system enclosed in a finite volume-and the vanishing at infinity leads to the gravothermal catastrophe, the expulsion of particles from a gravitational system. None of these issues appears in the 1D system where the potential is finite at zero separation and the potential is unbounded from above. For further discussion of these phenomena in d = 2, 3 dimensions see the review [17].
One-dimensional gravitating systems have principally been used to explore the approach to equilibrium and, more recently, to understand structure formation in cosmology. This review is divided into two parts: The first deals with equilibrium thermodynamics and statistical mechanics and the approach to equilibrium, while the second part focuses on applications to cosmology. We will see that there is a large body of literature for each topic which we will address in the appropriate sections.
In section 2.2 we introduce the Vlasov limit and discuss in some detail the solutions of the static Vlasov equation for a gas of particles moving along a line and interacting with 1D gravitational potentials. We cover both cases of a continuous system, and the 1D lattice gas of particles interacting by 1D gravitational potentials. The exact statistical mechanics of the 1D gravitational gas with N particles with free boundaries, was derived by Rybicki for both the canonical and microcanonical ensembles [18]. In the infinite volume limit of the continuous system we recover the results of Rybicki by an alternative approach.
In section 2.3 we discuss the approach to thermal equilibrium of 1D gravitational systems. Thermalization is a complex process, and numerical studies show that it proceeds in two stages, with the formation of so-called QSS that are very long lived. Understanding the long lifetime of QSS has been an elusive goal. A natural formalism for the description of these states is the Lynden-Bell statistics, which was introduced in astrophysics for modeling the phenomenon of violent relaxation. Using a dynamical systems approach, thermalization is related to a complete mixing in the phase space. We summarize the status of knowledge about the spectrum of the Lyapunov exponents of the OGS.
In order to understand the ergodic properties of the 1D system that result in long lived quasi stationary states, investigators turned to systems with a small number of particles. Section 2.4 discusses the dynamics of 1D gravitational systems with a few particles N, paying special attention to the cases N = 3, 4. Numerical studies show that the dynamics is mostly periodic and quasi-periodic, and the regions of periodicity are surrounded by regions of localized chaos. Numerical simulations suggest that as the number of particles increases, the size of the stable regions surrounding periodic orbits in the phase space becomes vanishingly small.
Recently, motivated by applications to cosmology, the statistical mechanics of the 1D system with periodic boundary conditions was also investigated in some detail [19]. Surprisingly a phase transition does occur in this case, in contrast to the unbounded system [18] where no phase transition occurs. This work is presented in section 2.5 which discusses the properties of 1D systems with periodic boundary conditions, motivated by applications to cosmological structure formation.
In the second part of the paper we consider applications of 1D models to the expanding Universe. Both observation and theory tell us that, on large scales, the Universe is homogeneous and isotropic and is undergoing a uniform expansion with a time dependent expansion factor a(t). Cosmological simulations are an essential tool for understanding the evolution of structure formation following the epoch of recombination. They are carried out in the comoving reference frame expanding with the Universe and characterized by constant density. Because 3D simulations are compromised by requiring a short range softening length and some type of long range cutoff, they are less than ideal for studying the self-similar or fractal aspects of cluster formation that appear to underlay structure formation in the Universe. Conversely, 1D models that are susceptible to precise evolution algorithms are perfectly adaptable to fractal analysis (see appendix B).
Rescaling the spatial coordinate according to the expansion and also the unit of time leads to generalized comoving coordinates. The existence of a stationary state provides an alternative derivation of the Friedmann equation [20,21]. Section 3.2 shows how different models can be formulated depending on its parameters. It was shown by Aurell and Fanelli that the 1D system could be regarded as a perturbation embedded in a 3D expanding model [22]. In section 3.3 we will show how this approach can be generalized to an arbitrary number, d, of dimensions and that a 'conservative' or frictionless system appears in the limit of large d. The first 1D applications considered cluster development in the Einstein-de Sitter cosmology that represents a matter dominated Universe. This application is considered in section 3.3 and is accompanied by a new set of high precision simulations including the real space and phase space particle distributions, two particle correlation function and density spectrum. The large scaling ranges that develop over time are useful for estimating fractal dimensions, and this is included in section 3.3. Since dissipative visible matter is a small part of the Universe, in section 3.4 we consider a two component model that includes both a dissipative and conservative component representing respectively Baryonic and the more prevalent dark-matter.
A difficulty of N-body simulations is that they do not provide a good representation of the sparse matter distribution in voids. To address this situation, 1D investigations were extended to include a continuum, collissionless, or Vlasov representation of the system in the comoving frame. As discussed in section 3.5, for the Einstein-de Sitter case it was shown that the fractal description is improved at low density. In particular the generalized fractal dimension exhibits the correct monotonicity requirement.
Three challenges for the standard ΛCDM cosmology are the absence of theoretical justifications for the assumptions of inflation, dark matter and dark energy. The Dirac-Milne cosmology provides an untested alternative where antimatter is prevalent and repelled by ordinary matter [23]. In section 3.6 we explore the evolution of a particular 1D version of the DM cosmology and compare it with the standard ΛCDM version. We will see that structure formation 'freezes out' after the current epoch in each scenario.

Generalities
The natural description of the OGS in the limit of a very large number of sheets N → ∞ is taken at fixed mass M = Nm. This approach leads to the Vlasov equation for the dynamics of the phase space density of the system.
Of particular importance are the solutions of the static Vlasov equation, which comprise, but are not limited to, the state of the OGS in thermodynamical equilibrium. The thermodynamical equilibrium solution is recovered when the momentum distribution is the Maxwell-Boltzmann distribution. The density of the gas in thermal equilibrium satisfies the system of equations, see [24]: where C is a normalization constant and the potential ϕ(x) is expressed in terms of the twobody interaction potential U(x − y) as an integral over the volume V occupied the system In this section, we present a detailed solution for the gas density and the thermodynamical properties of the OGS in a finite volume. Taking the volume to infinity recovers the well-known Rybicki solution [18]. In addition to the continuous OGS system, we discuss also the statistical mechanics of the one-dimensional lattice gas of particles interacting by linear attractive potentials.
The study of the approach to thermal equilibrium is one of the main motivations for the study of 1D gravitational systems, as toy models for understanding the similar relaxation process in three dimensions. The main advantage of the 1D systems is that their simulation can be performed exactly, due to the simplicity of the equations of motion for the individual particles. In particular, it is possible to compute exactly the particle coordinates between collisions.
Numerical studies of the relaxation in 1D gravitational systems showed a number of surprising properties. First, the relaxation process is slower than expected, although the precise scaling with the number of sheets N is still debated. Second, the relaxation proceeds through the formation of long-lived intermediate QSS, followed by eventual relaxation to the Boltzmann equilibrium states. Essentially they have been used to investigate the relaxation towards thermal equilibrium and its time scale. They revealed the existence of quasi-stationary states (QSS) in which the OGS appears to remain for a long time, and the violent relaxation mechanism proposed by Lynden Bell to account for the regularity of the light emitted by elliptical galaxies [25]. In Section 2.3 we summarize the different results obtained from numerical simulations, followed by a summary of violent relaxation theory, as well as an illustration and theoretical ideas regarding the QSS.
The nature of the QSS has been studied both numerically and theoretically. Their properties depend on the initial state of the gas. For a particular set of initial conditions satisfying the virial condition, the QSS coincide with the solution of the Lynden-Bell statistics, proposed in 1967 to explain the phenomenon of violent relaxation. For generic initial conditions, the QSS have a core-halo shape, which reflects the ejection of particles from the core, due to rapidly oscillating gravitational field acting on any given particle.
Another approach to the relaxation of the 1D systems uses dynamical systems theory. In this approach, the long-time evolution of the system is described by its Lyapunov exponents, which give the rate at which density fluctuations at different spatial scales grow or dissipate. We summarize the state of the art of the Lyapunov exponents for the N-particle OGS, and their scaling with N.
The study of the dynamics and ergodic properties of OGS with a small number of particles have provided insight into understanding the relaxation process. Here, we focus on the three body problem N = 3 which is the best studied case, and give a detailed discussion illustrating the evidence for chaotic behavior.

1D gravitational systems in thermal equilibrium
Although a numerical treatment of the ODE's (1) is possible, the complexity grows rapidly with N. A simpler analytical treatment becomes possible in the limit of a large number of particles N → ∞, m → 0. In this limit the system can be approximated as a continuous fluid, with density f(⃗ r,⃗ v, t) in the phase space, where⃗ r,⃗ v are the position and velocities of the particles in a small volume.
2.2.1. The Vlasov limit for the gravitational gas. Systems with short-range interactions simplify in the large N limit in the so-called thermodynamical limit. In this limit both the number of particles and the volume are taken very large, at fixed particle number density ρ = N/V. In this limit the surface effects vanish and the thermodynamical properties are dominated by the bulk contribution. A rigorous treatment of these systems is given by Ruelle in [26]. A wide class of systems of this type includes systems with hard cores and Kac type interactions, and were shown by Lebowitz and Penrose to have a well-defined thermodynamical limit [27].
For systems with long-range interactions and especially for gravitational systems, a more appropriate limit is the so-called Vlasov limit. Assuming that the system is enclosed into a fixed volume V, this limit is defined by where T V defines the so-called Vlasov temperature. In this limit the distribution function f satisfies the Vlasov equation where ϕ(⃗ x, t) is the mean-field potential given by the solution of the Poisson-Vlasov equations The coupled system of equations (9) and (10) must be solved simultaneously. We consider next the mean-field limit of a system consisting of a large number of interacting particles N in a finite volume V. This will be used to justify the Vlasov limit introduced in (7). This limit can be shown to correspond to the time-independent limit of the Vlasov equation, and is known in the literature as the static Vlasov limit [24].
Consider a system of N particles enclosed in a domain Ω of volume V, and interacting with the Hamiltonian where U(r ij ) are two-body potentials. The strength of the two-body interaction becomes weak in the limit of a large number of particles N → ∞. Assume that the system is maintained at fixed temperature T, which corresponds to the canonical ensemble. It can be shown [3,24,28] that for a wide range of potentials U(r), the free energy per particle f = F/N approaches a finite value as N → ∞ and furthermore this limit is given by the solution of the variational problem where ρ(r) is constrained as´Ω ρ(r)dr = 1. The two terms correspond to the interaction energy and the entropy per particle system, respectively, f = u − Ts. The solution of the variational problem ρ(r) is the one-particle distribution function and gives the probability that a particle is found in the volume element d⃗ r centered on ⃗ r. For a collection of solvable problems with various choices of interaction potentials in 1D, 2D and 3D, see [24]. This reference gives also precise conditions on the potentials U(r) for which the limit (12) exists.
The Hamiltonian of a system of N particles of mass m interacting by gravitational potentials is written as The canonical partition function is Z N = α e −βHgrav where the summation extends over all states of the system. We make several observations: • In the interaction term, replace one factor of m with M/N according to the scaling (7).
• The resulting Hamiltonian is linear in m, the particle mass. When substituted into the partition function, the factor of m can be absorbed into a redefinition of the temperature β V = mβ, (7). • Under the Vlasov scaling, the interaction Hamiltonian scales like O(1/N), and has the typical form of the mean field interaction (11).
Under the rescaling mN = M and temperature redefinition β V = mβ, we have βH grav = β V H grav with which admits a mean-field limit as N → ∞. The single particle density function ρ(x) is found by minimizing the free energy F = U − TS written as in (12). This shows that a self-gravitating system of N particles admits a mean-field limit as N → ∞. This is the static Vlasov limit. In lower dimensions the mathematical treatment of the Vlasov equation is greatly simplified and analytical results can be obtained in several cases. In the following paragraphs, we discuss the solution of the static Vlasov equation for the simplest case of the 1D gas of particles interacting by 1D gravitational forces.
Consider a one-dimensional gas of N particles enclosed in a finite volume [−L/2, L/2], and interacting by two-body attractive potentials The equilibrium thermodynamics of this system has been studied first by Salzberg [29], who considered a gas of N particles in a box of volume L, and interacting with potentials (15) with hard cores d (the particles are not allowed to come closer than d). The mass m is kept fixed. This leads to non-extensive thermodynamics: the total gas energy scales as U ∼ N 3 as N → ∞. Working in the isothermal-isobaric ensemble, the equation of state has been also derived, and has the form L = Nd + 2k B T/p. This is essentially the free gas equation of state corrected by the hard core p = 2k B T/(L − Nd). These properties are not realistic, and this is seen to be due to the large range properties of the 1D gravitational interaction, which require a more careful treatment of the thermodynamical limit.
The solution of the 1D gravitational gas in the Vlasov limit in an unbounded volume has been obtained by Rybicki [18], both in the microcanonical and the canonical ensembles. This paper considered a gas of N particles interacting by potentials (15) (without hard cores), which can move on the infinite line. The solution of the one-particle distribution function was found for arbitrary finite N, and the N → ∞ limit was considered at fixed total mass M = Nm and total energy E. This corresponds to scaling the particle masses as m = M/N. We present an alternative derivation of the result in [18] using the Vlasov limit, working in a more general setting of a finite volume.
As explained in the previous section, one factor of the mass m is absorbed into a modified temperature mβ → β V . Furthermore, scaling m such that the total mass is fixed M = mN, the interaction Hamiltonian is put into the mean-field form The particle positions are bounded to the interior of the box x i ∈ − 1 2 L, 1 2 L . The properties of this system in thermodynamical equilibrium in the canonical ensemble in the limit N → ∞ are given by the solution of the variational problem for the configurational free energy per particle f Q with g 2 = 2πGM, and the single particle density function ρ(x) satisfies the normalization con-dition´L The solution of the variational problem (17) for the single particle density function is given by the following result. Proposition 1. The single particle density distribution function of the gas of particles interacting by 1D gravitational interactions in thermodynamical equilibrium at temperature T is where δ is the solution of the equation (19) and the central density is The gas density profile has the form ∼ 1 cosh 2 ( 1 2 δx) . This was derived first by Camm [7], and is known as the Camm law for the galactic profile density. It has been shown to be well satisfied by the density profile of the galaxies along the direction orthogonal to the galactic plane.
In the infinite volume limit L → ∞, the parameter δ approaches a finite value which depends only on temperature This parameter determines the size of the gas cloud as ∆x ∼ 4 δ = 4 βg 2 and the central density is ρ(0) = 1 4 βg 2 . As the temperature increases, the size of the gas cloud also increases, and becomes infinitely diffuse in the infinite temperature limit. Conversely, as the temperature decreases the gas becomes more and more concentrated until it collapses to a single mass point in the zero temperature limit. We note that the infinite volume result for the central density ρ(0) = 1 4 βg 2 agrees with equation (2.37) in Rybicki [18] for the central value of the single particle distribution function. In our notations this is See the appendix A for the proof of proposition 1. The thermodynamical properties of the gas in the Vlasov limit are given by the following result. Proposition 2. The free energy per particle f = F/N = f Q + f kin of the gas with interaction V(x, y) = g 2 N |x − y| enclosed in a volume L in thermodynamical equilibrium has a configurational contribution and a contribution from the kinetic degrees of freedom In (22) δ is the solution of the equation (19). The energy u = U/N and entropy s = S/N per particle are given by u = f − T∂ T f and s = −∂ T f. Explicitly The first term is the contribution from the interaction energy u int , and the second term is the kinetic energy u kin . The plot of the interaction energy per particle u int vs 1 2 δL is shown in figure 2.
See appendix A for the proof of proposition 2.
The equation of state of the gas is The first term is the ideal gas equation of state, and the second term is a negative correction to the pressure, due to the linear attractive potential. The plot of the function Ψ(x) is shown in figure 2 (right panel). This has the small-x expansion There is a non-zero correction to the ideal gas equation of state even in the infinite temperature limit! This correction term decreases as the temperature is lowered. We note from (19) that this equation has solutions for δ for any values of L, T. The equation has two solutions ±δ with either sign. As seen from (22), the gas properties depend only on δ 2 , so the sign ambiguity is immaterial. The configurational free energy per particle has the form f Q = k B T(− log L + Φ( 1 2 βL)) depends on volume L through the combination δL which is a function of βL, as seen from (19). This means that the free energy is linear in N, the particle number, but not in L, the gas volume. This is a consequence of the long-range nature of the interaction.
From figure 2 (left) we see that lim x→∞ f u (x) = 1 which means that the interaction energy per particle u int approaches k B T for L → ∞. This is the infinite volume limit, considered by Rybiki [18]. In this limit the total energy per particle is which agrees with the result in equation (2.40) of Rybicki [18]. This result also follows from the virial theorem for the gas with linear interactions, which says that the kinetic energy and interaction energy are related as 2U kin = U int . Using u kin = 1 2 k B T reproduces the result (28) in the infinite volume limit.
We outline the derivation of the equation of state (25). The pressure of the gas is given by Recall the result for the configurational free energy per particle from equation (22), which can be written as Taking a derivative of f Q with respect to L gives The second line in (32) was obtained as follows: recall that δ depends on L through the equation (19). Taking a derivative of this equation with respect to L we get Substituting into the first line of (32) gives the second line. Evaluating the derivative Φ ′ (x) in (32) reproduces the final result (25) for the gas pressure. This completes the proof of the equation of state.
Let us consider the behavior of the 1D gravitational gas in the large temperature limit. We take the β → 0 limit of the exact solution for the thermodynamical parameters of the gas. An examination of (19) shows that in this limit we have δ → 0 as Substituting into (24) gives A similar result is obtained for the free energy (22) which can be written as with Φ(x) defined in (31). Using the small-x expansion Φ(x) = 1 The equation of state (25) is approximated, in the large temperature limit, as where we used the small-x approximation Ψ(x) = 1 3 x + O(x 3 ). Note the finite correction − 1 6 g 2 N to the ideal gas equation of state which survives even in the infinite volume limit. The correction term in (39) coincides with the interaction energy of a uniform gas with density ρ(x) = 1 L . In the large temperature limit we expect the gas density to approach a uniform density. The interaction energy per particle becomes The presence of this correction term is a feature of the long-ranged nature of the onedimensional gravitational interaction, which introduces a volume dependence of the interaction energy. This is linear in L, with a coefficient which is precisely the correction term in (39).
It can be shown that the gas pressure given by (25) is always positive p > 0. This is not obvious since it is the difference of two positive terms. This can be seen by noting first the inequality This follows from the inequality 1 2 δL ⩾ 1 2 βg 2 L which is obtained from the equation (19) for δ, and noting that Ψ(x) is a decreasing function.
Substituting the inequality (41) into the equation of state (25) this gives which follows from the inequality xΨ(x) ⩽ 1.
The stability of the gas against collapse requires that the isothermal compressibility is non-negative Taking the derivative of the pressure (25) we get The last factor is evaluated using (34), which gives and h(x) is defined in (33). The function Ξ(x) is defined as The plot of this function is shown in figure 3. It is negative, has the small x Taylor expansion and it approaches 0 as x → ∞. It reaches a minimum value of Ξ(x 0 ) = −0.032 at x 0 = 3.5.
Using the lower bound on Ξ(x), the compressibility is bounded from below as The multiplier is positive for 1 2 βg 2 L ⩽ 5.59. It is also positive for 1 2 βg 2 L ⩾ x 0 = 3.5, since in this region we have Ξ( 1 2 δL) > Ξ( 1 2 βg 2 L) from the increasing nature of Ξ(x) in this region, and the inequality 1 2 δL ⩾ 1 2 βg 2 L following from (19). We have further which follows from x 2 Ξ(x) ⩾ −1 which is obtained by numerical study of the function Ξ(x). The two regions 1 2 βg 2 L ⩽ 5.59 and 1 2 βg 2 L ⩾ 3.5 cover the entire real positive axis, such that we conclude that the isothermal compressibility κ T is always non-negative.
This implies that there is no Antonov instability for the 1D gravitational gas. A similar conclusion was noted in [30]. It is known that in three dimensions an isothermal gas with total energy E and mass M enclosed in a spherical container of radius L is unstable against gravitational collapse provided that L > 0.335GM 2 (−E) . This is the Antonov instability [14,31]. This phenomenon does not occur in the 1D gravitational gas.

1D gravitational lattice gas.
To our knowledge, there is no treatment in the literature of the 1D gravitational gas with hard cores in the Vlasov limit. Since the interaction potential is not singular as |x − y| → 0, it may seem that the introduction of a hard core will not change qualitatively the thermodynamical properties. This is intuitively clear at large temperature, where the distances between particles are large, and they do not come within distances of the order of the hard core. However, at small temperature we expect deviations from the no-hard core behavior because the hard cores limit the density of particles to a finite value.
A simple model which simulates the presence of hard cores in the interaction is a lattice gas. The gas particles occupy the sites of a uniformly spaced lattice, and at each site at most one particle can be present. The configurational partition function is a sum over all possible placements of the gas particles on the lattice sites.
Assume that the lattice has n sites, and the gas contains N particles. This system can be mapped to a system of n − N non-interacting bosons which can be placed on a set of N + 1 energy levels. The thermodynamical properties of the system have a striking resemblance to the Bose-Einstein statistics, which is reflected also in the thermodynamical properties of the system [32,33].
The equivalence to a bosonic system can be seen by denoting the positions of the N particles on the lattice by i 1 < i 2 < . . . < i N where i j ∈ {1, 2, . . . , n} are integers. Denote y j = i j+1 − i j − 1 the number of empty lattice sites between the consecutive particles {j, j + 1}. (The right-most particle requires special treatment, we define y N = n − i N .) The energy of the system can be written as where the ground state energy corresponds to the particles occupying a block of consecutive N lattice sites and the energy of the excited state is a sum of integer multiples of The variables y k satisfy the sum rule N k=1 y k = n − N which has a simple geometrical interpretation as the total number of empty lattice sites. We see that the {y k } N k=1 can be interpreted as the occupation numbers of the n − N energy levels ω, each configuration being counted only once.
Assuming the scaling g 2 = cn −2 , the thermodynamical properties of this system approach a well-defined limit as n, N → ∞ at fixed particle number density d = N/n. This limit corresponds to taking mn = constant as n → ∞, which is identical to the Vlasov limit mN = M defined above in (7).
The free energy and equation of state have been derived in [32] under this scaling, in the isobaric-isothermal ensemble. The thermodynamical properties of the lattice gas with linear attractive potentials are given by the following result [32] 5 Proposition 3. The free energy density f = F/n of the lattice gas with interaction ε ij = g 2 n 2 |i − j| and particle number density d = N/n is where π(d, T) is the solution of the equation and J(d, T) denotes the integral The equation of state is The lattice gas with linear attractive potentials and a constant external field was studied recently in [33]. This corresponds to the two-body interaction ε ij = 1 n 2 |i − j| + x n 2 (i + j) where x is the external field strength. A novel feature for this system is the appearance of a gas-liquid phase transition, for arbitrarily small field strength x.
The thermodynamical properties of the lattice gas with two-body linear attractive interaction potential ε ij = g 2 |i − j| can be solved also exactly for finite n, N, using a recursion relation for the canonical partition function in the lattice size n. Denoting Z n (N) the canonical partition function of a lattice gas with n sites and N particles interacting with potentials g 2 |i − j|, we can show [32] that it satisfies the recursion relation 5 Reference [32] considered a more general form for the 2-body interaction including also an universal attractive term εij = − 1 µn ξ + 1 µn 2 |i − j|. The pure gravitational case corresponds to ξ = 0. We denoted 1 µ = g 2 .
This follows by expressing the total interaction energy of the N particles as Using this equation one can express the energy of the lattice with n + 1 sites and N particles in terms of the energy of the sub-lattice with n sites as: (i) E n+1 (N) = E n (N) if no particle is present at site n + 1, or (ii) E n+1 (N) = E n (N) + g 2 (N − 1)(n + 1) if a particle is present at site n + 1. This gives the recursion relation (58) for the canonical partition function Z n (N). The recursion (58) can be easily solved numerically starting with the initial condition at n = N: The numerical results for the thermodynamical properties (equation of state) following from the recursion approach were found to agree well with the exact solution for lattices with n ∼ 100-200 sites [32].

Thermalization and quasi-equilibrium states
Experience with systems with short range interactions suggests that gravitational systems should approach thermal equilibrium in the infinite time limit. The situation is different for integrable systems, or close to integrable, such as the Fermi, Ulam, Pasta system [34], which show a pattern of energy transfer between different degrees of freedom, which slows down the approach to equilibrium. Theoretical arguments and numerical experiments with the OGS show a more complex pattern of the phenomena in the relaxation to equilibrium. Numerical simulations have been extensively used to explore the approach to equilibrium of an OGS and we review them here.
It was recognized early on that the relaxation to equilibrium in the OGS is slower than in the typical 3D gravitational systems. For the latter case, Chandrasekhar [35] showed that the relaxation time scales like ∼ N log N with the number of particles N. However some of the approximations used to obtain this scaling in 3D do not apply to the 1D system. Early studies by Hohl, Feix and collaborators [36][37][38] proposed that the OGS relaxes to thermodynamical equilibrium on a time scale of the order of N 2 t c , where N is the number of sheets and t c is the typical time it takes one sheet to cross the system. The precise definition of t c is author dependent and t c ∼ 4πG/ρ, where the density ρ is a mean density taken over a typical length. By studying the correlation function, Hohl and Broaddus [37] observed relaxation to thermal equilibrium in numerical experiments also on a time scale of t R ∼ N 2 t c .
This scaling was challenged as high performance numerical simulations became more accessible. Wright, Miller and Stein (WMS82) [39] studied the approach to equilibrium for a OGS with N = 100 particles, started in three initial states: (1) uniform rectangular distribution in phase space with a high virial ratio V R = 2⟨T⟩/⟨U⟩ ≫ 1; (2) uniform rectangular distribution with virial ratio V R = 1.0 (corresponding to the thermal equilibrium value) and (3) the isothermal distribution (as a control). For case (1), after running the simulation for ∼20 000t c , they concluded that the system did not approach thermal equilibrium, but instead transitioned rapidly into a metastable state with a different configuration characterized by a halo surrounding a cold dense core and remained there. They employed the Kolmogorov/Smirnov goodness of fit test to show with high probability that neither the position nor velocity distributions were sampling the equilibrium distributions. In case (2), after a run of 1000t c , the system retained a uniform distribution in phase space) while for case (3), considered as a control, the system remained in isothermal equilibrium.
These results were challenged in turn by Severne et al [40] who found evidence of a rapid approach to equilibrium, also as measured by the Kolmogorov-Smirnov statistic for the position and momentum distributions, for a different subset of initial conditions. In response, a study that emphasized the dependence of the thermalization time on the initial distribution in phase space was presented by the WMS84 group in [41]. Later studies of the correlation functions of the OGS showed that equilibrium was not established on this time scale for different initial states [42].
In order to settle this apparent disagreement, Gouda, Konishi and Tsuchiya (GKT thereafter) revisited the issue of the relaxation process in the OGS in a series of three papers [43][44][45]. Similar to previous studies, they simulated the evolution of the system started in the water-bag distribution, which is a uniform distribution in the phase space.
In [43] GKT show that there exist two characteristic time scales: (i) the microscopic relaxation time T m ≃ Nt c , which describes the time required for energy to be distributed among all particles, although the resulting distribution is not the thermal equilibrium distribution, and (ii) the macroscopic relaxation time T M , which is much longer and gives the time required for approach to thermal equilibrium. During this time the system is in a so-called quasi-stationary state QSS (see the next section for a more detailed discussion). These time scales were further investigated in [44], which revealed a more complex pattern, showing that the system switches between quasi-equilibrium states. This itinerant behavior was studied in [45].
To summarize their work, we can distinguish three main regimes in the relaxation process of the OGS. The first regime 0 < t < T m is the so-called collisionless phase, the initial particle distribution is slowly relaxed, on the time scale of the microscopic relaxation time. In the second regime T m < t < T M , the system stays in a long-lived quasi-stationary state. Finally, the socalled macroscopic relaxation time denotes the time scale for relaxation to thermodynamical equilibrium where the distribution approaches the microcanonical distribution. Nevertheless, up to five steps have been detected [45] showing that the situation is still not very clear.
A related system consists of two distinct species of particles with different masses. This bimodal mass distribution was first suggested by Hohl and Tilghman Broaddus [37]. It was employed to study the Lynden-Bell hypothesis and the development of the quasi-stationary state [38]. It was also used to study the approach to equilibrium [40]. The ratio of kinetic energy of each each mass species was employed as a progress variable to quantify the deviation from the fully relaxed equilibrium state. Equipartition can be taken as an additional indicator of thermalization. These simulations suggested a tendency towards equipartition but it was never achieved in the time scale considered. Later, with improved algorithm design and hardware, Yawn and Miller observed complete equipartition in a bimodal mass system over a very long time scale [46]. They also derived equations for the equilibrium properties of each species and showed that fluctuations are very long lived and correlations decay extremely slowly [47].
The persistence of structures in phase space (in fact holes) seems to prevent the system from reaching equilibrium. The influence and dynamics of these holes have been studied in Rouet and Feix [48] by means of the usual N-body code, while Mineau et al show that, in the Vlasov limit and except on average, the system does not relax to any function of the energy alone [49] (see section 2.3.4). An alternative method for studying the relaxation to equilibrium in the OGS was proposed in [50]. In this paper a stochastic diffusion model was proposed, and was used to study the exploration of the phase space by the system. The agreement of the model with the exact OGS system is studied by numerical simulation, and is found to be better at low energy.
Later studies by Joyce and Worrakitpoonpon [51] using larger systems with N ∼ 10 3 essentially confirmed this picture, i.e. that nearly virialized initial states relax towards equilibrium on a time scale T R ≃ Nt c . In contrast, using a different approach, Levin et al [52], found that the relaxation to equilibrium proceeds on a time scale t R ∼ N 1.8 using a limited set of initial conditions.
An important issue in studies of relaxation to equilibrium is the choice of the measure of distance to the equilibrium distribution. Natural choices are distributional measures such as the Kolmogorov-Smirnov goodness of fit test for the velocity and/or spatial distribution, taking into account the known equilibrium distributions (Boltzmann distribution for velocities, and the inverse cosh law for density). This was the measure used in [39,40]. A similar role is played by the ratio of kinetic energy between each species in a system with a bimodal mass distribution [40,46,47]. GKT employed the deviation of the time averaged single particle energies from their equilibrium values [43][44][45]. In [51] the authors employed the use of order parameters ϕ αβ defined by normalized moments of the phase space distribution. Levin et al [52] used a measure of the separation of the distribution in phase space between the QSS and final equilibrium.
A more recent result is due to Barnes and Ragan [12,13], who studied the relaxation of a one-dimensional self-gravitating system in the collisionless approximation. Working in the collisionless approximation has the advantage that the collisionless Boltzmann, or Vlasov, equation can be studied analytically in a perturbation theory approach. As a measure of the rate of relaxation they study the entropy growth rate. They identify a barrier to full relaxation in the form of time-independent modes, which imply that a separable equilibrium of the phasephase distribution cannot be reached either through phase mixing or by violent relaxation.
Motivated by cosmological considerations concerning cold dark matter, several authors have investigated the QSS formed from very cold initial conditions, specifically when all particles in the simulation start at rest. Most recently Tashiro has concluded that the density in the resulting QSS has a universal form that is a variant of that observed in two and three dimensions [53]. See also the related earlier work of Binney and Schulz et al [54,55] as well as Colombi and Touma [56]. Binney, in particular, explored the effect of discretization, or 'graininess', in the phase space due to the finite population of particles [54]. Graininess was also partially addressed by Romero and Ascasibar in related work by constructing ensemble averages from independently sampled initial conditions as a function of population [57], and also plays a role in the Schrodinger-Poisson formulation [58] where the effective Planck constant limits the possible structure produced by a pure Vlasov evolution.
We also mention briefly results obtained in the Hamiltonian mean-field model (HMF) which is a system of particles moving on a circle interacting by cosine potentials. This model is used as a toy model for studying the relaxation to thermal equilibrium which is more accessible to numerical simulation with larger numbers of particles. The study in [59] found a scaling for the relaxation time of the form t R ∼ N 1.7 .

Quasi-stationary State (QSS).
Numerical simulations of the dynamics of the OGS with initial conditions filling a finite volume in the phase space, demonstrated the appearance of QSS after the microscopic relaxation scale. These states do not have the spatial Boltzmann-Gibbs distribution expected for the thermal equilibrium states. The system spends a long time in these states, corresponding to the second stage of the relaxation noted in the previous section. Numerical studies by Tsuchyia, Gouda and collaborators [43][44][45] demonstrated the important role of the long-lived QSS for an understanding of the approach to equilibrium in the OGS.
We illustrate the pattern formation in the phase space in the relaxation process with a numerical simulation for a system of N ∼ 32 000 particles of equal mass initially distributed uniformly in a square in phase space (position-velocity). The initial state has virial ratio V R = 2.0. Figure 4 shows the time evolution of the phase system, showing separately the spatial density (left), velocity distribution (middle), and the phase space distribution (right). The initial distribution is uniform within a square in phase space. There is no physical difference between particles. They have just been labeled depending on their energy at T = 70, then have been color coded according to their label at each time.
The main feature of the simulation, the formation of a core-halo structure, can be observed here. It is clear that the particles that belong to the halo initially have the highest energy. Because of the differential rotation (particles of high energy rotate more rapidly than particles of low energy), the halo exhibits filamentation due to the higher rate of stretching. As time goes on, the finite number of particles is too small to keep following this phenomenon precisely. In the infinite particle (Vlasov) limit the filaments maintain their identity indefinitely. Because of the area preserving property of the phase space evolution, empty regions are trapped between the filaments. With a finite number of particles the larger trapped regions are apparent and have the appearance of holes or cavities in the distribution. Thus two rotating holes can be observed in this region. The core alone seems to quickly reach a stationary state characterized by the constant density in phase space retained from the initial condition. We emphasize that neither the halo, which has not yet relaxed here, nor the constant density core, represent the final, thermal, equilibrium derived by Rybicki.

Lynden-Bell statistics and violent relaxation.
The formation of the QSS has been related to the process of violent relaxation proposed by Lynden-Bell [25]. In the collisionless limit, the 2-body interactions between particles can be neglected, and the dominant interaction is with the mean field produced by the remaining N − 1 particles. Under certain conditions the relaxation proceeds through a rapid change of the mean field, and this has the effect of ejecting particles which form a halo surrounding a central core.
We give in this section a brief summary of the Lynden-Bell statistics for the collisionless dynamics of the Vlasov equation. This will be used in the next section to analyze the QSS which are the equilibrium states of this statistics, following [60]. The collisionless Boltzmann equation, or Vlasov equation has the form where Φ(x) is the gravitational potential. For 1D gravity this is with g = 2πGm.
In the collisionless limit the phase space density f is constant along trajectories in the phase space. This implies that non-overlapping phase space regions remain so at all later times. This is equivalent to a microscopic exclusion principle for the phase space elements. Let us divide the phase space into a large number of microcells of volume ω, and denote the phase space density η. The mass associated with each microcell is ηω. The state of the system can be specified the set of occupied microcells ω j and their densities η j .
The Lynden-Bell statistics describes the state of the system in terms of the coarse grained phase space distributionf, defined by combining the phase space densities η j of the microcells. The equilibrium coarse grained distribution is obtained by maximizing the entropy under the constraints of constant total energy and mass.
The result for the equilibrium distribution was obtained by Lynden-Bell and has the form (see appendix I in [25]) with µ j , β j Lagrange multipliers and j is an index summing over all microcell in a coarse grained space phase element.
In the non-degenerate limit, given byf j ≪ η j , the coarse grained phase space distribution turns into a sum of Maxwellian distributions with different velocity dispersions This corresponds to a superposition of thermalized Gibbs-Boltzmann states.
The relevance of the Lynden-Bell statistics for this problem was noted already in the early work by Hohl and Campbell [38] and Cuperman et al [61,62]. Along with Harten, they also studied the Lynden-Bell hypothesis in the collisionless limit, integrating the boundaries of water-bag distributions [63,64]. In the same period, Janin [65] use this approach to address the Jeans instability and Colombi and Touma [56] extended this method to investigate the form of the density for a cold initial state. Using N-body simulations, the Lynden-Bell hypothesis gained additional attention following studies by Yamashiro et al [66], Yamaguchi [67], Joyce and Worrakitpoonpon [51,60] and Teles et al [68], who performed extensive numerical simulations and demonstrated the connection to the QSS.
Joyce and Worrakitpoonpon [51] studied the relaxation of the OGS started in several initial configurations. For a range of simple water bag and cold initial conditions, numerical simulations show that the system evolves into two phases with very different time scales. The state evolves on a time scale characterized by the dynamical time t c (roughly defined as the crossing time of a particle through the system), to a QSS. Finally, on a much longer time scale scale, dependent on the number of particles, this approaches thermal equilibrium. They found the characteristic time scale for relaxation behaves, to a good approximation, as where f QSS is a numerical factor which depends on the initial condition, and N is the number of particles, but others don't necessarily agree (see [52] and below). The detailed shape of the QSS depends on the initial condition. Recall that in thermal equilibrium the average kinetic K and potential energy of the OGS satisfy the 1D virial condition If the initial state satisfies the virial condition, the relaxation is non-violent, and the shape of the QSS agrees well with the prediction of the Lynden-Bell statistics. On the other hand, initial states which deviate from the virial condition have a violent relaxation, with large oscillations of the kinetic and potential energy, and eject particles into a halo surrounding a central core-like distribution.
These phenomena were noted in early numerical simulations by Goldstein et al [61] and Lecar and Cohen [69], and were confirmed in more recent simulations by Levin et al [52] and Joyce and Worrakitpoonpon [51].
In the particular limit of a cold initial state with vanishing velocities, the virial ratio V R = 2K U vanishes. The system evolves to a quasi-stationary state with a distinct shape, which was studied by Binney [54], by Shultz et al [55] and more recently by Tashiro [53]. The Lynden-Bell theory is based on the collisionless approximation. Going beyond this approximation requires that one takes into account dissipation effects, due to viscous damping or due to inelastic collisions. The impact of these effects was studied by Joyce et al [70], who showed that similar QSS exist, in appropriately rescaled variables.
Another issue concerns the robustness of the QSS under the impact of external perturbations, which become relevant for non-isolated systems. This problem was studied by Joyce et al [71]. They find that under a wide class of perturbations the system approaches a different QSS.

Lynden-Bell statistics predictions for QSS.
We study here the shapes of the QSS in the one-dimensional Lynden-Bell theory, following [51]. See also the detailed discussion of the predictions of the Lynden-Bell statistics for one-dimensional states in [69]. Denote the density ρ(x) and the gravitational potential φ(x). The energy is The coarse grained phase space density is The Lagrange multipliers β, µ are determined from the total mass and energy constraints The spatial distribution of the QSS is obtained from the Poisson equation We study in some detail the solutions of the Poisson equation for the Lynden-Bell statistics (69). Define the function where Li n (z) := ∞ k=1 z k k n is the polylogarithm. Expressed in terms of this function the Poisson equation has the form This equation must be solved with the boundary conditions φ(0) = φ ′ (0) = 0. For illustration purposes we set the normalization constants β, µ are set to their values from figure 21 of [60] β = 0.0035, µ = 3315 and assume for simplicity g = 1, f 0 = 1.
The contours shown correspond to the core (filled ellipse) with ε ⩽ ε F and the halo ε h .

Core-halo states.
According to the Jeans theorem [3], the distribution of the (quasi-)equilibrium states of the Vlasov equation can depend only on integrals of motion. The simplest integral of motion is the local energy functional ε = 1 2 v 2 + ψ(x). This led Teles et al [68] to postulate a phase space density of the form This consists of two uniform distributions: the core with uniform density η − χ which contains states with energy up to the Fermi energy ε F , and the halo, with density χ, with energies up to the halo energy ε h . See figure 6 for a graphical representation.
Integrating over velocity v, this gives the spatial distribution function ρ ch (ψ) :=´dvf ch (x, v) in the form The Fermi energy ε F and the mass of the halo χ are determined from the normalization conditions for mass and total energy (68).
For example, starting with a rectangular waterbag state with uniform density gxm from unity measures the deviation of the initial state from the virial condition (65).
Substituting ρ(x) from (73) into the Poisson equation (69), this can be solved as in the previous section, yielding the potential function ψ and the density profile ρ(x). The resulting distribution depends on the parameter V R : for V R very different from 1 the spatial distribution consists of a core similar to the LB QSS states, and a halo, extending to large distances. The halo contribution decreases as V R approaches 1, which is the limit when the initial state satisfies the virial condition. Under this condition the predictions of the LB theory are satisfied with a good approximation.

Chaotic dynamics
We recall briefly the definition of chaotic behavior. Chaotic dynamics are characterized by three properties: (i) the particle trajectories are bounded and steady-state, but distinct from periodic and quasi-periodic trajectories; (ii) if the system is dissipative they converge to a set in phase space called strange attractor, which is not a simple manifold like a point, circle or torus, but has a complex (fractal) geometrical structure with a fractional Hausdorff dimension, whereas if the system is conservative they are ergodic and eventually fill the available phase space; (iii) they exhibit sensitive dependence on initial conditions, such that chaotic trajectories locally diverge away from each other, and small changes in starting conditions build up exponentially fast into large changes.
A more rigorous study of the relaxation to equilibrium can be performed using methods from the theory of dynamical systems. In this approach, relaxation is related to mixing in the phase space, with a time scale τ KS = 1/h KS given by the Kolmogorov-Sinai entropy h KS . This entropy is also related to the Lyapunov exponents of the system, which give the characteristic time scales for the decay of perturbations. We discuss this approach in this section.
The dynamics of the OGS appears to be different for small numbers of sheets N, and for large N. For sufficiently small N the system is not ergodic, and a segmentation of the phase space is observed into trajectories of well defined type, with windows of chaotic behavior. On the other hand for large N > 11, evidence has been presented of ergodic behavior, although doubts still persist Yawn and Miller [72].
We start by summarizing the status of the knowledge of the Lyapunov spectrum and its implications for the relaxation dynamics of the OGS. Then we discuss in some detail the two cases with N below and above the critical value 10.

The Lyapunov spectrum of the OGS.
The Lyapunov exponents of a dynamical system give information about the time scale of the damping of fluctuations on different spatial scales. They are obtained by studying the exponential rates of divergence between two nearby trajectories. We give here a brief overview. For detailed accounts see for example Posch and Hoover [73], Benettin et al [74], Shimada and Nagashima [75].
Denote D the dimension of the phase space of a dynamical system. There are D Lyapunov exponents and the set {λ i } D i=1 is called the spectrum of the Lyapunov exponents. For ergodic systems the spectrum of the Lyapunov exponents is independent of the initial condition in the phase space. For systems with equations of motion which are simplectic and time reversible the Lyapunov spectrum is organized in pairs of zero sum. This implies that it is sufficient to compute only the positive exponents. The sum of all Lyapunov exponents vanishes as a consequence of Liouville's theorem for the conservation of the phase space volume for symplectic systems.
Numerical algorithms for their calculation were proposed by Benettin et al [74] and Shimada, Nagashima [75], and are based on periodic re-orthonormalization of the perturbation vectors using a Gram-Schmidt procedure.
The phase space of the OGS with N sheets has dimension D = 2N, such that there are They appear in pairs of zero sum λ j + λ N−j = 0. We assume that they are ordered as 0 < λ N−2 < λ N−3 < . . . < λ 1 , such that λ 1 is the largest Lyapunov exponent.
Of particular interest are the largest Lyapunov exponent λ 1 and the largest non-zero positive exponent λ N−2 . The largest Lyapunov exponent is the rate for the fastest growth of a phase space perturbation, and is dominated by the fastest dynamical events. The perturbation vector associated with λ 1 is strongly localized in phase space. The perturbations associated with the higher order Lyapunov indices are less and less localized, until the smallest positive exponent is associated with collective modes involving all particles.
The Kolmogorov-Smirnov entropy is related to the spectrum of the Lyapunov exponents, more precisely the KS entropy is equal to the sum of the positive Lyapunov exponents h KS = [76]. A first study of the Lyapunov exponents for the OGS system was performed for 3 < N < 10 by Benettin et al [74]. They found an increasing trend for the KS entropy and conjectured a linear increase with the number of sheets h KS ∼ N. In view of the results of [72] these values of N are too small to allow the study of the approach to equilibrium.
This study was extended to OGS with larger numbers of sheets 10 ⩽ N ⩽ 24 by Milanovic et al [77], in the more general setting of a 2-body power-like interaction |x − y| α . For the OGS case α = 1 they computed the entire spectrum for N = 10, 16, 24. For 32 ⩽ N ⩽ 64 only the largest two Lyapunov exponents λ 1 , λ 2 have been computed. They confirm the decreasing trend of the largest Lyapunov exponent with N, but find a change of the N dependence of the KS entropy at about N = 10. For N > 24 they find that h KS grows approximatively linearly with N. The regime change at the critical value of N = 10 is in agreement with the results of Reidl and Miller [78].
Tsuchiya and Gouda [79] extended these results to larger number of sheets N ⩽ 256. Their computation of the Lyapunov exponents is still the state of the art for the OGS in the unbounded domain.
They find generally good convergence and stability of the Lyapunov exponents in the large time limit. The largest Lyapunov exponent scales like λ 1 ∼ N −1/5 , confirming the trend noted by [77] using a more restricted range of values of N. The scaling of the KS entropy is h KS ∼ N 4/5 , which is almost linear. See figure 7 for the plot of h KS /N ∼ N −1/5 . The characteristic time for the system to distribute itself over many states, as measured by the KS entropy, decreases as τ KS ∼ N −4/5 . This corresponds to the macroscopic relaxation time-relaxation to a quasi-equilibrium state, which does not yet correspond to the system reaching the equilibrium microcanonical distribution. Complete relaxation to the microcanonical equilibrium state occurs over the longer time scale of microscopic relaxation time and signals approach to ergodic behavior. Such relaxation was indeed observed numerically, and occurs over much longer time.
Tsuchiya and Gouda [79] associate the microscopic relaxation time with the inverse of the smallest positive non-zero Lyapunov exponent λ N−2 . In other words, the microscopic relaxation time is determined by the time scale of the weakest instability, which is determined by the smallest positive Lyapunov exponent. The rate of decrease of this exponent to zero as N → ∞ measures thus the rate of growth of the relaxation time of the system to an equilibrium state as the number of sheets N increases.
Numerical simulation [79] suggests that the Lyapunov exponent λ N−2 approaches zero as λ N−2 ∼ N −1 , see the dashed-dotted curve with empty triangles in figure 7 (at least over the range N ⩽ 256 considered in [79]). This suggests that the relaxation to the microcanonical equilibrium state has a characteristic time scale of the order O(N). Further investigation of this scaling to larger values of N would clarify whether this asymptotics holds for all N.

Small number of particles (N < 10).
The best studied case corresponds to N = 3, which was considered both under the simple non-relativistic Newton dynamics, and by taking into account relativistic effects. Also, the cases of equal and non-equal masses have been considered.
We consider in detail here the non-relativistic N = 3 case with equal masses. It was shown by Lehtihet and Miller [80] that the system of three impenetrable sheets is equivalent to a particle moving in two dimensions in a symmetric wedge of opening angle 2θ with θ < 45 • under the effect of an uniform gravitational field. The equivalent wedge system is a particular case of a gravitational billiard. This equivalence is obtained by a change of variables for the differences of coordinates The momenta p i are also expressed in terms of conjugate momenta p α,β as Expressed in coordinates (α, β), the OGS Hamiltonian (4) takes the form where the interaction term is proportional to The contours of constant value of the interaction energy V int (α, β) are shown in figure 8. This is a constant force field pointing towards the origin in each angular sector of a hexagon in coordinates (α, β). A particle started at some arbitrary point with an initial velocity will move on parabolic trajectories within each sector, feeling a force pointing towards origin. The equal mass case is special as the distinction between transparent and impenetrable sheets disappears. For this case the wedge angle is θ = 30 • . We discuss this case in some detail in the next section. Consider the dynamics of a particle moving in a vertical constant field g, and constrained to the interior of a symmetric wedge with opening angle 2θ = 60 • . The equations of motion are The energy is Following [80] we choose g = 1 2 . We will study the trajectories with energy E = 1 2 for which the energy constraint takes the form The fixed point trajectory has parameters where y is the coordinate of the collision point with the wedge, and (v x , v y ) are the velocity components at the collision points. The collisions are normal to the walls. The most general e = 1 trajectory can be described by two parameters which are chosen as (ε, δ), corresponding to initial parameters, just before collision with the left wall The fixed point trajectory (81) is reproduced by taking the initial conditions ε = 0, δ = 0 • . The allowed range of values for ε is [−3/5, 2/5]. The lower bound comes from the condition y > 0 and the upper bound from v 2 > 0. Close to the lower boundary, the particle starts very close to the wedge vertex. The angle δ = 0 • corresponds to normal initial collision. δ > 0 means that the particle points up from the normal when it collides with the left wall, and δ < 0 means that it points down.
The trajectories are conveniently examined by plotting the Poincaré surface of section (or simply Poincaré sections). They show the values of the velocity components (v ∥ , v 2 ⊥ ) at each collision, with the right/left walls. Following the conventions of [80] we denote x := v ∥ , z := v 2 ⊥ the coordinates of the Poincaré section plot.
The energy constraint e = 1 gives that the section is bounded within the curve which gives z ⩽ 1 − x 2 . This curve is shown as the blue curve in figure 9. Along this curve the particle hits the left wall very close to the vertex y = 0. We distinguish between the cases of positive and negative ε, the perturbation of the initial condition from the fixed point (81).

Negative ε < 0.
For ε = 0 the trajectory is the fixed point shown as the red dot in figure 9. As ε decreases to negative but small values, the Poincaré sections are closed curves which are the sections of KAM tori.
As ε becomes more negative, the initial collision point with the left wall approaches the vertex y → 0. This is approached as ε → −0.6. In this region the system starts to display chaotic behavior.
Let us examine in more detail trajectories which comes closest to the vertex. Figure 10 shows the Poincaré section corresponding to ε = −0.595, for which the trajectory comes very close to the wedge vertex. We observe an area filling region at the corners of the triangular shaped section. In this region the system displays confined chaotic behavior.
All T a -type trajectories are above the Γ 1 curve, and all T b -type trajectories are below this curve.
As the collision velocity decreases, the T b -type movement becomes more and more likely (consecutive collisions on the same wall). As mentioned above, all such trajectories are below the Γ 1 curve. Figure 11 shows the Poincaré sections for ε = 0.1, 0.2, 0.3 and figure 12 shows an extreme case ε = 0.39 where the velocity at the initial collision with the left wall is very small ∼0.01.
In conclusion, there are two types of one-step trajectories, which define maps T a , T b . The T a maps correspond to repeated collision with the same side of the wedge, while the T b maps correspond to collisions with the opposite sides in succession. There is stable and chaotic behavior associated with finite period sequences of T a and T b . These regions in phase space are surrounded by a region of global chaos which contains trajectories passing arbitrarily close to the wedge vertex, and are thus near triple collision trajectories in the N = 3 equal mass OGS.
Under the classification of [81], the trajectories for N = 3 are of three main types: annulus, for which each particle crosses the other two in succession, pretzel, for which two particles cross each other at least twice before either crosses the third particle, and chaotic trajectories, for which there is no general pattern.
Similar behavior is observed also for the unequal mass case corresponding to θ < 45 • . For θ = 45 • the motion is completely integrable, and a dense set of orbits in phase space are densely filled by periodic points.  The role of the pretzel-type trajectories has also been studied in [82]. A numerical study showed that two particles which are crossing can form a long-lived pair, referred to as a 'molecule', which is eventually destroyed by interaction with the remaining particle. The interaction with the third one is not predictable and the observed long lifetime of the 'molecule' may prevent the system from being ergodic. Notice that at each time the central particle can be displaced between the remaining pair without affecting either the momentum or energy. By randomly displacing it at regular times during a simulation and taking a long-time average, it was demonstrated that one recovers Rybicki's micro-canonical density and velocity distribution for the three particle system (see [82] for details).
The relativistic version of the OGS has been studied in [83] for N = 3. They find two broad categories of periodic and quasi-periodic motions (annulus and pretzel-type patterns), as well as a set of chaotic motion that appears in the region of phase space between the former trajectories. Thus the global structure of the phase space remains similar to that of the N = 3 non-relativistic system.
We mention in passing that the gravitational wedge has seen interest in experimental and computational physics. It has been employed in optical lattices (see [84] and the numerous citations) and a driven, dissipative version has been explored experimentally [85] and computationally [86]. The ergodic theory for a range of angles and a family of related systems has been explored by Wojtkowski [87]. In particular he proved that for a wedge half-angle greater that π/4 (what he refers to as a a 'fat billiard') there is a single ergodic component. An extension to a conical geometry has also yielded interesting results [88].
A similar picture holds for N > 3 but not too large: the N particle OGS can be mapped to the system of a single particle moving in N − 1 dimensions, in a linear potential whose equipotential surfaces are those of a N − 1 simplex. The N = 4 case was studied in detail in [81]. The trajectories can be again periodic and quasi-periodic, and of chaotic type. Using a Braid Group representation, the periodic and quasi-periodic trajectories can be classified into three main types A,B,C, of which A,B are similar to those for N = 3, and C correspond to a new pattern where 2 pairs of sheets cross each other in succession. (N > 10). In contrast to the situation for N < 10 where the phase space of the OGS is segmented into confined and unconnected sectors, one chaotic and at least one with regular dynamics (period or quasi-periodic), the situation for N > 10 appears to be different.

Large number of particles
Reidl and Miller [78] have presented evidence that for N > 10 there is no segmentation of the phase space and the OGS may be ergodic. Numerical simulations by Milanovic et al [77] confirmed this by demonstrating by numerical simulations that the momentum distribution and the particle distribution converge to the equilibrium microcanonical results for the critical case N = 10. For N > 10 they find that after a few million characteristic periods of oscillations, the system relaxes in a state with a well-defined maximum Lyapunov exponent, and a well defined KS entropy per particle h KS /N.
The convergence properties of the largest Lyapunov exponent λ 1 for N < 32 (figure 6 in [77]) show that even when reaching equilibrium, there exist regions in phase space that do not contribute significantly to the growth of small perturbations. An arbitrary phase space trajectory spends a considerable time in these regions. This is seen in simulations as drops in λ 1 during the averaging process, followed by slow recovery as the trajectory moves back into the chaotic region. This picture was corroborated in [46,47,50,72].

One-dimensional gravitational systems with periodic boundary conditions
Several extensions of the OGS have been proposed, motivated by applications to structure formation in cosmology, which replace the free boundary condition with periodic boundary conditions. The main model of this type is the one-dimensional static cosmology (OSC) model, introduced by Aurell et al [22,89,90] and studied further by Valageas in [91,92].
We describe here in some detail the Miller-Rouet model [93] which is similar to the OSC model but differs in the details of the implementation of the periodic boundary conditions. In contrast, to the OSC model where periodicity is imposed by adding an external potential, the Miller-Rouet model maintains translation invariance and implements the periodic condition by modifying only the two-body interaction potential.

Miller-Rouet model.
A model for one-dimensional gravitating systems with periodic boundary was proposed by Miller and Rouet in [93]. This corresponds to the following setup: the system is enclosed in the box [−L, L) and periodic boundary conditions are imposed, identifying the points −L, L. In addition, the system is assumed to be placed into the uniform background of a mass distribution. The model is appropriate for studying one-dimensional density fluctuations in a uniform mass distribution. A version of the model replacing the gravitational interaction with Coulomb potentials has been used to investigate single-component plasmas in [94].
In contrast to the OSC model where periodicity is imposed by adding an external potential, the Miller-Rouet model maintains translation invariance and implements the periodic condition by modifying only the two-body interaction potential.
The interaction potential in the MR model has the form This is the potential energy of a mass at position x due to the interaction with another particle at y and all its mirror images separated by the periodicity length 2 L. We recall briefly the derivation of this potential and its relation to one-dimensional gravitation. The potential V(x, y) is the difference of two terms: the sum of the contributions from mirror images V 0 (x, y), and the contribution of the uniform background of mass Φ(x) The interaction V 0 (x, y) gives the potential felt by a particle placed at y from a particle at x plus the infinite number of its mirror images, separated by 2 L in both directions The damping factor e −κ|x−y+2kL| with κ → 0 is introduced following [10] and renders the sum over mirror images convergent. The sum over mirror images can be evaluated in closed form with the result where the first term is the contribution from the n = 0 term in the sum, and the remaining terms are the contributions from the mirror images of the particle at x. Expanding (88) in the limit κ → 0 and keeping only the terms which do not vanish in this limit we get The first term is the original linear attractive interaction, and the second term is a quadratic repulsive interaction, which vanishes in the limit of a very large periodicity radius L → ∞. The physical meaning of this repulsive term is as follows. As two particles are separated by more than L, the attractive effect of their mirror images in the nearby cells overcomes the attractive interaction between them. This appears as a repulsive force when the distance satisfies |x − y| > L.
Finally we subtract the contribution of a uniform background of mass. This amounts to a interaction energy Φ(x) given by This is a uniform potential, independent of position. The effect of subtracting the uniform background contribution Φ(x) from (89) amounts to canceling out the (positive) constant term 4πGm 2 /(2κ 2 L). The remaining constant term − 1 3 2πGm 2 L is negative and finite.  [93]. Each particle feels the effect of the other particles in the primitive cell and of their mirror images reflected through the boundaries of the primitive cell. There is an infinite number of images, and their total contribution is appropriately regulated with an exponential damping of the one-dimensional gravitational interaction due to Kiessling [10].
The Hamiltonian of the system is As explained, the effect of the mirror images is to introduce a repulsive interaction term quadratic in the particles separation. The chaotic dynamics of the system of N = 3 particles was studied in [95] in order to determine whether the separation of chaotic and ergodic behavior observed for the free boundary case [80,83] holds also in the system with periodic boundary conditions. We summarize here some of the results obtained.
The interaction energy can be simplified by introducing rhombic coordinates. Their definition requires some care, due to the periodic nature of the system. Define the primitive cell [−L, L), with the boundary points identified. Particles crossing the boundaries reappear instantly at the opposite boundary.
Assume for simplicity that the center of mass is initially fixed at zero x 1 + x 2 + x 3 = 0. The dynamics simplifies by introducing instead of (x 1 , x 2 , x 3 ) the rhombic coordinates α, β defined as and β which is defined in a similar way by replacing Expressed in the rhombic coordinates the interaction term in the Hamiltonian (91) takes the form The chaotic dynamics can be investigated by constructing Poincaré surfaces of section. Recalling that the phase space of the N = 3 system is three-dimensional, fixing one coordinate gives two dimensional Poincaré plots, which are easily interpreted visually. The additional constraint imposed is α + β = 0, which corresponds to crossing of particles 2 and 3. The resulting Poincaré plots are shown in figure 13, following [95], for two values of the total energy H = 0.226 (panels A-E) and H = 0.624 (panels F-H).
In contrast to the three body gravitating system with free boundary conditions, the Poincaré plots show dependence on the total energy. We start by analyzing the small total energy case H = 0.226. The Poincaré section is shown in figure 13(A). We note three major stable regions in the central and the upper left and right portions of the plot. The fractal region, located in the lower part of the plot, consists of self-similar sets of nested stable islands and includes infinite 'period-N' orbits that are surrounded by quasi-periodic orbits between which narrow chaotic regions exist. Figure 13 This picture is similar to that observed in the system with free boundary conditions, discussed in [80,83]. This can be understood by noting that for low energies, the inter-particle separations are very small, and therefore, the quadratic terms in the potential may be disregarded.
As the total energy increases, the gravitating system starts deviating from the free-boundary behavior. In the free-boundary version of the gravitating system, changing energy does not bring about a change in the structure of the phase-space, it just scales the phase space. However, in the periodic version, as energy is increased, the phase space gets more chaotic. However, small stable islands start appearing, as seen in figure 13(F) for H = 0.624, which grow and finally engulf the chaotic region. At higher energies, the particles are able to cross between rhombic planes (or equivalently, adjacent cells in the periodic 3-body system). Plots of primitive cell evolution and motion on the rhombic plane are shown for the central P2 orbit in figure 13(F) in figures 13(G) and (H) respectively. Another set is shown for a different P2 orbit (lying at the bottom of the horseshoe) in figures 13(I) and (J) respectively. Also note that periodic orbits form 'closed loops' (for the periodic boundary conditions, this means that after a finite number of strands on the rhombic plane, the trajectory will simply repeat on top of each other), quasi-periodic orbits result in 'bands' on the rhombic plane that trajectories will never move out of, and chaotic orbits lead to unpredictable trajectories on the plane. Of course, conservation of energy may set a boundary beyond which the particle will never go, in which case, the chaotic trajectory will fill in the allowed region as time progresses). As one moves away from a fixed point (a 'period N' trajectory, N ∈ I) on the Poincaré plot, quasi-periodicity takes over and the strands on the rhombus starts spreading into bands.

Lyapunov spectrum for 1D systems with periodic boundary conditions.
A new method for computing the Lyapunov spectrum was developed and applied to the study of the 2component plasma with periodic boundary conditions in [94]. This was extended also to onedimensional gravitational systems with periodic boundary conditions in [97]. We discuss here in some detail the case of the gravitational system.
The paper [97] presents a method for computing the entire spectrum of Lyapunov exponents λ i . The initial conditions are chosen such that the sheet velocities are Gaussian distributed, while their positions are chosen randomly. The system is allowed to relax over a sufficiently large time such that the Lyapunov exponents converge to well defined values. The number of sheets was taken as 5 ⩽ N ⩽ 20.
The results for the Lyapunov spectrum satisfy the expected properties: (a) the sum of all exponents converges to zero, (b) the Lyapunov exponents appear in pairs of opposite sign λ i = −λ 2N−i+1 , and (c) four Lyapunov exponents are consistent with zero λ N−1 ≃ λ N ≃ λ N+1 ≃ λ N+2 ≃ 0. The spectrum depends on the energy per particle.
Here T is the temperature and ρ the gas density. This result is important in cosmology where it gives the conditions under which density inhomogeneities can grow and form dense objects.
The assumption of equilibrium of the infinite homogeneous system is not completely justified, due to the absence of pressure gradients which would balance the infinite gravitational problems. In order to proceed, one must accept the existence of the equilibrium, and postulate that the gravitational potential is determined by the density fluctuations. This assumption is the so-called Jeans swindle, see section 5.1 in [3] for a discussion. The assumption has to be checked for each case, and ensure that it leads to consistent predictions.
One solution to the Jeans swindle was proposed by Kiessling by introducing an exponentially damped regularization of the gravitational interaction [10]. The range of the exponential damping can be taken to infinity at the end of the calculation, and physically meaningful results are obtained for quantities for which this limit exists and is finite.
The simplest setting where this phenomenon can be studied is the OSC model of Fanelli and Aurell [89]. This consists of N particles on a line, interacting by 1D gravitational potentials 2πGm 2 |x − y|. The particles are restricted to a finite region [0, L], and continuity boundary conditions are imposed ρ(0) = ρ(L), but not the continuity of the first derivative.
The properties of this system have been studied by numerical simulation, and analytically in the Vlasov limit by Valageas [91,92]. Both static and dynamical equilibrium properties have been obtained. The system is found to have qualitatively different properties below and above a critical temperature. The equilibrium state above this temperature has a uniform density, while below this temperature the gas has a non-trivial density.
The OSC model breaks translation invariance, due to the boundary conditions used. This problem disappears in a related system, considered by Miller and Rouet [93], which is discussed in the next section. [93]. This model and its derivation were presented above in section 2.5.1. We repeat here only the main features, and discuss its thermodynamical properties.

Thermodynamics of the Miller-Rouet model. A model for one-dimensional gravitating systems with periodic boundary was proposed by Miller and Rouet in
The system is enclosed in the box [0, L] and periodic boundary conditions are imposed, identifying the points 0, L. In addition, the system is assumed to be placed into the uniform background of a mass distribution. The model is appropriate for studying one-dimensional density fluctuations in a uniform mass distribution.
The interaction potential in the MR model has the form 7 The thermodynamics of the MR model was studied in [19]. This system has a more complex behavior than the one-dimensional gas enclosed in a box described in the previous section. Above a critical temperature T c1 the gas density is uniform, while below T c1 the density develops a spatial density with one minimum. Decreasing the temperature further there is an infinite sequence of critical temperatures T c1 > T c2 > . . . at which further states become accessible, with more minima within the primitive cell.
We give an overview of these phenomena, starting with the gas density, which is given by the following result. (95) in thermodynamical equilibrium at temperature T satisfies the Lane-Emden equation

Proposition 4. The single particle distribution function of the gas of particles interacting by the two-body potentials
This is normalized aŝ We would like to solve the equation (96) with the constraint (97), for given temperature T. It is convenient to introduce the new unknown function y(x) defined by Expressed in terms of this function, the differential equation (96) reads with the normalization constraint We impose periodic boundary conditions At this point we note the similarity of the equation for the single particle distribution function in the Miller-Rouet model (99) with the corresponding equation in the OSC model [91], see equation (12) in [91] (up to a redefinition y(x) → −βψ(x) and rescaling x/L → x). There is also a slight difference, as the boundary conditions (101) are more constraining than the corresponding boundary condition in [91]. In particular, we require y(0) = y(1), which is not imposed in the OSC model [91]. As a result, although the qualitative properties of the solution is similar in both models, the details of the solution are different.
The qualitative behavior of the solutions of the equation (99) can be visualized in an intuitive way by writing it as where α 2 := 2βg 2 . Expressed in this form, the problem can be visualized in terms of a dynamical analogy: this is the same as Newton's equation of motion for a particle of mass 1 moving in the potential V(y). See figure 14 for a plot of V(y). The total energy is conserved Using this dynamical analogy it is easy to understand the qualitative behavior of the solutions of the differential equation (102). The equation (102) has always the trivial solution y(x) = 0, which corresponds to the particle sitting at rest at the bottom of the potential well V(y). In addition to this trivial solution it can have oscillatory solutions, corresponding to the particle moving in the potential V(y), starting at some non-zero value y(0) ̸ = 0 with a positive or negative initial speed y ′ (0), and then performing one full oscillation or several oscillations before returning to the starting point y(1) = y(0) with the same velocity y ′ (1) = y ′ (0) at time 1. The movement of the particle is spanned by y L ⩽ y(x) ⩽ y R , where y L < 0, y R > 0 are the turning points at which the particle speed vanishes. They are related by energy conservation to the initial position and speed as V(y L ) = 1 2 [y ′ (0)] 2 + V(y(0)) = V(y R ). The simplest solutions correspond to trajectories where the particle performs one full oscillation before returning to y(0) at 'time' x = 1. There are also solutions for which the equivalent particle performs two, three, etc oscillations. We will denote the solution corresponding with k oscillations the kth mode. They are obtained for different values of the α parameter: for 2π < α < 4π there is solution k = 1, for 4π < α < 6π there are two solutions with k = 1, 2, and so on.
The higher order solutions are related to the k = 1 solution as The analysis presented above gives the following qualitative behavior of the gas density as the temperature is lowered. In the infinite temperature limit T → ∞ we have α → 0 and the gas density is constant ρ(x) = 1. As the temperature is lowered, the density ρ(x) remains constant until we reach α = 2π when one non-trivial solution for y 0 appears. This point corresponds to temperature Compared to the critical temperature in the OSC model [91], this is smaller by a factor of 1 4 . This is due to our boundary condition y(0) = y(1) which is not imposed in [91]. However, the result for T c1 has the same dependence on model parameters as in the OSC model, see equation (11) in [92] which gives T c1 = 2g 2 π 2 in our notations. Note that in this reference 2πG is denoted g.
As the temperature is lowered below this point, non-trivial solutions with inhomogeneous gas density appear. They are translated versions of the basic solution ρ 1 (x) = exp(y 1 (x)). ρ 1 (x) has a maximum at x = 1/2.
As the temperature is lowered further, we reach the point α = 4π, corresponding to temperature T c2 = g 2 ML 8π 2 . Below this temperature there are two solutions for y(x) corresponding to k = 1, 2. In addition to the k = 1 solution we have another solution with k = 2, which has oscillatory density behavior, and has two maxima/minima within the box. In general there is an infinite sequence of critical temperatures at which new solutions appear, given by α = 2nπ, with n = 1, 2, . . .
The thermal properties of the gas are described in terms of the free energy F of the gas. This can be found in closed form and is given by proposition 2 in [19]. We discuss here some of the main properties of this function.
Around the critical point T c1 (α just above 2π), the free energy per particle can be approximated as where f 0 = T(log(N/L) − 1 + 1 2 log h 2 2πT ) is the free energy per particle in the homogeneous density phase, below the critical temperature, and we denoted α 2π = Tc1 T and x ≡ T/T c1 . From this expression it follows that the free energy difference f − f 0 and its first derivative with respect to temperature vanish at T = T c1 , while the second derivative has a jump from 0 Since the difference f − f 0 vanishes for T > T c1 , this implies that the free energy and its derivative are continuous at T = T c1 while its second derivative has a jump. We conclude that the phase transition at T = T c1 is a second order phase transition.
The specific heat is discontinuous at the critical point, and drops as the temperature increases and crosses the critical point T c1 . Above the critical point T c1 the specific heat is constant and equal to c V = 1 2 , and below the critical point it takes the value lim T→Tc1−ϵ The behavior of the specific heat close to the critical point was studied further in [19].

Concluding discussion
It seems appropriate to close the section on the approach to equilibrium by revisiting the three questions asked in the Conclusions section of the Reidl and Miller [42]: (a) Has thermalization ever been observed in a one-dimensional self-gravitating system? Several studies presented evidence for relaxation to the equilibrium state in the OGS [44,51,52]. Joyce and Worrakitpoonpon [51] studied in detail relaxation to the equilibrium distribution, as measured by the order parameters ϕ 11 , ϕ 22 , defined as expectations over the phase space of the general form ϕ αβ = ⟨|x| α |v| β ⟩ ⟨|x| α ⟩⟨|v| β ⟩ − 1.

(b) If thermalization occurs, what is the source of dynamic hyperbolicity (instability)?
Experience with systems with small numbers of particles suggest that thermalization is related to three-body collisions. These collisions introduce chaos in the N = 3 system, as discussed in section 2.

(c) If thermalization occurs, what is the time scale for relaxation?
Joyce and Worrakitpoonpon [51], working with OGS systems with N ∼ 10 3 particles, find evidence that the relaxation time to equilibrium is linear in the number of particles t R ∼ O(N), similar to that found earlier by Tsuchiya et al [44] and [40]. A similar scaling follows from the two-step relaxation mechanism involving a QSS in the intermediate state, where the relaxation time is given by equation (64), although the factor f QSS could introduce also mild dependence on N and the initial state. In the language of dynamical systems this time scale is reflected in the N scaling of the smallest positive Lyapunov exponent λ N−2 . Extrapolation in N of numerical simulations presented in [79] (see section 4.3) gives a similar N scaling. Levin et al [52] study the relaxation from a core-halo state to thermodynamical equilibrium. They present numerical simulations with N = 500, 750, 1000 sheets which give t M ∼ O (N 1.8 ), similar to the scaling ∼O(N 2 ) found by [37]. However, the metric used in this paper N(x, t) the spatial density and N ch (x) the core-halo density, is different from the metrics used by the other studies mentioned above. It is not obvious that their N scaling should be the same. In view of this apparent discrepancy, further study of this question is warranted.
In addition to the two-stage relaxation mechanism discussed at the beginning of section 2.3, evidence of a more complex pattern, hinting at transitions among different quasi-equilibrium states, was pointed out by Tsuchiya and Gouda [45], which refer to it as itinerant behavior. This phenomenon is still poorly understood and more work is required to understand the details of the relaxation to equilibrium.
We mention an interesting connection to atomic physics. The linear attractive interaction typical of 1D gravitational systems appears in other problems unrelated to self-gravitating systems. It was noted in [99] that in a gas of neutral cold atomic gas, the atoms interact with linear attractive potentials when placed in the radiation field of a quasi-resonant laser beam. This phenomenon has been studied experimentally in a cold Strontium gas. This allows the laboratory study of this type of systems, and of their non-equilibrium dynamics.

Generalities
Observations show that, on large scales, the Universe is composed of galaxies, galaxy clusters, superclusters, and large voids [21]. This hierarchical structure has led scientists to speculate that there is an underlying fractal structure to the distribution of mass in the Universe [100,101]. Numerous, but not all, investigations support this hypothesis [102,103]. While not the subject of this review, we mention that an important study concluded that there is an absence of fractal structure in the local Universe and homogeneity already at about 70 Mpc [104]. However more recent investigations take issue with this conclusion and provide evidence of robust fractal structure [105] and the absence of homogeneity on much greater scales [106]. To demonstrate the emergence of fractal behavior requires the ability to construct dynamical simulations that are precise on many length scales [107]. In the past this has proved challenging for three dimensional evolution codes. While less realistic, event-driven one-dimensional simulations can be carried out with high precision and have have provided valuable insights into the development of cosmological self-similar, fractal structures [108][109][110].
Although not the focus of this work, it is worth noting that, following recombination, dynamical considerations suggest that the first structures to form were approximately highly flattened, one-dimensional sheets referred to as 'Zeldovich pancakes' [20]. Thus, there is a possible connection with 1D cosmological models and the real Universe. However, since the pancakes are short-lived on a cosmological time scale, one should not take this connection too seriously, especially since it is the development of structures on cosmological time (ie, 1D 'galaxies') that interests us here.
In the previous section we saw that, when the temperature is lowered below a critical value, fluctuations about the mean density become unstable and the system undergoes a phase transition and reverts to an inhomogeneous state. This can be considered an equilibrium representation of cosmological structure formation. In actuality we know that the Universe is expanding. From relativistic cosmology we know that the expansion factor, a(t), is a solution of the Friedman equation. According to the cosmological principle, the Universe is homogeneous and isotropic on large scales. To understand how structure formation-ie, the creation of galaxies, galaxy clusters, and voids-occurs as the Universe expands, cosmologists focus on a finite segment of the Universe that is larger than any structures of interest, but small enough that Newtonian dynamics is adequate. It is useful to consider the dynamics in a frame of reference that keeps up with the expansion, ensuring that the mean density is constant. This is accomplished by working in the comoving frame. Since there is nothing special about the particular segment under consideration, it is typically assumed that it obeys periodic boundary conditions. It is thought that alternative choices of boundary conditions will have a stronger influence on the development of large scale fluctuations, which we want to avoid. However other boundary conditions have also been studied and yielded similar results [89,111,112]. From section 2.5 above, we recall that we have established a mathematical framework for accounting for periodic boundary conditions [93]. We will see below that one-dimensional models have provided numerous insights into the properties of structure formation in an expanding Universe.
The first one-dimensional models of cosmological evolution were carried out by Rouet and Feix, and, following Fanelli [89], we shall refer to them as 'RF' models [111]. In its simplest form, we have seen that a one-dimensional gravitational system can be thought of as a system of parallel mass sheets of infinite extent. Each sheet experiences the (constant) gravitational field of the remaining sheets (an electrical analogue would be parallel plate capacitors that could be viewed as concentric spherical capacitors with an arbitrarily large radius). In the RF model, only gravitational forces were considered. To take into account the cosmological expansion factor, a(t), it was assumed that, as time progresses, the surface density of each sheet decreases as 1/a 2 . We can introduce a coordinate, r, perpendicular to the surface of each sheet and consider an assembly of say N sheets in a cell of length L(t). Due to the expansion, L(t) is increasing as a(t).
To maintain a constant average density, the equations of motion were transformed into the comoving frame expanding with a(t). For simplicity, specular reflection of particles at the system boundaries was initially assumed. Later it became possible to asses the fields in the presence of periodic boundary conditions exactly and these were employed in all subsequent work [93]. Since only gravitational forces were considered, it was assumed that a(t) ∼ t 2/3 corresponding to a matter dominated cosmology [21]. In addition, the time was also re-scaled with b(t) ∼ t β . While the power law governing a(t) is determined by the expansion, the choice of the exponent β is free. It is shown below that β = 1/2 gives a rescaled equation of motion with time-independent coefficients. For this system it was possible to carry out a very precise N-body simulation by computing the next crossing time between each pair of mass sheets and updating the positions and velocities accordingly. The crossing times were obtained by analytically solving a cubic equation (see B.1). Because of the high numerical precision it was possible to follow the system for very long times and study the resulting fractal density and phase space particle distributions.
It was later demonstrated by Fanelli and Aurell [89] that the parallel sheet system could be considered to be a planar perturbation embedded in an isotropic, spherical three dimensional Universe with only a small change in the equations of motion. For the case where the forces are purely gravitational, the evolution could also be determined numerically by solving for the crossing times (see appendix B). However in this case a fifth order equation is obtained so it is commonly referred to as the quintic or Q model. Moreover, it was shown that the embedding concept could be extended to arbitrary dimension d to describe a family of models [113]. Also, in that framework it was possible to introduce the standard cosmological description of expansion including, for example, the cosmological constant and density parameters [114].
Because of the high precision with which the 1D simulations with point-wise particles can be evolved, they are ideal for studying fractal geometry. Most of these studies have focused on the configuration space [108,110,112,115,116] but the extension to the phase space has also been considered [116][117][118]. The results suggest that multi-fractal behavior is present as well [113,119]. However, the representation of the matter distribution is limited at low density by the granularity of the particles. To gain greater insight into multi-fractal behavior, a continuous or Vlasov representation of the system has also been investigated [120].
The 1D systems have been used to investigate and compare gravitational clustering in several other cosmologies besides the matter-dominated cases mentioned above. The Dirac-Milne cosmology [23] (hereafter DM) assumes a Universe where there are equal amounts of matter and antimatter, which separated before annihilation could eliminate the latter. To maintain the separation, it is assumed that the antimatter repels both itself and regular matter. An advantage of this cosmology is that it does not require the assumptions of inflation, dark matter or a cosmological constant. Recently, the 1D model was employed to investigate a version of the DM cosmology where the antimatter component was uniformly distributed as a background. A comparison with a 1D version of the ΛCDM model was also carried out [114].
In the following we will consider in some detail each of these investigations. For mathematical economy we will introduce a formalism where the 1D system of sheets is treated as a perturbation embedded in an isotropic system of d-dimensional shells. Where it is appropriate we will introduce the standard cosmological constants that acquire their usual meaning for the case d = 3.

Expanding Universe in comoving coordinates
We consider a spherically symmetric Universe in d dimensions, where all quantities depend on a single spatial variable r, plus the time t. We shall use Newtonian gravity, but include the presence of a finite cosmological constant. The equation of motion of a spherical shell in such a Universe is as follows: where E r = −∂ r ϕ is the gravitational field, c is the speed of light, and Λ is the cosmological constant. The gravitational field is a solution of the Poisson equation where ρ(r, t) is the matter density. We now transform the equation of motion to the generalized comoving variables, denoted by an over-caret and defined by the transformations where a(t) is the dimensionless scale factor of the Universe. The velocity transforms as where the dot stands for time differentiation with respect to t. Using these transformations, the equation of motion becomes whereÊ(r,t) is the scaled gravitational field. As the density must scale asρ(r,t) = a 3 (t)ρ(r, t) in order to preserve the total mass, we scale the gravitational field asÊ(r,t) = a 2 (t)E r (r, t), so that Poisson's equation remains invariant in the scaled variables: In equation (114), it is convenient to choose a scaling that satisfies b 4 = a 3 , so that the coefficient in front of the gravitational field is time-independent. However, other choices can be found in the astrophysical literature [121]. This choice yields: Most cosmological models assume a uniform distribution of matter at large scales. Hence, it is useful to establish whether the comoving equation of motion (116) admits stationary solutions that are uniform in space. A stationary solution in the comoving variables corresponds to an expanding solution in real space, with expansion factor a(t).
For a constant densityρ = ρ 0 , Poisson's equation (115) can be solved exactly in a d-dimensional space to give the gravitational field:Ê r = −4πGρ 0r /d. Substituting into equation (116) and considering an equilibrium configuration (ie, setting all time derivatives to zero), we obtain the condition: We define the density parameters, i.e. the densities normalized to the critical value ρ c = 3H 2 0 /(8πG), as: Rearranging terms, equation (117) becomes: For d = 3, this is nothing but the second Friedmann equation for a non-relativistic Universe without radiation pressure. The first Friedmann equation can be obtained by integrating equation (119) once, which yields: where Ω T = Ω M + Ω Λ and Ω K = 1 − Ω T is the spatial curvature energy density. All in all, we have shown that the Friedmann equations can be recovered in a fully Newtonian context simply by requiring that a homogeneous steady-state solution exists in comoving space.
Substituting the Friedmann equations (119) and (120) into the comoving equation of motion (116), one gets, after some cancellations: Subsequently, we consider locally planar perturbations embedded in this expanding Universe and denote the corresponding comoving coordinatex. In this locally planar system, Poisson's equation can be written in a 1D form: ∂ xÊx = −4πGρ(x,t). In order to use periodic boundary conditions, we incorporate the first term on the rhs of equation (121) into the 1D Poisson's equation, yielding the following system of equations: Note that, while for d = 1 the Poisson equation (123) is consistent with the equilibrium solution ρ = ρ 0 , for d > 1 we had to artificially remove a factor 1/d in front of the equilibrium density on the rhs of equation (123). This is a small artefact that occurs because of considering 1D perturbations embedded in a 3D Universe. Three cases are particularly interesting: • Einstein-de Sitter Universe: Ω M = 1, Ω Λ = 0, and Ω K = 0 (flat-space Universe), for which equation (119) leads to a(t) ∝ t 2/3 . In this case, all coefficients in equation (122)  In the remainder of this section, we will illustrate how the above 1D gravity models can be used to study structure formation in a cosmological context, with particular emphasis on the Einstein-de Sitter Universe.

Einstein-de Sitter model
Historically, the Einstein-de Sitter model was adopted to perform the earliest cosmological simulations. It corresponds to a matter dominated Universe and represents the critical case separating an infinitely expanding Universe and a periodically evolving one. Rescaling methods, such as those used above to introduce the generalized comoving coordinates, are well-suited to tackle this kind of evolution [122,123]. Moreover, the Einstein-de Sitter solution is selfsimilar, the scale factor going like t 2/3 . This has some importance as the time evolution of the system never stops. In particular we will show that the clustering process which is observed in the simulations (provided the excitation wavelength is small enough) continues as long as the physical limit of the system does not play any role. Therefore it is a good candidate for performing a fractal analysis. As indicated before, this model corresponds to Ω M = Ω T = 1, Ω Λ = 0. With these values, equation (118) give where ω 2 J0 is the Jeans frequency defined at present time. Then equation (122) reads: Depending on the value chosen for d, equation (125) defines a family of models [113]. The RF model is obtained by taking d = 1, while the Q model requires d = 3. Obviously it is also possible to consider d = 2 for which the particles of the physical system are expanding cylinders. Let us consider the source of instability that initiates cluster formation. Note that as d is increasing, the friction coefficient term is decreasing. So in more than three dimensions, d could be greater than 3, and taking d going to infinity will reduce the friction to 0. Then the system dynamics is conservative in the comoving frame. This case has also been investigated with numerical simulation [108]. In that case, the one-dimensional perturbation of the rescaled system gives exactly the opposite of a one component plasma, that is attracting particles in a repulsive fixed and neutralizing background. Taking into account the similarity of the two systems, the dispersion equation reads which indicates that the system is unstable for wavelengths l > 2πv T /ω J . For such a system, the Jeans length is given by the limit value l J = 2πv T /ω J . If the system is initially perturbed with a wavelength greater than l J , it dominates quickly as it is the most unstable. Then we observe the formation of clusters the size of which is driven by l J , so by the initial thermal velocity v T . Even in the case where the friction does not vanish, the same behavior is observed. The Q and RF models and the model without friction have been studied extensively using numerical simulation [90,111]. Two principal approaches have been employed to build the initial conditions. The earliest simulations employed a water bag distribution where the particles are equally spaced and the velocities are chosen by uniformly sampling a symmetric interval [124]. Later on, the Harrison-Zeldovich approach was employed where a power law was assumed for the Fourier representation of the density [112,125]. In particular, in [108], the choice ρ k ∼ k 3 was selected because of the correspondence with ρ k ∼ k in three dimensions that is close to the observed values given by studies of the cosmic microwave background (CMB) [108,113,119,120].
In the following, as an example, we will focus on the RF model. Figure 15 shows the time evolution of the density and the distribution in µ-space (position, velocity) for N ∼ 2.5 × 10 5 particles enclosed in a periodic box. The units are normalized such that the field created by a single particle is ±1/2, and the average density is 1 (so the length of the system is equal to the number of particles). Accordingly, the time is given in these units, andt = 0 should be considered as the 1D analogue of the cosmological recombination time. The particles are initially localized so that the density spectrum goes as k 3 and the velocities are chosen to select the growing mode [90,108].
As the system evolves, a hierarchical structure formation is observed with a bottom up scenario. The initially homogeneous system fractures into groups of particles, which fracture again to form subgroups, and so on. This process continues as long as the boundaries do not have any effect. It is almost achieved att = 30 where only three large clusters remain at that time. That is too few to claim that these structures are uniformly distributed on the system size. To obtain better insight concerning the structure formation att = 30, successive zooms are plotted in figure 16. From top to bottom, each picture gives a zoom of the region defined by a rectangle drawn in the phase space picture just above it. We observe a lack of mixing in the overall box by attributing a color to an ensemble of nearby particles att = 10. The particles retain their color at each time in figure 17. It is shown that the clusters only gather nearby particles. Both the density spectrum and the correlation function then give a more global picture of the self similar distribution of the particle positions. Figures 18 and 19 show the time evolution of, respectively, the density spectrum and the correlation function. As time goes on, a scaling range appears and increases in each. In log-log plots the density spectrum yields a slope of −0.48 whereas the correlation function exhibits a slope of −0.56.
Note that the power spectrum shows two distinct scaling ranges separated by a peak at the crossover, say k c (t). At long wavelengths, i.e. for k < k c , the slope in the log-log plot retains the power-law associated with the initial condition. This represents low-density regions where the particles' trajectories have not undergone crossings with each other. On the other hand, the scaling region for k > k c manifests the effects of cluster formation. This power-law behavior results from the scale-invariant properties of the clustering phenomena which we have discussed earlier and provides a 1D version of a fractal Universe. It is natural to consider that the crossover wave vector k c (t) represents the typical inverse scale for the largest clusters formed at that time. As the non-interacting regions decrease in size, k c (t) moves to smaller and smaller wave vectors so the initial conditions only influence the very largest scales. The determination of the cross-over scale has been discussed in [108,110,126]. The first quantitative indications that scale-free, or fractal, behavior was present in the distribution of galaxies came from studies of the galaxy-galaxy pair correlation function. This required the ability to determine proper as well as angular inter-galaxy separations from large scale surveys, such as the sloan digital sky survey (SDSS) [127]. Here as well, in our 1D simulations, power-law decay is a manifestation of fractal geometry. However, as a result of Baryonic Acoustic Oscillation, additional structure shows up in the astrophysical observations [21]. In addition to the power law decay of the correlation function and power spectrum, the existence of a single scaling exponent suggests that the mass distribution is becoming fractal on larger scales as time progresses.
Box counting is the most popular approach for analyzing the properties of a fractal structure and here we use it to determine the generalized fractal dimension D q [107,128] of the particles in configuration space. In order to do so, the system (0 ⩽x ⩽ L) is covered with boxes of length ℓ of decreasing size: ℓ = L/2, ℓ = L/4, and so on. Then the quantity I(q, ℓ) = i (N i /N) q (ℓ) is computed for each value of ℓ, where N i is the number of particles contained in the ith box and q is an integer which is intended to give more weight to high-density (when q > 0) or low-density regions (q < 0). The generalized fractal dimension is defined as for a given value of q ̸ = 1. In the limit q → 1, equation (127) allows us to define the information dimension D 1 : Figure 17. For a small spatial slice of the simulation given in figure 15, clusters are painted in different colors att = 10 and then drawn at different earlier times.
One also defines D 0 , called the box counting dimension and D 2 the correlation dimensions.
The following relation holds between the index n of the power spectrum P(k) ∼ k n and the correlation function ξ(x) ∼ x −γ [108]: A log-log plot of ln[I(q, ℓ)]/(q − 1) as a function of ln(ℓ) gives a scaling range of slope D q as seen in figure 20. The box counting dimension D 0 has been computed both in configuration space and phase space (in that case, the µ-space is covered with rectangles). Because the data of real observations show a more complex structure than a uniform fractal, the generalized dimension has also been computed using size-oriented methods, such as the common box counting method or the point-wise method. In these methods the parameter q allows one to distinguish the high density (q ≫ 1) and low density regions (q ≪ 1). However these methods show a non increasing D q spectrum, which is not what should be expected [129,130]. At t = 30, figures 18-20 exhibit a strong scaling range. They give the values n = −0.48, γ = 0.56 and D 2 = 0.44 which agree quite well with relations (129). These exponents have also been computed for different models and initial conditions [108,110] Figure 20 shows that the scaling ranges are not so well defined as q decreases and have to be carefully chosen [131]. Consequently Shiozawa et al explored alternative methods of computing the generalized fractal dimension introduced by Grassberger et al [132] and extended by van de Water and Schram [133], based on equal mass partitions [134]. They are more sophisticated than standard methods based on equal size boxes, in the sense that they do not give directly the D q spectrum. Nevertheless, as for the previous method, a scaling range of a curve in a log-log plot has to be determined. This scaling region has to be chosen carefully. All these methods have been carefully tested on arrays of particles distributed according to well known fractals such as the Binomial Multiplicative Process or the generalized Cantor Set [107]. The result is that, at least for these two fractals, they are able to provide insight regarding the D q spectrum especially for negative values of q.
The usual box counting method was originally used, both in configuration space and µspace to determine the box counting dimension D 0 for the RF and the water-bag initial conditions. Later the box counting as well as the point-wise dimension was used to give a picture of the overall spectrum D q for the RF and Q models for an initial power law power spectrum [108]. As discussed above, the correlation dimension value D 2 is also recovered by the slope taken from the density power spectrum and the correlation function itself (at the same time). Nevertheless, these methods give an increasing D q as a function of q which is unacceptable [129,130]. To obtain better insight on the negative interval values of q, the two mass-oriented methods previously described have been used. All the methods give the same results for positive values of q, that is almost D q ∼ 0.4, independent of q and close to the description of a monofractal (see figure 25). But they differ for negative values of q for which a decreasing D q is now observed and suggests the existence of multi-fractal geometry [113]. This difference of behavior suggests that the low density regions are not well treated by methods based on partitions of equal size and that a Vlasov approach may overcome this problem.

Two-component model
By observing the distribution of luminous matter in galaxies and clusters of galaxies we are able to extract information concerning the fractal nature of the mass distribution in the Figure 19. Time evolution of the correlation function of the mass distribution for the simulation given in figure 15. The red arrows define the region over which the slope is calculated att = 30.
Universe. However, according to current theory, the visible matter only accounts for roughly 20% of the total matter, and we have no way of directly observing the dark matter component, although its presence can be inferred from gravitational lensing and galactic rotation curves. An important question then arises: What is the fractal geometry of the total matter distribution of the Universe or, equivalently, is the geometry of the dark matter component different from that of luminous matter?
Theoretical cosmology has not provided us with an analytical prediction of the geometry of structure formation and investigators have resorted to computational simulations to gain insight into these questions. Until recently, structure formation has been modeled by a single component system which does not distinguish between the components-they are essentially dark matter simulations. An essential difference between luminous, or baryonic, matter and dark matter is the presence of dissipation in the former through radiation and other processes. In order to explore possible differences between each matter component in a 1D context Shiozawa et al constructed a two-component self-gravitating model where one of the components can lose energy [119]. Of course, in the comoving frame, there is already an apparent dissipative contribution created by the friction term that arises from the transformation from inertial to comoving coordinates. However, what is needed is a dissipation mechanism that only applies to one of the matter components.
In the 1D dynamics dissipation can be conveniently included by introducing energy loss in one component whenever a pair of 'matter sheets' cross. In common with inelastic collisions that occur between material objects in the laboratory (steel ball bearings, for example), one can introduce an effective restitution coefficient at each sheet-crossing. A velocity-dependent collision coefficient κ was introduced analogous to a restitution coefficient. Following each crossing of the luminous component, the velocity of each sheet in the pair is reduced by a factor κ. The velocity dependence is given by κ = exp(−c|v 1 − v 2 | 3/5 ) where v 1 and v 2 represent the velocities of two colliding particles or sheets. The coefficient c was chosen arbitrarily in the simulations so that the trajectories of luminous matter particles are substantially different from dark matter particles without forcing them to collapse too quickly. The luminous particles lose more energy when the velocity difference between the two crossing elements is large.
The initial conditions were chosen following the 1D version of the Harrison Zeldovich prescription as in the simulations discussed earlier. Simulations were carried out with 10 5 particles and a fifth of them were chosen to represent dissipative, or luminous, matter. To illustrate what happens, in figure 21 the evolution of the positions of a system with 300 particles is tracked just for clarity, so 60 particles are subject to the velocity dependent collision law. At the onset the tracks are nearly vertical, indicating that the velocities are small and nearly constant. As crossings occur, fluctuations in the local field develop and the tracks change direction. As anticipated, the development of self similar clustering occurs with dark matter forming a 'halo' around a core of luminous matter particles. In figure 22 the distribution of matter in phase space at three selected times in the evolution is displayed as well as a plot of the potential, kinetic and total energy showing energy loss starting after about ten time units. Strong clustering is apparent here as well. Employing a simulation with 10 5 particles, the positional data was subjected to a complete mutifractal analysis. In order to avoid the known difficulties associated with equal size decomposition (see above), the mass-oriented k-neighbor and fixed k methods were employed, where k is the number of neighbors separating a pair of points. For a discussion and comparison of these methods with known fractals see Shiozawa et al [134]. In general, at a given time, three scaling ranges were observed as a function of ln k. The highly clustered matter contributes strongly at small values of k, whereas for very large separations, the system appears homogeneous. An intermediate range represents the transition from the clustered to the homogeneous scale.
It was possible to extract generalized dimensions D q at specific times from the analysis, and specific values are provided in table 1. The table indicates that, in contrast with equal size decompositions, here we find that D q is a decreasing function of q as it should be. Let us just focus for now on D 0 , and D 2 , the box counting and correlation dimensions.
The values of D 0 , and D 2 for dark matter at later times agree with earlier studies of single component systems using the more traditional equal size decomposition [108]. In contrast, as time progresses we see from the table that the luminous matter has consistently lower generalized dimension than the dark matter, suggesting denser, more compact clusters. The similarity of D q for negative q values for each component suggests that the low density 'voids' have no difference in structure.  At the beginning of the simulation, fractal analysis correctly shows that the model Universe is homogeneous on all scales. In common with 3D simulations, however, as time progresses, in-homogeneity continuously grows and the scale at which the distribution becomes homogeneous expands over time. We see that each type of matter follows a similar evolution but with significant quantitative differences. As we would naturally expect, due to dissipation, luminous matter is concentrated at the core inside clusters. This gives rise to lower fractal dimensions in the positive range of q. Due to the energy loss during the collisions, luminous matter starts to coalesce first. In contrast, the void regions show no significant difference in structure between the two types of matter as the long range force is primarily responsible for void formation. The structures continue to increase in size and the self-similar patterns, with well-defined fractal dimensions, persist over time. The difference in fractal dimensions within clusters is a manifestation of the bias of the luminous matter distribution against the dark matter distribution.

Vlasov dynamics of structure formation
A large amount of work on the formation of complex gravitational structures in our Universe has been accomplished using Newtonian 1D models [108,110,112,113,135]. However, Nbody simulations suffer from an intrinsic undersampling of the phase space because of the finite number of particles used in the codes. As the number of bodies is extremely large, one should instead solve the Vlasov equation-ie, the N → ∞ limit of the N-body model-coupled to the Poisson equation for the gravitational field. This is a formidable computational task, both in terms of run duration and memory storage, particularly for situations where intricate phase space structures develop over time, which is generally the case for cosmological simulations. Indeed, in order to solve the Vlasov equation, one needs to mesh the entire phase space, which is 6D in the most general case.
Fortunately, today's computer resources make it possible to envisage direct numerical simulations of the Vlasov-Poisson equations on a phase-space grid, if not for the full 3D case, at least in one spatial dimension. Compared to N-body codes, Vlasov codes exhibit a much lower level of numerical noise, particularly in regions of low matter density. Hence, they should constitute a valuable tool to study structure formation and in particular to understand the role played by large regions of underdensity in the Universe (so-called cosmic voids). A detailed comparison of the pros and cons of N-body (particle-based) and Vlasov (grid-based) codes is provided in appendix B. The present section is devoted to the presentation of cosmological results obtained with a 1D Vlasov approach, focussing on the Einstein-de Sitter Universe [120].
From the equation of motion (125) with d = 1, one defines a phase spacex,v, wherev ≡ dx/dt, and a probability distribution function F(x,v,t). The latter obeys the following Vlasov equation, with periodic boundary conditions in the scaled spatial coordinatex: which corresponds to the equations of motion (122) with Ω M = Ω T = 1 and Ω Λ = 0. 8 The distribution function is pushed in time using a split-operator scheme [136]. Interpolations on the grid are performed using an accurate finite-volume algorithm [137] that preserves the positivity of the distribution function. For the results presented below, we used 32 000 points in the spatial coordinate and 1000 points in velocity space. The initial condition is a cold Maxwellian, with a small variance in velocity space. The initial matter density ρ(x) is so constructed as to display a power-law spectrum of the type: |ρ k | 2 ∼ k 3 , where k is the wave number inx-space. This spectrum is compatible with that of an expanding Universe following an early inflation stage [20].  Figure 23 shows the phase-space distribution function at a later time. The distribution function displays a hierarchical structure at different scales, with small clusters orbiting each other to form larger clusters, which in turn also revolve around each other. This hierarchy is at the basis of the fractal structure also observed with N-body simulations.
Next, we show the wave number spectrum of the matter density |ρ k | 2 (figure 24). Rather quickly, a decreasing power-law spectrum builds up, with a slope approaching −0.5, which is close to the value observed for N-body simulations [116]. The range of the power-law region (k min , k max ) increases with time, with k min getting smaller and smaller while k max remains roughly constant. The steep cutoff at k > k max is due to numerical diffusion.
The clustering of the phase-space (figure 23) and the power law observed in the density spectrum (figure 24) point to an underlying fractal structure of the matter distribution, as was the case for the N-body simulations [118]. Box counting is the most standard method to analyze the properties of a fractal structure and here we use it to determine the generalized fractal dimension D q [128]. The approach is similar to the above but with the replacement of the number of particles in a given cell by the included mass, m i . As above, the system (0 ⩽x ⩽ L) is covered with boxes of length ℓ of decreasing size: ℓ = L/2, ℓ = L/4, and so on. Then the quantity I(q, ℓ) = i m q i (ℓ) is computed for each value of ℓ, where m i =´( i+1)ℓ iℓ ρ(ξ)dξ/m tot is the proportion of mass contained in the ith box, m tot is the total mass, and q is an integer which is intended to give more weight to high-density (when q > 0) or low-density regions (q < 0). As before, the fractal dimension is defined as ( where ω J is the Jeans frequency at θ = 0. Reprinted with permission from [120]. Copyright (2016) by the American Physical Society.
To improve the statistics, the value I(q, ℓ) was averaged over 1024 realizations. It can be proven [138] that D q should be a monotonically decreasing function of the exponent q. However, N-body simulations showed that, while D q displays the expected (decreasing) trend for q > 0, it is an increasing function for q < 0 (see open circles in figure 25). Since negative values of q over-represent low-density regions, this behavior was attributed to the poor sampling of these regions, where the number of particles is small and the statistics is thus very noisy [118].
Vlasov codes, by sampling the entire phase space with the same accuracy irrespective of the matter content, should provide better results precisely in such low-density regions. This is indeed what we observe in figure 25. For positive values of q, which are dominated by largedensity regions, the Vlasov and N-body results coincide; in contrast, for negative q the Vlasov curve (open squares) is basically flat (within the accuracy of the simulations) and definitely greater than the corresponding N-body results.
In order to compare to the N-body result, we have introduced an artificial cutoff in the density ρ obtained from the Vlasov simulations. Values of ρ that are below a certain threshold are set to zero. In figure 25, we show the results for four values of the threshold, ρ th = 10 −6 , 10 −4 , 10 −2 , and 1. It is clear that, by increasing the cutoff level, the Vlasov results progressively move towards the N-body results.
These findings nicely prove that the incorrect behavior of D q observed in N-body simulations obtained from the box-counting approach was indeed due to poor sampling of the low-density regions. In summary, the Vlasov approach confirmed the N-body results for highdensity regions, and extended such results to regions of low matter density. Therefore, in addition to employing mass-oriented methods for computing multifractal properties, Vlasov codes should be recommended for future studies focussing on cosmic voids, where the matter density is extremely thin.

Structure formation in the ΛCDM and Dirac-Milne cosmologies
The Einstein-de Sitter model discussed above (sections 3.2 and 3.5) predicts a Universe whose expansion is slowing down. Indeed, as the scale factor behaves as a(t) ∝ t 2/3 , the expansion rate goes asȧ(t) ∝ t −1/3 , and thus decreases with time. However, since the late 1990, several observations point to the fact that the Universe expansion is actually accelerating, or at least not decelerating. As a consequence, the Einstein-de Sitter model was progressively abandoned in favor of other cosmologies that account for such accelerating expansion.
To date, the generally accepted standard model of cosmology, known as ΛCDM, comprises an unusual mix of baryonic matter (ordinary nuclear matter) constituting less than 5% of the total mass-energy content of the Universe, cold dark matter (CDM, ≈ 25%) and dark energy (in the form of a cosmological constant Λ, ≈ 70%). The Friedmann equations corresponding to different cosmological models have been described in section 3.2, see equation (122). It is interesting that the various models differ only in the fictitious friction term present in the respective Friedmann equations, a term that results from the different expansion rates. In the case of ΛCDM, the rate of expansion goes through different phases: a rapid exponential acceleration in the very early Universe, followed by a slowing down in the matter-dominated-epoch, then a further exponential acceleration when the cosmological constant becomes dominant. The latter transition occurs roughly at the present epoch.
Given that the standard ΛCDM model features several unobserved 'dark' components, some cosmologists have been looking for alternative descriptions of the Universe history. An important family of such models are the so-called 'coasting' cosmologies [139], for Table 2. Interactions between matter (+) and antimatter (−) particles in the Dirac-Milne universe.

Type of matter
Type of matter Interaction which the scale factor grows linearly with time a(t) ∝ t, hence displaying no acceleration nor deceleration. The first of such coasting models, a Universe void of matter, was proposed by Milne [140] long ago. More recently, Benoit-Lévy and Chardin [23] proposed an alternative matter-antimatter symmetric Universe, which they named the 'Dirac-Milne' Universe, where the gravitational interaction between matter and antimatter, and between antimatter masses themselves, is repulsive (see table 2). This Universe, analogous in its gravitational behavior to the Dirac electron-hole system, avoids annihilation between matter and antimatter domains after cosmological recombination. The effect of gravity on antimatter is currently being investigated with several experiments developed at CERN which measure the acceleration of antihydrogen atoms free-falling in the gravitational field of the Earth. The first results of the Gbar [141], ALPHA-g [142], and AEgIS [143] collaborations are expected within the next few years.
Here, we summarize some recent results that compare gravitational structure formation in the Dirac-Milne Universe with that occurring in the Einstein-de Sitter and ΛCDM cosmologies [114,144]. The results were obtained with an N-body code that solves the equations of motion (122) for the various cosmologies. Typically, the simulations were performed with N = 2.5 × 10 5 particles.

Dirac-Milne vs. Einstein-de Sitter.
In figure 26, we show the 1D power spectra of the matter density for the EdS and Dirac-Milne cosmologies. At t = 0, both spectra follow a power law |ρ k | 2 ∼ k p with p = 2. The initial evolution is clearly linear, with each mode growing independently as a result of the Jeans instability. As time increases, nonlinearities become dominant and the spectrum takes a similar shape in both the Dirac-Milne and EdS cases. At low wave numbers, the spectrum is still of the type |ρ k | 2 ∼ k 2 , a remnant of the initial condition; at intermediate wavenumbers, the spectra follow a power-law spectrum with negative exponent, slightly steeper for Dirac-Milne (p = −0.78) compared to EdS (p = −0.67). This power-law region is a signature of hierarchical clustering in the phase space, as was observed in earlier works [108]. For even larger wavenumbers, the spectra are again flat, as one reaches the limit of resolution compatible with the number of particles used.
The spectra display a peak at wavenumber k peak , which represents the inverse of the typical cluster size. The position of this peak against time is shown in figure 27. In both cases, the peak initially moves towards smaller and smaller wavenumbers, somewhat faster for the Dirac-Milne case. For longer times (>10 9 y) the behaviors diverge: for EdS, k peak (t) keeps moving to larger scales, eventually reaching the size of the computational box, while for Dirac-Milne it saturates at a constant value. At the present epoch, the typical cluster size in comoving coordinates (k −1 peak ) is almost two orders of magnitude larger for the Dirac-Milne Universe compared to the EdS case. We also stress that such cluster size is determined primarily by nonlinear effects occurring during the matter-dominated epochs. 3.6.2. Dirac-Milne vs. ΛCDM. Finally, using the same N-body code, we performed numerical simulations comparing structure formation in the Dirac-Milne and ΛCDM cosmologies [114], assuming for the latter Ω M = 0.3 and Ω Λ = 0.7. Here, distances (in the comoving frame) are measured in physical units, and time is expressed in terms of the cosmological redshift z = (1 − a)/a, where a(t) is the scale factor. We set the initial condition at recombination (z = 1080), but note that this redshift does not correspond to the same epoch in the Dirac-Milne (≈14 million years) and ΛCDM (≈380 000 years) cosmologies, see [144]. As above, the initial matter density fluctuations follow a power law spectrum P 1D (k) = |ρ k | 2 ∼ k p , with p = 3. The 3D spectrum P(k) then satisfies: P 1D (k)dk = P(k) 2πk 2 dk, yielding P(k) = P 1D /(2πk 2 ), and hence behaves as the standard Harrison-Zeldovich spectrum P(k) ∼ k.
The power spectra at z = 0 for the Dirac-Milne and ΛCDM cases are displayed in figure 28. The horizontal axis was scaled so the peaks are located at 0.018h Mpc −1 , where h = H 0 /(100 km s −1 Mpc −1 ), as in the observed SDSS spectrum [145]. 9 The shape of the two spectra are virtually identical, pointing to the fact that, although the Dirac-Milne and ΛCDM   [114]. Copyright (2020) by the American Physical Society. universes go through very different evolution stages, the outcomes at z = 0 are very similar (and compatible with observations). In addition, the expected slopes for long wavelengths (P ∼ k) and short wavelengths (P ∼ k −3 ) are recovered by the Dirac-Milne simulations. We also note that this power law behavior, which is mainly due to nonlinear effects after recombination, suggests a self-similar matter distribution and the existence of a robust fractal dimension [108,119].
The evolution of the peak wave number of the power spectrum (figure 29) shows that structure formation is initially faster in the Dirac-Milne Universe. However, for scale factors a ≈ 0.25 or larger (cosmological redshift z ≈ 3, or ≈ 3.5 billion years after the Big Bang), k peak saturates around a value corresponding to ≈ 50 h −1 Mpc. This indicates that structure formation has stopped at a similar epoch for both universes.
Overall, the above results show that precise cosmological estimations can be obtained using 1D models. Comparisons between the ΛCDM and Einstein-de Sitter cosmologies show that Figure 29. Wave number k peak corresponding to the peak of the power spectra for the Dirac-Milne and ΛCDM universes as a function of the scale factor a(t). A striking feature of these simulations is that k peak evolves in time, describing the nonlinear evolution in both the Dirac-Milne and ΛCDM models, whereas in the standard ΛCDM analysis the usual assumption is that the nonlinear evolution represents only a small correction, while the peak position is fixed at k peak ∼ 0.018 hMpc −1 corresponding to the mode entering the horizon at matter-radiation equality. Reprinted with permission from [114]. Copyright (2020) by the American Physical Society. the latter is ruled out by current data on the large-scale structure of the Universe. Alternatives, such as the Dirac-Milne cosmology, fare better as long as structure formation is concerned, although more studies are needed to evaluate their validity against other cosmological data, particularly the CMB spectrum.

Concluding discussion
In this section, we discussed the use of 1D methods to study various cosmological scenarios. Although obviously not realistic, as our Universe possesses three spatial dimensions (not to mention the 10 dimensions of string theory), 1D models present some distinct advantages for cosmological applications.
First, the Newtonian interaction is very simple in 1D: it displays no divergence at short range (hence no need for an artificial cutoff as in 3D) and is constant in space at long range (actually, at all distances), again needing no long-range cutoff. This peculiar form of the interaction allows one to integrate exactly the N equations of motion of the particles (which, in 1D, are in fact infinite plane sheets) in between two subsequent crossings between the particles themselves. For the RF model the crossing times can be computed analytically: hence, by moving from one crossing event to the next, one can solve exactly the N-body dynamics, the only source of error originating from the finite number of decimals available to code real numbers in the computer. For the Q model, it is necessary to solve for the crossing times numerically, but the overall scheme going from one crossing event to the next remains the same. The resulting exact N-body code is described in some details in the appendix B. This type of approach is distinct from the more usual molecular dynamics approach, where the equations of motion are solved in an approximate fashion using some time-stepping techniques (Runge-Kutta, Verlet, . . .) that are exact only up to a certain order in the time step ∆t.
Further, the low dimensionality of the models is also advantageous for the development of grid-based codes, which integrate directly the Vlasov equation on a phase-space mesh, see the appendix B. Although potentially more accurate, these codes imply a higher computational cost, which is prohibitive in 3D as the corresponding phase space is 6D. In this review, we presented results obtained with both types of codes. Grid-based codes were more costly, but particularly efficient to evince subtle structures in low-density regions, which are difficult to access for particle-based codes.
3D codes often suffer from poor sampling, because of the limited number of particles that can be used. Hence, whenever accurate statistical results are required, 1D models can be advantageous. For instance, this is the case when one wants to study the (multi)fractal nature of the mass distribution in the Universe. In order to extract the value of the fractal dimension from the code, it is necessary to have good statistical sampling over several decades, which is easily feasible in 1D, but much less so in 3D.
Some tension between the 3D physics and the 1D models arises when one deals with the Universe's expansion. Indeed, the expansion itself is intrinsically 3D (e.g. in spherical coordinates): to obtain a 1D model, one considers a thin shell around some radius R and then assimilates this shell to a planar slab. However, this procedure-which amounts to considering 1D perturbations embedded in a 3D Universe-is not entirely consistent, and the resulting 1D Poisson equation has to be slightly modified ad hoc. This yields the so-called 'Quintic' (Q) model.
Alternatively, one can start from a fully 1D expanding Universe. In this case, the whole procedure is exact, but one has to assume that the surface mass density of each sheet decreases as 1/a 2 in order to mimic the 3D expansion, where a is the expansion factor. This yields the so-called Rouet-Feix (RF) model. Historically, the RF model was the first 1D cosmological model considered. In practice, the Q and RF models are very similar, and differ only in one of the coefficients in the equations of motion of the sheets.
Cosmological simulations are usually performed in comoving coordinates, i.e. coordinates that follow the Universe expansion with the time-dependent scale factor a(t). In the comoving coordinates, the Universe is static: however, due to the Jeans instability, this static state is unstable and after some time it develops intricate hierarchical structures in the form of galaxies, clusters of galaxies, and superclusters. As stated above, 1D codes enable us to resolve these complex gravitational structures with such precision that their statistical properties can be thoroughly investigated, in particular their fractal or multifractal dimension and the evolution of the power spectrum of the mass density P(k).
Here, we showed some examples of such results obtained for several 'universes' characterized by different scale factors a(t). Most simulations were performed within the Einstein-de Sitter Universe, which is a critical case where the expansion is self-similar and a(t) ∝ t 3/2 . While the exact N-body code is useful to evince gravitational structures in regions of large density (where many particles are present and thus the statistics is good), the grid-based code provided yet-unattained insights into the behavior of dilute regions such as cosmological voids.
We also discussed structure formation for two other cosmological models: the ΛCDM model, which includes a cosmological constant and is the standard scenario of modern cosmology; and the so-called Dirac-Milne Universe, which is an alternative scenario where matter and antimatter coexist in equal amounts, and the latter possesses negative gravitational mass, leading to a linear scale factor a(t) ∝ t. Both universes display very similar features as regards to cosmological structure formation.
Finally we mention two other recent approaches for investigating gravitational collapse: Firstly, the study of singularities with a Lagrangian-based dynamical approach, which originates from the use of Burgers' equation for the study of turbulence in fluid dynamics. Recently, there have been extensions of the Zeldovich approximation to investigate the development of secondary singularities in the mass flow [146]. This an analytic approach to studying the mechanism for clustering. Secondly, we mention the recent application of wavelet analysis to the dynamics of gravitational clustering, both in 1D using the Zeldovich approximation and in 3D employing large scale N-body codes [147]. It will be interesting to see what new insights develop from these methods.
In summary, 1D gravitational models have proved very useful to explore the major outstanding questions of today's cosmology and constitute a valuable tool alongside more standard 3D simulations.

Conclusion
We hope that the reader is convinced by now that one-dimensional models of Newtonian gravity provide a useful laboratory for investigating concepts of current interest. The observation that stellar motion perpendicular to the plane of a highly flattened galaxy could be modeled by 1D gravity provided some initial motivation. Starting in the 1960s, the ability to precisely simulate the evolution of these toy systems has distinguished them from their more realistic 3D counterparts. Because of the long range of the gravitational force, these systems do not admit to standard theories of statistical mechanics, both in 1D and 3D. Their investigation requires non-extensive thermodynamics and reveals ensemble inequivalence. Moreover, in the infinite particle limit, these systems obey the Vlasov (or collisionless Boltzmann) equation, so their relaxation time grows indefinitely with N, the number of particles. We have shown here that this dependence is still open to question after all these years. About forty years ago investigations by two groups (see [39,148]) showed that relaxation depended preferentially on the initial condition and, in spite of intense investigation, this picture still remains today. As a result of filamentation in µ-space (ie, the single-particle phase space), relaxation of the Vlasov equation would require some type of coarse graining or, alternatively, a time dependent spectral analysis, and recent work has shown progress in this direction.
In order to understand the nature and source of ergodic behavior in these systems, investigators turned to small systems which were amenable to nonlinear analysis. The equivalence of a wedge-like geometry with the three-particle system led to laboratory experiments using optical lattices and inelastic collisions. It was shown that the µ-space consists of both stable and unstable regions identified with periodic orbits and the unstable regions dominate with increasing N. The study of Lyapunov spectra provided great insight into the relaxation process. Analogous studies were undertaken on relativistic versions of the three and four particle systems. The mathematical structure of the phase space is still an open question. Recently the role of dissipation in the relaxation process has also been investigated.
The hypothesis that the evolution of the Universe after recombination proceeded by the early formation of Zeldovich pancakes surely stimulated interest in 1D models of cosmology. This arose from the failure of the more commonly invoked spherical collapse model to explain the size of structures in the present Universe [21,149]. Nonetheless much insight has been gained from the spherical collapse model and it also continues to play a role in exploring different cosmologies [150,151].
By appropriately scaling both position and time, it was shown in the 1980s that it was possible to formulate a workable version of the 1D system in a comoving coordinate frame incorporating the expansion of the Universe. The ability to precisely explore the system evolution made it possible to capture the underlying fractal nature revealed in the simulations as time evolved. The investigations were refined by improved algorithms and the addition of periodic boundary conditions. The multi-fractal aspects of the evolution have only begun to be explored. The system was later extended to include two components-one dissipative representing baryonic matter and one conservative representing dark matter. It was shown that the fractal properties of each component differ. Recently the 1D system was used to compare different cosmological scenarios. In particular, a recent comparison of the standard ΛCDM and less conventional Dirac-Milne cosmologies demonstrated that both exhibit 'freeze out' in the clustering process on similar time scales. Work is continuing in this arena. The application of wavelet analysis and the evolving application of Lagrangian dynamics in the configuration space may provide additional insight in the future. In conclusion, it is unlikely that the original investigators of the 1D gravitational system could have anticipated all of these applications and developments. Clearly the story still has not been completed.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.
Finally, δ is obtained from the normalization condition This equation gives (19).

A.2. Proof of proposition 2
The proof follows closely the derivation of thermodynamical properties of a one-dimensional gas with long-range interaction in [19], see section 3 in this reference. The configurational contribution to the free energy is the extremal value of the functional f Q [ρ] given in (17). It is convenient to put it into a simpler form as where J =´L /2 −L/2 ρ(x) log ρ(x)dx, and λ is the Lagrange multiplier introduced in (A.1). The Lagrange multiplier λ has a physical interpretation as the chemical potential of the gas. They are given explicitly by The density function h(x) = log ρ(x) and its derivative at the boundary are given by the solution of the Lane-Emden equation (18) h We start by presenting the derivation of (A.7). This is obtained by multiplying the Lane-Emden equation (A.2) with 1 2 ρ(x) and integrating over x. This expresses the interaction energy as a one-dimensional integral Substituting into f Q gives the result (A.7).

B.1. Exact N-body code
In one spatial dimension (1D) and planar geometry, the dynamics of N bodies interacting through Coulomb or Newton forces can be solved exactly. Indeed, in 1D the 'particles' are actually infinite planes of uniform surface mass density. In that case, the field created by a particle is −sign(x − x i ) 2πGµ i , where x i the position of the particle i and µ i its surface mass density. Unlike in the 3D case, the field has no divergence at the particle position x i (only a discontinuity), so that there is no need to regularize (smooth) the field at short range, and the particles are allowed to cross each other (see figure B1). For a collection of N particles, the field perceived by any one particle is piece-wise constant in space and only depends on the number of particles situated to its right and to its left. Thus, as long as the particle does not cross one of its neighbors, the field is constant in time and its trajectory can be computed exactly.
More precisely, the equation of motion of particle i is given by equation (125): with the one-dimensional field N + (x) and N − (x) being the number of particles located respectively on the right and on the left of x i . As these numbers remain constant, so does the fieldÊ x , which consequently has a step-like shape. It is convenient to normalize the time to the Jeans frequency by introducinḡ t = ω j0t , and to take the mean distance between two particles equal to unity (so that N = L).
In these units, we have to solve the following equation of motion (overcarets are omitted for simplicity of notation): In between two crossings, the field E i experienced by particle i is constant and the solution of equation (B.1) reads: x i (t) = A i exp(r 1 (t − t i0 )) + B i exp(r 2 (t − t i0 )) − E i v i (t) = A i r 1 exp(r 1 (t − t i0 )) + B i r 2 exp(r 2 (t − t i0 )), (B.4) where t i0 is some initial time for which x i (t i0 ) = x i0 and v i (t i0 ) = v i0 , r 1  For each pair of particles that see a new neighbor (because one of their neighbors experienced a crossing), or at initial time for all pairs of neighbors, we have to compute their crossing time t c for say particle i and its right neighbor j = i + 1. At t = t c , we have which leads to For the RF model, equation (B.11) yields a third-order algebraic equation, which reduces to a second order equation after particles i and j have just collided (at t i0 = t j0 , x i0 = x j0 , so t c = t i0 is an obvious solution). For the Q model, equation (B.11) yields a fifth order equation. For the simple gravitational system (without friction nor background) crossing times are determined by the resolution of a second order equation.
We note that, as the 1D gravitational field is constant and does not decrease to infinity as the inverse square of the distance, long-range collective effects are even more important in 1D than they are in 3D. Thanks to these special features of the 1D Coulomb and Newton forces, an exact algorithm can be set up, which works according to the following steps: • Compute the crossing times between all neighboring pairs of particles; • Select the smallest crossing time and the corresponding particle pair; • Advance these two particles until the computed crossing time; • Recalculate the field seen by the particles and compute their velocities at the last crossing time; • Go back to initial step.
Note that only the field seen by the two particles that cross each other is to be recalculated (the other particles keep the same number of particles to their right and to their left, so that the total fields they perceive are unchanged) and only three new crossing times are to be recalculated at each crossing, i.e. those of the two particles that have just crossed and their direct neighbors on the right and on the left. The information needed to compute the successive crossings is the initial condition for each particle after a crossing, ie the position, velocity, and field felt by each particle at the last crossing time.
Finally, if the system is confined in a box (with either periodic or hard edges), the box edges may be modeled as fictitious particles. In the case of a periodic box, an Ewald sum which takes into account the field of the replicas could be computed [93]. As already mentioned, periodic boundary conditions are a common choice to mimic infinite systems, but they act as a lowpass filter and erase wavelengths longer than the system size. The field at a given point takes into account all the particles in the box and all its replicas. The non-decreasing nature of the 1D Newton field prevents the sum of all these replica contributions to converge. A way to obtain the convergence was suggested by Kiessling [10], who proposed to screen the field by an exponential factor whose length scale is set to infinity at the end of the calculation.
The same purpose may be achieved by polarising the boundaries: as soon as a particle leaves the system on one boundary, this side receives the (gravitational) charge of this particle, while another particle enters from the opposite side and leaves an opposite charge on that boundary. Hence, the particles do not experience any change of their respective fields, even though the number of particles at the right and the left have changed.
In this type of algorithm, the particles advance two by two at the rhythm of the crossing times. If the system is frozen at any one instant in this process, the particles are located at different times that correspond to their last crossing time with their neighbor. In order to store the state of the system at a given moment, one just needs to advance the set of all particles until that time, which is necessarily between two crossings.
With the above algorithm, the trajectories are known exactly (or, rather, with the precision allowed by the number of bits with which are coded floating numbers), crossing-time after crossing-time. Double or quadruple precision is required to avoid the risk of miscalculating a crossing time, which would result in the particles being no longer ordered according to their position. This is a crucial criterion that must be respected by this type of code.
The heart of the algorithm is to keep the order relation between the positions of the particles, so as to search for the smallest crossing time over all the N − 1 crossings of the N particles with their neighbors. However, it is not necessary to sort all the crossing times, as only the smallest one is relevant for the algorithm. To this purpose, two different techniques have been proposed in the past, namely the fast table and the heap-sorting technique. on. At the initial time of the simulation, the table is empty and all the N − 1 crossing times have to be computed. If the crossing time between the particles i and i + 1 belongs to the duration covered by the table, then the index i of the particle is placed in the box corresponding to the time interval during which the crossing will take place ( figure B2 gives an example of a system of 8 particles).
Once this procedure has been performed for the N − 1 crossing times, one just needs to check the first box of the table: if it is empty, there is no crossing between T and T + ∆T, then we look at the next box, and so on. If one box contains a particle index, then the corresponding crossing is processed, three new crossing times are calculated, and the indexes of the respective particles are assigned to one of the cells of the table (if the corresponding crossings take place during the interval of time covered by the table). Beforehand, it is advisable to erase from the table, if they are there, the indexes of the three particles whose crossing times have become obsolete. When the last box in the table has been processed, the table is empty again and ready to process the events occurring over the next M∆T duration.
If the table contains enough boxes, each box is the seat of only one particle index at a time. Otherwise, one needs to shift the indexes of the particle that have the largest crossing time to the right of the table, so as to have only one index per box. If a box is already occupied with an index, it also has to be shifted to the right until an empty box is reached. Alternatively, when a crossing is treated, it is necessary to check if backward shifts have to be performed. These shifts slow down the sorting process and should be avoided by having a sufficiently large table, i.e. boxes whose time interval is sufficiently short. A table that is 75% empty seems to be an optimum (an even higher rate would result in exploring an almost empty table, thus unnecessarily occupying memory space). The algorithm scales as O(N) for the filling of the table and as O(1) for the search of the smallest crossing time.

B.1.2. Heap-sorting technique.
For sorting using the heap-sorting method [152], a tree with N − 1 nodes is built from top to bottom, each node having two (or more) branches. The nodes contain the N − 1 crossing times of the N particles with their neighbors. Initially, the tree is empty and is filled from above (ie, from its root), crossing time after crossing time. If the crossing time of a child node is smaller than that of the parent node, then the values are interchanged, and the smaller crossing time climbs up one level in the tree. At the end of the fill, the smallest crossing time is at the top of the tree. At each crossing, one must readjust the positions of the three recalculated crossing times. In order to allow the updating of the tree, this algorithm requires two complementary tables to associate a node with a particle index, and vice-versa. The complexity of the algorithm (once the tree has been constructed) is of order O(log N) per crossing.
An alternative technique of the same complexity consists in starting from the bottom of a tree comprising (N − 1)/2 'leaves' (a leaf is a node with no ramifications to the bottom). Initially, as the tree is empty, one compares two by two the crossing times of consecutive particles (for example crossings between particles (i, i + 1) and particles (i + 1, i + 2)), and the label of the particle which experience the smallest crossing with its right neighbor is record at the bottom of the tree. As i goes from 1 to N, there are N − 1 crossings, so that (N − 1)/2 leaves are necessary. The process is repeated to fill the next level of the tree and so on to reach the top of the tree, which will contain the index of the particles that have the smallest crossing time with their right neighbor. Taking N = 2 k + 1, the tree has k levels for N − 1 nodes. At each new crossing, there will be 2 k new comparisons at most, and k at the least. It is possible to reduce this number by comparing not only two but rather three or more crossings at each level. An equivalent question has been raised for the heap sorting method, and three appears to be an optimum [152].

B.2 Molecular dynamics and particle-in-cell (PIC) codes.
The methods discussed above can be applied whenever it is possible to integrate exactly the equations of motion. When this is not the case, they can be integrated numerically on a time grid with time step h. If two particles cross each other, it is necessary to recalculate the field. The numerical error depends on the order of the time-stepping integrator (usually, the Störmer-Verlet method of order 4 or the velocity-Verlet method of order 3). This technique is similar to the 'molecular dynamics' methods currently employed in other areas (condensed matter, biophysics . . .). It is less precise than the exact algorithm described above, but avoids the cumbersome computation and ranking of the crossing times.
Finally, it is worth mentioning the particle-in-cell (PIC) method, whereby the gravitational field E(x, t) is computed on a fixed mesh using Poisson's equation: ∂ x E = −4πGρ, where ρ is the mass density. The N particles obey the standard equations of motionẍ i = E(x i , t). A typical PIC code operates along the following steps: • From the discrete distribution of the N particles, compute the smooth mass density ρ, using some regularizing kernel; • Using Poisson's equation, compute the field E(x, t); • Using the computed field and a time stepping algorithm, integrate the equations of motions from t to t + h; • Update the particle positions x i (t + h); • Go back to initial step.
It must be noted that the PIC algorithm does not solve the full N-body problem (as does the exact code described above), but rather a mean-field approximation to it. This should be evident from the fact that the field employed in the equations of motion is not the exact (piecewise constant) field, but a smoothed approximation to it. In this approximation, two-body and higher-order correlations are neglected and only the mean (self-consistent) field is retained.

B.3 Mesh-based (Vlasov) codes.
Numerical simulations of self-gravitating systems are generally based on N-body codes, such as those described above in B.1 and B.1.3, which solve the equations of motion of a large number of interacting particles. However, this approach may suffer from poor statistical sampling in regions of low density, where the distribution of particles is insufficient to accurately sample the mass density. Instead, Vlasov codes work by covering the entire phase space (x, v) with a uniformly spaced grid and evolving in time the particle probability distribution function f(x, v, t). They are able to provide an accurate description of the distribution function in all regions of the phase space, but their computational cost is higher, because the entire phase space needs to be meshed.
Vlasov codes work by solving numerically the Vlasov equation: where the gravitational field E obeys the equation: and ρ =´fdv is the mass density. Just like the PIC codes described in the B.1.3, Vlasov codes also solve a mean-field approximation of the exact N-body problem. The time-stepping technique used to solve equation (B.12) and depicted in figure B3 is based on a splitting algorithm [136], which treats separately the free-streaming term ∂f ∂t + v ∂f ∂x = 0, (B.14) and the acceleration term ∂f ∂t + E ∂f ∂v = 0, (B.15) in the Vlasov equation (B.12). The solution from time t n to time t n+1 can thus be obtained in three steps, corresponding to the solution of the free-streaming term (B.14) over half a time step, then the solution of the acceleration term (B.15) over a full time step, and finally again the free-streaming term (B.14) over half a time step: where f * and f * * denote intermediate solutions. The gravitational field E is computed from Poisson's equation (B.13) just before equation (B.17). Using the above symmetric scheme, the method is second-order accurate in the time step ∆t. We note that each term in equations (B.16)-(B.18) gives rise to a constant shift in either position or velocity space. In their numerical implementation, these shifts require the interpolation of the distribution function in phase space, which can be performed by means of different schemes (cubic splines, finite volumes, fast Fourier transforms . . .). For the results reported in section 3.5, we employed a numerical technique based on a finite-volume algorithm, in which the mass distribution is assimilated to a phase-space 'fluid' [137]. The scheme performs a detailed balance of the fluid entering and leaving each phase space cell: in this way, the total mass´´fdxdv is conserved exactly. Using a slope limiter, the finite-volume method can be made to avoid spurious negative values in the distribution function.
Recent simulations of self-gravitating systems in more than 1D using Vlasov codes have employed adaptive meshes [153] and higher-order time stepping techniques [154].