1. Introduction
Understanding the mass distribution in self-gravitating systems has remained a longstanding challenge that has intrigued researchers for over seven decades [
1,
2,
3]. The resolution of this enigma holds the potential to illuminate many theoretical puzzles, such as the physical mechanisms underlying the regularities observed in the light profile of elliptical galaxies and the mass distribution in dark matter halos.
Classical simulations conducted by Navarro, Frenk, and White [
4,
5] produced density profiles of dark matter halos that current theories struggle to explain. Numerous references [
6,
7,
8,
9,
10,
11,
12,
13] explore the scope of this problem, extending from the foundations of statistical mechanics to the large-scale evolution of the universe [
14].
In contrast to systems with short-range interactions, which equilibrate through collisional processes, long-range systems relax primarily through collective effects arising from the interactions of individual particles with the entire system. The long-range interacting systems exhibit intriguing features such as: ensemble inequivalence [
15], temperature inversion effects [
16], and the emergence of long-lived quasi-stationary states [
17]. These unique characteristics set them apart from their short-range counterparts and make them fascinating subjects for study and exploration.
In the thermodynamic limit, interparticle correlations are absent in systems with long-range interactions, and the
N-particle distribution function,
, can be factorized into a product of one-particle distribution functions,
,
where
are the generalized coordinate and momentum of a particle in the phase space. This factorization reduces the
-dimensional phase space to the phase space of a single particle. The evolution of the one-particle distribution function is governed by the Vlasov equation [
18]:
where
denotes the mean-field potential associated with the pair potential
.
The evolution of
in the phase space is similar to that of an incompressible fluid—maintaining a constant local density along the flow. Furthermore, the Vlasov dynamics have an infinite number of invariants known as the Casimirs [
19,
20,
21,
22]—any local functional of the distribution function is a Casimir invariant. In particular, a hyper volume of any density level
of the initial distribution function is a Casimir invariant of the Vlasov dynamics,
The total mass
M and the energy of the system
are also conserved quantities:
In the case of non-neutral plasmas,
M would refer to the total charge of the system [
23,
24]. Finally, both linear and angular momentum are also conserved quantities. However, for symmetric initial distributions studied in the present paper, their initial values can be set to zero.
Unlike the collisional Boltzmann equation, which has the Maxwell distribution as the global attractor, the Vlasov equation does not have an attractor. The relaxation to the stationary state of the Vlasov equation proceeds through the process of filamentation. In fact, on a fine-grained scale the relaxation never stops—only on the coarse-grained scale can we say that the system has reached a stationary state. Any distribution of the form
, where
is the single-particle energy, is a potential candidate to describe the stationary state (
). However, predicting
a priori the specific distribution to which a system will relax starting from an arbitrary initial condition remains a challenge [
1].
Long before the empirical work of Navarro, Frenk, and White [
4,
5], Lynden-Bell (LB) introduced a very elegant statistical theory of relaxation in systems with long-range interactions, with particular emphasis on self-gravitating systems [
7]. This theory became known as the Theory of Violent Relaxation. By discretizing the initial particle distribution into multiple levels within the context of Vlasov dynamics, LB assumed that the evolution of the density levels is ergodic and that there is a complete mixing of the one-particle distribution function in the phase space. Under these assumptions, LB argued that on a coarse-grained scale, the system will evolve to the maximum entropy state allowed by the conservation of mass, energy, and the Casimir invariants. In particular, for an initial distribution of one level—the waterbag distribution—the coarse-grained stationary distribution function is found to be,
where
and
are the Lagrange multipliers determined by the conserved quantities (Equations (
4) and (
5)), and the potential is calculated self-consistently by solving the Poisson equation,
where
is a constant associated with the dimensionality of the system and
represents the Laplacian operator in the
d-dimensional space [
25].
The form of the distribution proposed by Lynden-Bell (LB) (
6) exhibits similarities with the Fermi–Dirac distribution. However, it is essential to emphasize that we are dealing with classical “particles”. The exclusion principle associated with fermions, in this case, has a parallel with the constraint that the evolution of the distribution function, according to Vlasov dynamics, is analogous to that of an incompressible fluid, meaning that different levels of the distribution function cannot overlap at any point in the phase space.
At the core of LB’s theory is the fundamental assumption that violent relaxation induces a rapid phase-space mixing of the levels of the initial distribution function. On a coarse-grained scale, the one-particle distribution function is expressed as an average over the density levels present inside the macrocell, as depicted in
Figure 1.
We stress that the equilibrium distribution function proposed by LB exists only at the level of macrocells, as there is no increase of entropy at the microcell level—fine-grained Gibbs entropy is one of the Casimir invariants of the Vlasov dynamics. This is similar to what occurs in short-range systems when analyzed in terms of the
N-particle distribution function in a
-dimensional phase space. The increase of Gibbs entropy takes place only at the coarse-grained macrocell level since at the microcell level, the Gibbs entropy remains constant due to Liouville’s theorem. In both cases, this problem can be circumvented by using Boltzmann entropy, which is defined at the coarse-grained level [
26,
27]. Nevertheless, it is worth emphasizing that LB’s statistical approach is equivalent to the traditional method of maximizing the Gibbs entropy, provided that all the Casimirs are taken into account [
22].
When attempting to apply Lynden-Bell’s statistical theory to real-world systems that involve multiple levels
l of the initial one-particle distribution functions, a significant challenge arises due to the highly nonlinear nature of the problem. In this case,
l equations, similar to Equation (
4), are required, one for each level, which substantially increases the computational complexity. Consequently, obtaining solutions under these general initial conditions becomes impractical for large
l, making it challenging to achieve accurate and reliable results. Moreover, the system’s sensitivity to initial parameters further compounds the difficulties, hindering the discovery of the final equilibrium state. As a result, the application of Lynden-Bell’s theory has been primarily confined to single-level distribution functions [
28]. While this restricted set of initial conditions offers some valuable insights into certain aspects of self-gravitating systems, its generality is limited—preventing us from making conclusions about what will happen with much more complex continuous initial distributions. To address these challenges, the primary objective of the present paper is to propose a Monte Carlo (MC) approach that efficiently leads to the equilibrium state corresponding to the maximum of LB entropy, with all the relevant constraints taken into account.
2. Monte Carlo Algorithm
We start by discretizing the initially continuous distribution
into
l levels [
29], where each level
i has a density
, with
. We designate level
as the empty cell, with
. Thus, the initial distribution can be expressed as:
Here,
represents the number of microcells containing level
i inside the macrocell
j, centered on
. The total number of macrocells is defined as
, where
is the number of position grids, and
is the number of momentum grids. At
, each macrocell centered on
has only one level
, with the value closest to that of
. Since each microcell can contain only one density level, we have a bound
, where
is the number of microcells in a macrocell, see
Figure 1.
The key aspect of the Monte Carlo (MC) technique is to evolve the system to an equilibrium state consistent with the macroscopic constraints, such as the total energy, mass, momentum, and the conservation of the total volume occupied by each density level. In equilibrium, the system can explore all microstates that are consistent with the constraints. To force the system towards equilibrium, MC moves must respect the detailed balance condition:
where (
o) and (
n) refer to old and new configurations corresponding to two different microstates.
is the equilibrium probability that the system is in a given microstate. As can be seen in
Figure 1 and
Figure 2, all microstates on the energy shell are equally probable, so that
for any microstate.
The transition probability
consists of two steps:
where
represents the proposal of a trial move from one microstate (
o) to another (
n), and
denotes the acceptance probability associated with deciding whether to accept or reject this trial move.
In the MC simulation, various trial moves are possible, one simple example being the exchange of levels between two microcells. In this case, it is important to note that
since we can always swap levels between any two arbitrary microcells (on the energy shell) (see
Figure 2). As a result, the acceptance probabilities for both transitions are equal:
meaning that a swap between one or more pairs of microcells is always allowed, provided it conserves the total energy and other macroscopic constraints.
Such a brute force approach is very inefficient since only microstates that are confined to the energy shell can be moved. Most moves will be rejected because they do not conserve the total energy of the system. To enhance the efficiency of such microcanonical MC methods, Creutz [
30] proposed to relax the constraint of exact energy conservation—allowing the total energy of the system to slightly fluctuate close to its initial value. This is achieved by introducing an additional degree of freedom for a system to store energy. Creutz called this additional degree of freedom a “demon”. The first move is only allowed if it lowers the energy of the system. The excess energy is then stored inside the demon. In the following move, the difference
is examined. If the trial move lowers the energy, it is accepted and the energy gained from such a move is stored inside the demon. On the other hand, if
, the move is accepted
only if the demon has sufficient energy in store to compensate for the system’s energy gain. If it does, the move is accepted, and the quantity
is subtracted from the demon’s energy store. The MC dynamics is illustrated in
Figure 2.
For application to LB statistics, a simple MC algorithm described above is still very inefficient. The position of each microcell has to be stored in an array and track must be kept of all the swaps between the different microcells. Clearly, since the energy and the coarse-grained distribution function are calculated on the level of macrocells, the details of which microcell inside the macrocell is occupied by the specific density levels are irrelevant. The only relevant information is how many density levels of each type are present inside the macrocell. Based on this observation, we now introduce a more efficient approach for performing MC at the level of macrocells.
As before, the phase space is divided into
macrocells, each containing a fixed number
of microcells that can either be occupied by a certain level
i of the initial distribution or left empty, as illustrated in
Figure 1. Trial moves now involve adding or removing levels within a specific macrocell, as demonstrated by the processes
and
of
Figure 2. Only process
is relevant, since, in process
, the proposed exchange is between equal levels and does not affect the coarse-grained distribution function.
To illustrate the algorithm, let us consider an exchange of two distinct levels, denoted
and
, between two different macrocells, labeled
and
, as shown in processes
of
Figure 2. In this trial move, we will remove a level
a from the macrocell
A and place it into macrocell
B, and similarly, we will remove a level
b from the macrocell
B and place it into macrocell
A. This is summarized as follows:
Note that the redistribution of levels within the same macrocell does not affect the energy of the system, since it is calculated only at the level of macrocells. The total degeneracy
W for a given distribution of levels
across all the macrocells is, therefore,
We should observe that differently from LB’s original work, we adopt Gibbs counting, and treat the same density levels as indistinguishable [
7]. In practice, this does not affect any of the final results, and makes the combinatorics simpler.
The probability of finding the system in its old configuration, characterized by the macrocell
A with
, and the macrocell
B with
, will be proportional to:
while after performing the swap, the probability of finding the system in its new configuration, characterized by the macrocell
A with
, and the macrocell
B with
, will be proportional to:
Note that the choice of density levels that are being proposed for a swap is completely unbiased, performed with equal probability among all the levels that are present within the two macrocells. Therefore, the probability of choosing a trial move
. Substituting this into the detailed balance equation, Equation (
9), we find that the acceptance probabilities in the forward and reverse directions must satisfy:
Since the maximum acceptance probability is one, following Metropolis [
31], we conclude that:
In this new approach, the trial moves within the same macrocell, depicted by process
I of
Figure 2, are no longer necessary since their mixing and degeneracy are taken into account exactly by the combinatorial factors in Equations (
15) and (
16). This represents a significant gain compared to the original MC method. To summarize: we randomly select two out of the
macrocells in the system, along with two levels from the available levels within these macrocells. Next, we generate a random number
, uniformly distributed between 0 and 1, and then check if
. If this condition is met, we accept the swap; otherwise, we reject it. Note, the choice of trial moves remains symmetric,
, as we continue to choose the levels for a swap independent of their occupation of the macrocell. Although this MC algorithm represents a significant improvement over the brute force MC, the fact that
implies that many of the proposed moves will be rejected. This unnecessary inefficiency arises when we encounter cases where
. In such situations, there are probably other levels (different from
a and
b) that are more appropriate for a swap. This suggests that the choice of proposed moves should not be completely random, but biased [
32] toward levels that have larger occupations within the macrocells
A and
B.
To implement this bias, we now select the levels based on their occupation in the macrocell, rather than choosing them uniformly as we did before. The probability of choosing level
a inside the macrocell
A is taken to be proportional to its occupation within the macrocell
, and similarly for level
b inside the macrocell
B. The biased choice then:
Note that now,
since the choice in the reversed direction is:
To make this discussion clearer, consider the central macrocell of
Figure 2, and let us call it
A. Now, let us associate colors to represent different levels: yellow for level 1, red for level 2, green for level 3, blue for level 4, pink for level 5, and empty spaces for level 6. Therefore, in macrocell
A, we have the following distribution:
,
,
,
,
, and
, totaling
.
To randomly select a level within macrocell
A, we generate a random number between 0 and
. If the generated number falls between 0 and 1, we choose level 5 (pink). If it falls between 1 and 3, the chosen level will be 3 (green). If it falls between 3 and 9, the chosen level will be 4 (blue). Finally, if the number falls between 9 and 25, the chosen level will be 6, representing an empty space. This is similar for macrocell
B. The probability of selecting a level
a in macrocell
A and a level
b in macrocell
B is then given by Equation (
19). Substituting Equations (
19) and (
20) together with Equations (
15) and (
16) into the expression for the detailed balance Equation (
9), we obtain:
Therefore, as long as the levels proposed for a swap are chosen according to their occupation inside the two macrocells, the swap move is always accepted—assuming, of course, that the move either lowers the energy of the system or that the demon has enough energy in store to compensate it [
30]. This makes the MC algorithm extremely efficient. Below, we summarize the MC algorithm:
Start the system in a given initial configuration and set the demon energy to 0.
Calculate the total energy.
Choose two distinct macrocells (the first one can be chosen following some order, and the second one randomly) labeled A and B.
Select two distinct levels proportionately to their presence inside the two macrocells (biased selection). This is done by generating two random numbers between 0 and , and determining which intervals defined by and they fall into.
Perform the trial moves and calculate the change in the total energy of the system:
If , accept the move and update the energy of the reservoir as . Otherwise, reject the move and go back to step 3.
Go through all macrocells once (this defines one MC step).
Recalculate the potential
by numerically integrating Equation (
7) and return to step 2.
We start by testing our MC algorithm on simple one- and two-level waterbag distributions for which the LB entropy function can be maximized exactly. We undertake this investigation within the framework of the one-dimensional self-gravitating model (ODSGM), a model that has been extensively studied in the field of stellar dynamics since the seminal works of Lecar [
33] and Hohl [
34,
35], up to the more recent works of Miller [
36,
37,
38]. It has also found applications in cosmological models explored by Joyce and collaborators [
39,
40,
41]. Despite its simplification compared to real three-dimensional gravity, the ODSGM already contains the main aspects of the gravitational problem—long-range potential [
28,
42], collective motion, particle-wave interactions, which lead to collisionless relaxation [
23,
43,
44], etc. For this study, we opt to use the notation
instead of
to specify the position of a particle in the phase space. A notable practical advantage of the ODSGM is the absence of singular particle–particle interaction since the two-body potential is
, facilitating both theory and simulations [
39,
40].
We begin by studying the equilibrium state to which ODSGM relaxes from an initial single-level waterbag distribution:
where
represents the Heaviside function,
and
. The symmetry of the distribution results in a null total linear
momentum, with the only conserved quantity beside the total mass
being the total energy
. The LB equilibrium is obtained by numerically solving Equations (
6) and (
7) with constraints given by Equations (
4) and (
5). A perfect agreement between the exact (numerical) solution and the MC simulation is demonstrated in
Figure 3, validating the algorithm presented above. We next perform a similar calculation for a two-level waterbag distribution:
where
,
,
and
. Once again, we obtain perfect agreement with simulation results;
Figure 3 and
Figure 4.
We next look at the one-level and two-level distributions, with zero density at the center, see
Figure 5. For the one-level distribution with the hole at the center we use:
where
,
and
. For a two-level waterbag:
where
,
,
,
, and
. We observe that for these initial distributions, the LB equilibrium states do not have holes in the center—instead, we see that the central region is most densely populated according to LB theory (
Figure 6). Again, we observe a perfect agreement between the numerical LB entropy maximization and our MC simulations.
Considering that the validity of our algorithm is now fully established, we are now poised to leverage this approach for the analysis of continuous initial distributions, for which the numerical solution of the LB equations is no longer possible. As a demonstration of the applicability of the MC method, we consider an initial distribution that has a parabolic profile in both position and velocity:
with particles confined in a rectangle of
and
. To perform MC simulations, this distribution is discretized into
levels, randomly generated from the range of 0 to
, where
corresponds to the highest value of the initial distribution Equation (
26). The results of the equilibrium particle distribution obtained using our MC algorithm are presented in
Figure 7 and
Figure 8.
It is interesting to compare the equilibrium distributions predicted by the LB theory with the results of molecular dynamics (MD) simulations. In the MD simulations, we start by distributing
N particles according to Equation (
26) and then evolve the system using Newton’s equations of motion until a stationary state is reached.
Figure 7 and
Figure 8 present both the initial and final stationary particle distributions calculated using LB theory (MC simulations) and MD. We observe that while the initial distributions are identical, the final stationary state to which the system evolves is very different from the predictions of LB theory. While LB theory produces a slowly decaying tail in the particle density distribution function—see
Figure 7,
Figure 8 and
Figure 9—MD has a very sharp decay consistent with the resonant core-halo theory [
45]. Furthermore, the MD distribution in the core region clearly shows incomplete relaxation—the particle density remains depleted in the core, while LB theory predicts a complete population inversion in which higher density levels predominantly occupy the core region. These findings perhaps are not very surprising in view of the earlier work on gravitational systems and on systems with long-range interactions in general [
46,
47]. Still, there was an expectation that the increased complexity of the multilevel distribution might help the system to relax to the LB equilibrium. Instead, we see that relaxation remains incomplete.