Violent relaxation in one-dimensional self-gravitating system: deviation from the Vlasov limit due to finite-$N$ effect
Tirawut Worrakitpoonpon





[Corresponding author: ]worraki@gmail.com
Institute of Science, Suranaree University of Technology, Nakhon Ratchasima 30000, Thailand





We investigate the effect of a finite particle number N on 
the violent relaxation leading to the Quasi-Stationary State (QSS)
in a one-dimensional self-gravitating system.
From the theoretical point of view, we demonstrate that
the local Poissonian fluctuations embedded in the initial state
give rise to an additional term proportional to 1/N
in the Vlasov equation. This term designates the strength of
the local mean-field variations by fluctuations.
Because it is of the mean-field origin, we interpret it differently
from the known collision term in the way that it
effects the violent relaxation stage.
Its role is to deviate the distribution function  
from the Vlasov limit, in the collisionless manner,
at a rate proportional to 1/N while the 
violent relaxation is progressing. This hypothesis is tested
by inspecting the QSSs in simulations of various N.
We observe that the core phase-space density 
can exceed the limiting density deduced from the Vlasov
equation and its deviation degree is in accordance
with the 1/N estimate.
This indicates the deviation from the 
standard mean-field approximation of the violent relaxation process
by that 1/N term.
In conclusion, the finite-N effect has a significant contribution
to the QSS apart from that it plays
a role in the collisional stage that takes place long after.
The conventional collisionless Vlasov equation might not be able 
to describe the violent relaxation of a system of particles properly
without the correction term of the local finite-N fluctuations.






Violent relaxation in one-dimensional self-gravitating system: deviation from the Vlasov limit due to finite-N effect
    Tirawut Worrakitpoonpon
    
=====================================================================================================================





§ INTRODUCTION
 
Long-range interacting systems 
constitute a major branch in statistical mechanics and 
attract much interest due to their complex nature
(see <cit.> for review).
In general, an out-of-equilibrium state
evolves first through a stage of violent relaxation in which
the dynamics is governed by a self-consistent mean field
with a time scale of the order of dynamical times. 
This stage is otherwise known as the collisionless evolution 
due to the lack of the collisional interaction between the particles.
That collisionless dynamics during that stage 
can be described by the Vlasov equation.
At the end of the violent relaxation, 
the system settles into a the Quasi-Stationary State (QSS).
A descriptions of the violent relaxation and
the resulting QSS have been proposed by 
Lynden-Bell (hereafter LB) <cit.>,
initially aiming to explain the inconsistency between the
light distribution of an elliptical galaxy and
the isothermal sphere profile that is the 
thermal equilibrium analogue for the three-dimensional
self-gravitating system. 
That work estimated the violent relaxation time scale 
by the rate of change of the mean field 
and proved that the time scale was of the order of
dynamical times regardless of the system's physical parameters.
This allows a galaxy to be established as a QSS within
a time scale shorter than the age of the Universe.
The statistical description of the 
QSS differs from the thermal 
equilibrium which represents a state of
maximized Boltzmann entropy. Since then, 
the LB theory has became a central pillar 
in the field of astrophysics and the studies of systems governed
by long-range interactions.


While the Vlasov equation addresses the collisionless evolution 
of the long-range interacting systems in the fluid limit, 
the violent relaxation is mostly investigated using systems of particles
in which the time integration of the particle dynamics is simple
compared to the direct integration of the Vlasov equation.

However, the Poissonian fluctuations embedded in particle systems
must be properly considered. A finite particle number N 
constitutes the collision term proportional to 1/N 
that arises from the two-point correlation function
in the BBGKY hierarchy
<cit.>.
It introduces a collisional effect 
to the QSS and drives it to thermal equilibrium with
a time scale increasing with N as verified in
the one-dimensional self-gravitating system
<cit.>
and the Hamiltonian Mean-Field (HMF) model 
<cit.>. 

Although the knowledge of the finite-N effect on
the long-term evolution to the thermal equilibrium
is well established theoretically and numerically,
its impact on the preceding violent relaxation
is not well investigated. Studies of the HMF model 
have mentioned the finite-N effect on the QSS with 
a deviation of the QSS from the prediction based on
the fluid limit 
<cit.>.
The finite-N effect has also been addressed
for a system with attractive power-law interaction
and the correction term to the Vlasov equation
was required in order to describe the particle dynamics
properly <cit.>.
The N-dependence of the QSS is interesting 
since it was not described by the original
LB theory which considered the fluid analogue.
Although the finite-N effect was not taken into account,
the LB theory was shown to be applicable in
a limited number of situations. In most cases, the systems underwent
incomplete relaxation that prevented them to achieve the LB configuration
<cit.>.

The main purpose of this study is 
to better understand, theoretically and numerically, 
the involvement of the finite-N effect in the violent
relaxation and the resulting QSS in the one-dimensional 
self-gravitating system. The article is organized as follows. 
In Sec. <ref>, the analysis of the Vlasov
equation that incorporates the local Poissonian fluctuations 
is performed and the relevance to the finite-N QSS 
is discussed.
Sec. <ref> details of the 
initial condition for the simulations, the units, and 
the simulation method. In Sec. <ref>,
the numerical results and the comparison with
the theoretical framework are presented.
Finally, Sec. <ref> concludes this study.






§ ANALYSIS OF THE VLASOV EQUATION WITH LOCAL FLUCTUATIONS
 

We recall the Vlasov equation that governs 
the collisionless evolution of the distribution function 
f(𝐱,𝐯)

    df/dt = ∂ f/∂ t+
    𝐯·∂ f/∂𝐱+
    𝐚·∂ f/∂𝐯 = 0

where 𝐚 is the self-consistent mean-field acceleration,
an acceleration from all phase-space elements of the system,
derived from the gradient of the mean-field potential Φ as

    𝐚(𝐱) = -∇_𝐱Φ.

By definition, Φ can be written in terms of
the integral of f as

    Φ(𝐱) = ∫∫φ(𝐱,𝐱')f(𝐱',𝐯')
    d𝐱'd𝐯'

where φ (𝐱,𝐱') is the pair potential
at 𝐱 from the mass element at 𝐱'.
According to the LB analysis
<cit.>, there are two types of system
distribution functions: the fine-grained and the coarse-grained
distribution functions. The former one represents
the phase-space density of a system segmented into
smallest divisible phase-space elements, or microcells,
each of which is so small that we can approximate its density
to be uniform and non-zero inside when occupied and
zero otherwise. The latter one is the average of
the phase-space density inside a macrocell, a cell that consists of
multiple microcells. The statistical equilibrium
is attainable on the coarse-grained level, whereas the fine-grained
level never attains because fine-grained entropy does not increase.
Throughout this section, all distribution
functions are the coarse-grained ones.
In simulation aspect, the measured distribution function
and the other profiles derived from it can be considered as 
coarse-grained properties as they
are computed with a finite bin size.

If f(𝐱,𝐯) denotes the point-wise
one-particle distribution function, Eq. (<ref>)
is known as the Klimontovich equation. 
In the case of Newtonian gravity,
Eq. (<ref>) is coupled with the Poisson equation
and these equations are called
collectively the Vlasov-Poisson equations.

For an initial state consisting of N particles,
we express the system distribution function f in the form

    f=f̅+δ f,

where the first term corresponds to the distribution function
in fluid limit and the second term represents
the local variation due to the Poissonian fluctuations
(see <cit.> for review).
The variation term vanishes by ensemble average,
i.e., <δ f>=0, so <f> converges to f̅.

In statistical sense, the ensemble average of
a quantity or a function can be obtained by averaging it
over infinite number of realizations
starting with initial conditions sampled from the same
distribution function.
In practice, the ensemble average over hundreds of
realizations is found sufficient for reliable
distribution functions of the QSSs.

The fluctuations of f lead to the fluctuations of the 
mean-field acceleration δ𝐚 accordingly,
i.e., 𝐚=𝐚̅+δ𝐚
where 𝐚̅ is the mean-field acceleration 
from f̅.
From Eqs. (<ref>) and (<ref>), 
we can deduce that δ𝐚 is 
an integral of δ f
which vanishes after averaging in the same way as δ f.
All variation terms are scaled by 1/√(N) as they
originate from the Poissonian noise. By substituting
the perturbed quantities into Eq. (<ref>)
and performing the ensemble average, the perturbed
Vlasov equation becomes

    <df̅/dt> = 
    -<δ𝐚·∂δ f/∂𝐯>
    ∼1/N,

where

    df̅/dt=
    ∂f̅/∂ t+
    𝐯·∂f̅/∂𝐱+
    𝐚̅·∂f̅/∂𝐯.


Although the ensemble average eliminates all terms involving  
first-order fluctuations in Eq. (<ref>),
a non-vanishing second-order term proportional to 1/N remains,
which is of the mean-field origin and embedded in the initial state.

This analytical result may appear at first analogous to 
the one obtained
from the quasi-linear analysis from which the 1/N term
is recognized as the collision term
<cit.>.
In this work, however, we interpret the role of this term and its 
intrinsic properties differently.
Conventional analysis expresses the two
perturbative terms on the right-hand side in forms of
the spatio-temporal Fourier transforms of δ f, 
and the ensemble average of their product corresponds to
the spatio-temporal two-point correlation function. 
This is the result derivable from the standard
Lenard-Balescu equation <cit.>.

We have an alternative view on that term starting from a
simple dynamical perspective that δ𝐚
which involves the integral of δ f over the entire
phase space, designates how large the local fluctuation of
mean-field acceleration is at any location.
The product δ fδ𝐚 in that term
can then be regarded as the local fluctuations of accelerations
amplified by local variations δ f, which
can also be interpreted as the strength
of the deviation from mean-field approximation.

As the seed of this term, namely δ𝐚,
is of mean-field origin, we hypothesize that it is associated with
the collisionless relaxation to the QSS rather than
the collisional relaxation to thermal equilibrium.
Broadly speaking, this term can be considered as the mean-field
correction term in the Vlasov equation due to finite N.

The violent relaxation incorporating this term
can be described as follows. In a finite-N system, there
exists the mean-field fluctuations with amplitudes proportional
to 1/√(N) in addition to the mean field from the fluid limit.
This combined mean field 
self-consistently governs the early dynamics altogether
and it deviates the distribution function
f̅ from the Vlasov limit
at a rate proportional to 1/N.
The N-dependence arising from the proposed term, which
is of mean-field origin, is unlike the N-dependent
collisional effect that involves the
two-point correlation and comes into play much later.

The proposed mechanism is therefore a distinct 
process from the collisional evolution.
Because the violent relaxation time scale
does not depend on N, the accumulated deviation at the end
is expected to scale with 1/N.
In other words, such effect manifests the N-dependence 
of the QSS which is not foreseen by the original LB theory.

The 1/N term can otherwise be regarded
as the mean-field correction term of the Vlasov equation 
when considering finite-N systems.

As for the collision term, this fluctuation term
vanishes as N→∞. In that limit, the evolution
follows strictly the Vlasov equation, i.e.,
df/dt=df̅/dt=0. This happens when one directly solves
the Vlasov equation for the evolution of f, either
numerically or semi-analytically
<cit.>.
Note that this estimate of the 1/N term
is valid for any long-range interacting systems as long as
the Poissonian fluctuations are present.
Although applying the perturbative analysis to the
Vlasov equation yields the analytical results
that are analogous, in a broad sense, to those
involving the collision term from the BBGKY hierarchy,
the scope of the
applicability of our analysis is more limited.
Most importantly, the violent relaxation, which
is our focus here, leads to the QSS
that differs greatly from the initial state in most of
the cases and the system
evolves rapidly in a short time scale.
The full perturbative analysis therefore fails to
be valid. As a matter of fact, our analysis
only provides an estimate of the rate of the departure
from the Vlasov limit starting from the non-equilibrium
initial condition, via the analysis that yields lastly
Eq. (<ref>).
We do not provide a full description of how the
perturbation (or δ f) evolves, nor how the
final QSS appears.
On the contrary, the collisional term drives
the QSS slowly and quasi-statically through the series of the
dynamical equilibria in the course of the thermalization
and the system eventually attains thermal equilibrium.


In a confined system such as the HMF model, the term δ f 
can be assumed to take a form of the standing wave. 
A recent work by <cit.> performed the 
quasi-linear analysis in that way and investigated
its effect on the QSS. The results were compared 
with N-body simulations and it was found that the analysis
was applicable in a certain range of energy. 
Otherwise, a large deviation from the theory was found.
It may appear that the starting point 
and the purpose of this work are similar to ours:
the perturbative form 
of the distribution function is applied 
to the Vlasov equation and 
its consequence to the QSS is examined, 
but the in-depth details are different.

The major difference is that
they applied the analysis to the regime in which the
QSS did not differ much from the initial state.
In our case, we start from
the far out-of-equilibrium initial state.
Hence, our analysis does not aim to describe the
detailed phase-space density distribution of the QSS
via δ f, but we
only describe how far the global properties of
the phase-space density
of the QSS can deviate from the Vlasov limit.

Another difference is that their work regarded 
the perturbations in the fluid limit
while in our study, those perturbations represent
the local finite-N fluctuations.
Therefore, the resulting QSS from our analysis
is the QSS with the finite-N effect involved
and the test with the N-body simulations
is straightforward.

 




§ SIMULATION DETAILS, INITIAL CONDITION AND UNITS




We test the hypothesis of the 1/N term that arises from
the above analysis in the N-body simulations
of one-dimensional self-gravitating system 
of particles of identical mass m. This model
is also known as the self-gravitating sheet model.
The Newtonian gravitational acceleration 
from the particle j located at x_j
acting on the particle i at x_i reads

    a_ji = gmx_j-x_i/|x_j-x_i|

where g is the gravitational constant. 
The force expressed in Eq. (<ref>) is 
the mutual attractive constant force regardless of separation.
This facilitates the numerical integrations of the equations of motion 
of particles as they can be carried out in the exact way up to
the crossing where the forces on the encountering
particles have to be redetermined.
The most optimized algorithm to handle this is
the heap-based algorithm in which the crossing times
of all pairs of adjacent particles are stored and managed
in a data heap <cit.>.
The accuracy is such that the deviation  
of the total energy is comparable to the machine error
at each crossing. At the end, the deviation of the total
energy is less than 10^-6%.

For the initial condition, we adopt 
the rectangular waterbag configuration in phase space 
in which the initial position and velocity of particles 
are randomly distributed in the range of (-x_0,x_0) 
and (-v_0,v_0), with uniform probability. 
Thus, the uniform initial phase-space density (or f_0) 
inside the waterbag reads

    f_0 = M/4x_0v_0

where M=mN is the total mass and the phase-space density
is zero outside. By this construction, the phase-space
distribution is a two-level function 
which facilitates the analysis as f_0 
is the integral of motion and it also indicates
the maximum attainable phase-space density according
to the LB theory. We choose to parameterize  f_0 by
the initial virial ratio b_0≡2T_0/|U_0|,
where T_0 and U_0 are the initial kinetic and potential
energies of the system, respectively, as 

    f_0=√(M/2gL^3b_0)

where L=2x_0 is the initial system size.
Parameterizing the initial state by b_0 helps to 
identify easily how far the initial state is from
the virial equilibrium. Throughout this study,
we use a system of units where the gravitational constant g, 
the total mass M,
and the system initial length L are unity. The unit of time
t_d is chosen to be 

    t_d = √(L/gM)

which corresponds to the free-fall time of the b_0=0 system.
After all particles are randomly sampled
in phase space following the prescriptions above,
we adjust all positions and velocities by the constants 
so that the center of mass is static at the origin.
Relocating the center of mass in this way is equivalent
to translating the waterbag in phase space while 
the initial degree of randomness, which is our central interest,
is preserved.


Because we study the subtle effect from the random microscopic 
Poissonian fluctuations, which vary between realizations,  
to the macroscopic QSS, the simulation with 
a considerable number of realizations for each 
b_0 and N is preferable to 
obtain the dependable results. For this reason, 
we simulate 1000, 500, 250, 150, 100 
and 50 different 
realizations for N=1000, 2000, 4000, 8000, 16000 and 
32000, respectively, and the ensemble average 
can be applied if necessary. 



§ NUMERICAL RESULTS





 §.§ Waterbag evolution and phase-space density of QSS






In this section, we consider first of all the evolution
of the waterbag. Shown in Fig. <ref> is
the evolution of the waterbag in phase space for
b_0=0.1 and N=32000, drawn from a single realization.
The initial waterbag evolves to the QSS which,
as seen at 300  t_d, consists of a
dense core and surrounding diluted halo achieved
from wound filaments seen in the snapshots before.
This kind of structure is known as the core-halo structure.
The attainment to this structure can be described
by the Vlasov dynamics as follows.
In brief, the core phase-space density reflects
the density of the initial waterbag, namely f_0, while the 
filament densities do not. This is because filaments
become finer while they are winding, 
lowering their coarse-grained phase-space density with time
even though the fine-grained density is preserved to f_0.
On the contrary, the core does not stretch and wind,
so its phase-space density remains at f_0.
The Vlasov description for the collisionless
relaxation is deemed to be valid until filaments
become discontinuous due to limited resolution from finite N
(see <cit.> for review).
We observe that filaments are still intact until 80  t_d.
At 180  t_d, they start to dissolve, but their trace 
in the halo is still observable.
At 300  t_d, filaments are more assimilated.
We also spot a number of persistent holes in the halo, which emerge
between the core and the winding filaments.
The presence of phase-space holes does not significantly
affect the analysis because they are situated well outside
the core, whereas we will examine the core
to test our hypothesis of the 1/N term.




From the visual inspection of the
waterbag evolution, we continue with an inspection
of the phase-space density distribution of the QSS.
Shown in Fig. <ref> is the 
ensemble-averaged phase-space density as a function of the energy 
per unit mass f(ϵ) at 300  t_d for b_0=0.3 and 0.8,
and for various N.
To compute f(ϵ), we first determine the energy distribution
of particles F(ϵ) such that ∫ F(ϵ) dϵ= M. 
Then, the conversion of the distribution from the energy domain
(or F(ϵ)) to the phase-space domain (or f(ϵ))
is carried out by the relation

    F(ϵ) = g(ϵ)f(ϵ)

where g(ϵ) is the density of state 
at energy ϵ which can be computed from

    g(ϵ) = ∫_Φ_0^ϵ2√(2)/a(Φ)√(ϵ-Φ)dΦ.

In Eq. (<ref>), a(Φ)=dΦ/dx 
is the magnitude of the acceleration.
The lower integral limit Φ_0 corresponds to
the minimum potential energy which is located at x=0.
For simplicity, we can set Φ_0=0
without loss of generality. We determine the conversion
factor g(ϵ) by the numerical integration
of Eq. (<ref>). The ensemble-averaged
f(ϵ) is obtained by averaging f(ϵ) from
each individual realization, accompanied with the error
bars estimated from the standard deviation of the mean.

In the plot, we provide f_0 for comparison.
The advantage of inspecting the energy distribution 
in phase-space domain 
is that it can be directly compared to f_0 
which is the conserved quantity in the collisionless evolution
and indicates, in principle, the maximum attainable
phase-space density. Thus, the degeneracy of the system 
and the deviation from that limit can be analyzed.
In addition, the robustness
of the ensemble-averaged f(ϵ) is
tested by the bootstrapping method and the results
are reported in Appendix <ref>.

Because we aim to inspect the 
subtle properties of the QSS, to ascertain that
the state at the end appropriately represents the 
QSS is crucial. On the one hand, we assure that
the violent relaxation is put to an end well
before that time. According to past studies,
it was demonstrated that the violent relaxation that was
probed by the virial ratio terminated within
100  t_d <cit.>.
On the other hand, the fact that the collisional effect,
which also originates from the finite-N fluctuations
in the initial state according to the BBGKY hierarchy, 
is negligible for our conclusion is equally important. 
In a past research work <cit.>,
the collisional relaxation to thermal equilibrium
was tracked by a set of order parameters which were
defined to be zero at the
thermal equilibrium and non-zero at the QSSs.
The evolution of order parameters suggested that
the systems reached the QSS within 100  t_d,
in line with other studies, and remained in the
QSSs for a long period of time. 
With N=800, the order parameters
saturated well in the QSSs and they did not
start to evolve toward the thermal equilibrium 
earlier than 10^4 t_d.
In this study, we choose to inspect the systems with N
higher than that case and the framed time scale is well
before the collisional time scale.
It is worth noting that the BBGKY
hierarchy does not signify the total absence of the 
collisional effect in the time frame considered.
The collisionality is also involved since the start
as with our proposed 1/N term.
Nevertheless, we can safely neglect the collisional
effect for our results because, according to past studies,
it becomes significant at much later time.
Another point to be reassured is that,
albeit the presence of persistent phase-space holes 
(see Fig. <ref>), f(ϵ) is well stationary
in the last 100  t_d, which is the time frame of
consideration for further analysis.


From the f(ϵ) plot, the QSS exhibits
two distinct components: the inner part has f(ϵ) 
close to f_0 while the outer part at higher ϵ
has much lower density. This form corresponds to the core-halo
structure as shown in Fig. <ref> and the
violent relaxation leading to such structure can be explained
by the Vlasov equation. An alternative explanation for the
fact that the core ends up in the degenerate state,
which is the state of lowest possible energy, is given
by the theory of the parametric resonance
(see <cit.> for review). 
That theory was initially developed for explaining
the halo formation in an ion system
<cit.>, but the adaptation to
self-gravitating systems was also considered
<cit.>.
During the violent relaxation,
quasi-periodic oscillation of potential field
caused by a relaxing system can be in resonance
with some particles. As a consequence, those particles
gain the energy and are displaced further away,
constituting the halo.
The energy transfer halts when the core attains 
the degenerate state whose density is limited by f_0.
The parametric resonance in one-dimensional
system was validated by the oscillation of the
position of a test particle that was in phase with the system
oscillation, so it was able to gain the energy and become part
of halo <cit.>.
The fact that the core ends up in the LB degenerate state
has been verified in past works on the self-gravitating sheet model
<cit.>
and on other long-range interacting systems
<cit.>,
and our frame of work is based on that fact.
Because the filaments are still intact
at ∼ 80  t_d, it is suggested that the parametric resonance
takes place at the later time.

The continuum distribution function of the LB degenerate limit
is the two-level step function given by

    f(ϵ)=
    f_0,    if Φ_0 < ϵ < ϵ_f
    
    0,    otherwise

where Φ_0 is the minimum energy and ϵ_F is
the Fermi energy depending on the mass and f_0.

With a close inspection on the core phase-space density, 
we find that f(ϵ) can reach the value above f_0  
and this excess is more noticeable for lower N.
If N is higher, the core f(ϵ)
is more uniform and closer to f_0, albeit with fluctuations,
and f(ϵ) falls more sharply 
to the halo level beyond the core energy limit.
This implies that the core is closer to the LB
degenerate f(ϵ).
In fluid limit, the Vlasov equation
suggests that the core distribution function takes
the degenerate form (<ref>).

That the deviation from the degenerate distribution
function is more evident if N is lower can be
explained that the fluctuations of mean-field
accelerations are higher in magnitude. This causes larger
deviation of core particle trajectories from Vlasov
mean-field trajectories.
The phase-space density excess in Fig. <ref> 
is out of the scope of the
conventional violent relaxation theory because,
according to the description by that theory, 
f_0 marks the maximum allowed phase-space density.
We hypothesize that its origin is from the 
additional term in Eq. (<ref>).
It turns out that such term causes the QSS to be
more concentrated than f_0 in phase space, specifically
in the core part. Note that our study adopts
a simple waterbag that consists
of particles of identical mass and has a single uniform
density, so a single core emerges.
In a more complicated case, multiple cores
can arise <cit.>.

We remark that the phase-space distribution
for b_0=0.8 and N=32000 exhibits a plateau at
high energy which contradicts the Penrose criterion
from which a single maximum is allowed in the stationary state.
We speculate that
the emergence of this component is attributed to
the parametric resonance that favors a certain range
of the particle energy. These particles 
therefore constitute an overpopulated component at
the high energy range.
Another possibility is that, in continuity
with the plot of phase-space configuration in Fig. <ref>,
the underpopulated region is attributed to the voids
that harbor no particles.
Local maxima in the halo
energy range were also spotted in past works
<cit.>
but this component was found to be meta-stable as
it disappeared when the system reached the thermal equilibrium 
<cit.>.



 §.§ N-dependence of core phase-space density





In continuity with the observed excess of the phase-space
density of the core in Sec. <ref>, we will inspect
it to a greater extent in this section. First of all,
we verify if this excess 
is robust in the QSS and it is really achieved 
from the violent relaxation governed by mean field
with fluctuations. To do so, we define the 
core phase-space density f_c to be the averaged
f(ϵ) in the range of energy where f > 0.99f_0,
and the time evolution of 
f_c for b_0=0.3 and for various N is 
plotted in Fig. <ref>.
We observe that f_c increases from f_0
most prominently in the first 100  t_d 
which is during the violent relaxation.
After 200  t_d, f_c saturates at the QSS.
A system with low N tends to depart more from 
f_0 whereas for N=32000, f_c remains
close to f_0 all the time.

This plot ascertains, firstly, that the QSS with
the core phase-space density exceeding f_0  
is a robust structure as a consequence of 
the violent relaxation. 
Secondly, the departure from the theoretical 
limiting value in the fluid limit during the
violent relaxation is more noticeable in
a lower N system. These are in line with the analytical 
result in Sec. <ref>.


The deviation of f_c from f_0 can be
investigated in a more quantitative way by defining 

    Δ f ≡f_c-f_0/f_0

to be the phase-space density excess of the core and
the plot of the magnitude of the ensemble-averaged
Δ f (or |<Δ f>|) for different b_0
and N is shown in Fig. <ref>.
The time average is also performed in
the time window of width 50  t_d
from 250-300  t_d to suppress the time fluctuations.
The size of the error bars corresponds to the averaged mean errors
in the same time window.
For each N, |<Δ f>| tends to be higher when b_0
is lower, which implies that a further out-of-equilibrium
initial state leads to a QSS that
is further away from the Vlasov limit.
A waterbag of lower b_0 
collapses to the center more violently
because of lower support from velocity dispersion.
Thus, local fluctuations in the initial state
gather more concentratedly in the core,
leading to a more pronounced
effect from the mean-field fluctuation term.
Despite the variation in magnitude with b_0,
we observe the decrease of |<Δ f>|
with N for all b_0. To verify if the decrease is
in accordance with the 1/N estimate, we perform the
curve-fitting with the 1/N^γ decay, where
γ is a positive real number, for all b_0 and
the results are depicted in Fig. <ref>.
The line designating the 1/N decay is also put
for comparison. The curve-fitting with
the decreasing function that is not restricted
to 1/N allows us to test the validity of the
1/N estimate properly and, if it is about
to deviate from that estimate, to evaluate
how large the discrepancy is.

In all cases, the best-fitting γ is not
far from 1. This affirms 
the reasonable applicability of the 1/N fluctuation term in
the perturbed Vlasov equation
that we proposed in Sec. <ref>.
While it is desirable if we keep increasing N and inspect
further this finite-N effect, we foresee some limitations
by doing so. For N=32000, the deviation is of the order of a
part per thousand or even below, and the deviation
should be weaker if N is higher.
Therefore, the capability to measure it properly is doubtful.

In past studies, the disagreement between
the theoretical prediction of the QSS
based on the continuum limit and 
the N-body simulations has  
been remarked in the HMF model (see, e.g.,
<cit.>).
Because that model is confined in a finite spatial domain, 
the stability analysis of the continuum 
distribution functions for
the most probable state as the QSS is feasible. 

In our case, because a complete stability 
analysis to obtain the prediction of the QSS
is not simple, we choose
to infer from past studies that the core of the QSS
follows the LB degenerate state. By our choice of 
initial conditions, the degenerate distribution function
has a constant phase-space density equal to f_0 and 
we adopt it as the expected value. We not only find
the discrepancy between the theory and the simulations
as the past studies of the HMF model did, but we also 
attempt to understand better the origin of that
discrepancy and its N-dependence. We are able to 
explain the cause of that deviation by the hypothesis of
the 1/N term. If this 1/N term is generic 
for any long-range interacting systems,
we speculate that the deviation of 
the measured parameters in the N-body system
from the prediction in the fluid limit might scale
with N in the same way as our result.









 §.§ Core density profile and its agreement with LB theory







To inspect the core in more details, the ensemble-averaged
density profile of the core ρ (x) is plotted in
Fig. <ref> for some selected b_0 and N.
We define the core to be the component below ϵ_c
at which f(ϵ) falls to 0.99f_0.
We provide the corresponding density profile of the
LB degenerate limit for comparison, which is numerically
calculated from the distribution function (<ref>)
via the Poisson equation. The plot of the LB density profile
in Fig. <ref> is obtained by averaging 
the LB profiles calculated from  
the core mass of each realization.
The bootstrapping test for the robustness of the
core density profile is also performed and
the results are reported in Appendix <ref>.

From the plot, the deviation from the LB degenerate limit
is remarkable when one inspects the density profile
although |<Δ f>| is as small as a few percents
or sub-percent. For b_0=0.3, the simulated
core is evidently more concentrated than the
degenerate limit counterpart.
These two profiles become closer to each other
if N increases from 1000 to 32000. 
For b_0=0.8, the over-dense core is also observed for
N=1000 but the discrepancy is not as large as in 
the b_0=0.3 plot for the same N. 
As N reaches 32000, 
the simulated and theoretical profiles almost coincide.
That the core density profile becomes closer to the LB
degenerate limit for higher N is in coherence
with the f(ϵ) plots in Fig. <ref>.
The plots in Fig. <ref> are another evidences
of the influence of
the proposed 1/N term that deviates the core from
the continuum degenerate limit.

In past studies, a reasonable agreement between
the simulated QSS and the LB profile has been captured
for b_0 close to 1
<cit.>.
In this work, we prove further that the agreement 
with the LB theory, albeit delimited to the core, 
improves if we increase N. That tendency 
is found even if b_0 is considerably far from 1.
Note that a better agreement for a higher N system
may not be realized if we consider the entire system
because the situation is not simple.
This is because the underlying parametric resonance
that brings the core to the degenerate state,
on the other hand, decouples the halo from the core.
Thus, the entire system is not fully mixed.
This scenario is recognized as the incomplete relaxation
which prevents the QSS to take the LB form <cit.>.
Another possibility of the incomplete
relaxation is that the systems consist of particles
of different masses. This leads to the mass segregation
which also prevents the system to
be fully relaxed <cit.>.




Finally, we revisit the ansatz of Teles et al (2011) 
<cit.> which proposed 
that the core-halo structure can be 
well fitted by the three-level distribution function.
In that ansatz, the core phase-space density is fixed to 
f_0 while the constant halo phase-space density,
the core energy limit and the halo energy limit 
are determined self-consistently by the mass and energy
constraints, and the energy of the furthest particle.
Otherwise, the distribution function is zero.

Our results make evident that the simulated f(ϵ) 
is more complex than the three-level configuration 
(see Fig. <ref>). 
We note the intermediate region where f(ϵ)
decreases from the core density to the halo density
that spans a considerable range of ϵ,
and we also note that the profile of this region
varies considerably with N.

Furthermore, we have inspected the density profile of
the core that is restricted to be within 
the near-degenerate
central region and a poor agreement with 
the LB degenerate
distribution function has been observed in most cases.

It is nevertheless worth testing that ansatz
for the simulated core whether the agreement improves
if we escalate the core energy limit
from the value adopted in past sections while 
the core phase-space density is still fixed to f_0.
To do this, we shift the core cutoff energy
to the intermediate region at which f(ϵ) 
falls to 0.1f_0, i.e., ϵ_c,10,
in the attempt to average the over-dense core and 
a part of less dense intermediate region.

The density profile of the core which is now 
constituted from the particles with energies below 
ϵ_c,10
and the corresponding LB degenerate profile 
are plotted in Fig. <ref> for b_0=0.3 and 
0.5, both of which comprise 8000 particles 
that are comparable to the original work 
of Teles et al (2011) <cit.>.
We observe a good agreement with the LB degenerate
density profiles of these central components.
Although we do not investigate further the 
full core-halo distribution function,
our result affirms the applicability of
the two-level function, with modified energy limit,
to the central component.

The agreement found in the b_0=0.5 plot 
is in accordance with <cit.>
whereas we prove that the agreement 
can be found even though b_0 is significantly below.

We note a density tail beyond the
LB profile and it is the part that exhibits
an obvious inconsistency. We speculate that it originates
from the fluctuations of the energies of particles
near the cutoff energy as the particles merely inside
the degenerate profile can possibly have
the energies higher than the cutoff energy,
which are then excluded from the consideration,
and vice versa for the particles merely outside.

Despite this good agreement with numerical results, 
we may have a different opinion because
the importance of the 
finite-N effect might be overlooked.

We demonstrate that such effect causes 
the evident variation with N of the core.
This implies the violation of the collisionless 
Vlasov equation due to the 
finite-N effect that should be of concern.



§ CONCLUSION


We have investigated the significance of
the local finite-N fluctuations embedded in the
initial states of systems of particles
to the QSS that marks the end of a violent relaxation.
Based on the conventional mean-field description of 
the violent relaxation, the particle number N has 
no influence on either the violent relaxation 
time scale or the statistical description of the QSS.
We have identified the finite-N contribution in
the Vlasov equation as the local mean-field
fluctuation term proportional to 1/N. 
This term is  
intrinsically different from the collision term
since it is of mean-field origin unlike the collision
term that involves the two-body interaction.
Therefore, this term is effective in the violent
relaxation in the way that it deviates the system
from the Vlasov limit at a rate proportional to 1/N.
This implies the N-dependent QSS in the way that
the deviation degree from the Vlasov limit is scaled by 1/N.

That hypothesis is tested in the simulations of one-dimensional
self-gravitating systems starting from a simple waterbag initial 
condition of constant phase-space density f_0.
The phase-space density distribution
is a two-level function: it is equal to f_0 inside
the waterbag and 0 elsewhere. In this way,
f_0 can serve as a reference of the degeneracy.

We focus on the core of the QSS which, 
according to past studies, typically attains 
the LB degenerate state as a result of the 
parametric resonance during 
the violent relaxation. 

The LB degenerate limit designates the state 
of lowest possible energy in the fluid limit and
it is unique to f_0 and mass.

We find that the core phase-space density
can exceed f_0 which marks the maximum attainable
value according to the conventional violent relaxation theory. 
By a proper measurement of the phase-space density
excess, we find that the deviation is 
in agreement with the 1/N estimate.  

A larger deviation from the Vlasov limit in a system of lower N
is further affirmed by the inspection of the 
core density profile where the agreement between the
simulated core density profile and the
LB degenerate density profile improves for higher N. 

In conclusion, we underline the importance of the correction term
in the Vlasov equation when dealing with N-body systems.
Although our study provides a reasonably precise description
of the 1/N term that is applicable to the core, the detail
of the entire QSS remains to be resolved.





This research has received funding support from the 
National Science, Research and Innovation Fund (NSRF) 
via the Program Management Unit for Human Resources 
& Institutional Development, Research and Innovation
(PMUB) (grant number B05F640075), and partly by
Thailand Science Research and Innovation (TSRI)
via Suranaree University of Technology (SUT) 
(grant number 179349). Numerical simulations
are facilitated by HPC resources of the National
Astronomical Research Institute of Thailand (NARIT).
The author also thanks useful comments from anonymous reviewers
and C. Herold in improving the manuscript.




apsrev4-2.bst










47
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Campa et al.(2009)Campa,
  Dauxois, and Ruffo]campa_etal_LRreview_2009
  author author A. Campa, author T. Dauxois,  and author S. Ruffo,  @noop  journal journal Phys. Rep.  volume 480, pages 57 (year
  2009)NoStop
[Lynden-Bell(1967)]lynden_bell_1967
  author author D. Lynden-Bell, @noop  journal journal
  Mon. Not. R. Astr. Soc. volume 136, pages 101 (year 1967)NoStop
[Braun and Hepp(1977)]braun+hepp_1977
  author author W. Braun and author K. Hepp,  @noop  journal journal Commun. Math.
  Phys. volume 56, pages 101 (year 1977)NoStop
[Chavanis(2010)]chavanis_2010_jstat
  author author P. H. Chavanis, @noop  journal journal J.
  Stat. Mech. Theory Exp. volume 5, pages
  19 (year 2010)NoStop
[Fouvry and Bar-Or(2018)]fouvry+bar_or_2018
  author author J.-B. Fouvry and author B. Bar-Or, https://doi.org/10.1093/mnras/sty2571 journal journal Mon. Not. R. Astr. Soc. volume 481, pages 4566 (year
  2018)NoStop
[Miller(1996)]miller_1996thermal
  author author B. N. Miller, https://doi.org/10.1103/PhysRevE.53.R4279 journal journal  volume 53,  pages R4279 (year 1996)NoStop
[Tsuchiya et al.(1996)Tsuchiya, Gouda, and Konishi]tsuchiya_et_al_1996
  author author T. Tsuchiya, author N. Gouda, and author T. Konishi, https://doi.org/10.1103/PhysRevE.53.2210 journal journal  volume 53,  pages 2210 (year 1996)NoStop
[Joyce and Worrakitpoonpon(2010)]joyce+worrakitpoonpon_2010
  author author M. Joyce and author T. Worrakitpoonpon, @noop  journal journal J. Stat. Mech. Theory Exp. volume 10,  pages 10012 (year 2010)NoStop
[Yamaguchi(2003)]yamaguchi_2003
  author author Y. Y. Yamaguchi, @noop  journal journal
  Phys. Rev. E volume 68, pages 066210
  (year 2003)NoStop
[Jain et al.(2007)Jain,
  Bouchet, and Mukamel]jain_et_al_2007
  author author K. Jain, author F. Bouchet,  and author D. Mukamel,  @noop  journal journal J. Stat. Mech.
  Theory Exp. volume 2007, pages 11008
  (year 2007)NoStop
[Yamaguchi et al.(2004)Yamaguchi, Barré, Bouchet,
  Dauxois, and Ruffo]yamaguchi+barre+bouchet+dauxois+ruffo_2004
  author author Y. Y. Yamaguchi, author J. Barré, author F. Bouchet,
  author T. Dauxois, and author S. Ruffo, @noop 
  journal journal Physica A volume 337, pages 36 (year 2004)NoStop
[Antoniazzi et al.(2007)Antoniazzi, Califano, Fanelli, and  Ruffo]antoniazzi+califano+fanelli+ruffo_2007
  author author A. Antoniazzi, author F. Califano, author D. Fanelli,  and author S. Ruffo,  @noop  journal journal Phys. Rev. Lett.  volume 98, pages 150602 (year 2007)NoStop
[Campa and Chavanis(2017)]campa_chavanis_2017
  author author A. Campa and author P.-H. Chavanis, https://doi.org/10.1088/1742-5468/aa668e journal journal Journal of Statistical Mechanics: Theory and
  Experiment volume 5, pages 053202
  (year 2017)NoStop
[Santini et al.(2022)Santini, Giachetti, and Casetti]santini_et_al_2022
  author author A. Santini, author G. Giachetti, and author L. Casetti, https://doi.org/10.1088/1742-5468/ac4516 journal journal J. Stat. Mech. Theory Exp. volume 2022, eid 013210 (year
  2022)NoStop
[Gabrielli et al.(2014)Gabrielli, Joyce, and Morand]gabrielli+joyce+morand_2014
  author author A. Gabrielli, author M. Joyce, and author J. Morand, https://doi.org/10.1103/PhysRevE.90.062910 journal journal  volume 90,  eid 062910 (year 2014), https://arxiv.org/abs/1408.0999 arXiv:1408.0999 [cond-mat.stat-mech]
  NoStop
[Yawn and Miller(1997)]yawn+miller_1997
  author author K. R. Yawn and author B. N. Miller, https://doi.org/10.1103/PhysRevLett.79.3561 journal journal  volume 79,  pages 3561 (year 1997)NoStop
[Yawn and Miller(2003)]yawn+miller_2003
  author author K. R. Yawn and author B. N. Miller, https://doi.org/10.1103/PhysRevE.68.056120 journal journal  volume 68,  eid 056120 (year 2003)NoStop
[Arad and Lynden-Bell(2005)]arad+lynden_bell_2005
  author author I. Arad and author D. Lynden-Bell, @noop  journal journal
  Mon. Not. R. Astr. Soc. volume 361, pages 385 (year 2005)NoStop
[Yamaguchi(2008)]yamaguchi_2008
  author author Y. Y. Yamaguchi, @noop  journal journal
  Phys. Rev. E volume 78, pages 041114
  (year 2008)NoStop
[Joyce and Worrakitpoonpon(2011)]joyce+worrakitpoonpon_2011
  author author M. Joyce and author T. Worrakitpoonpon, @noop  journal journal Phys. Rev. E volume 84, pages 011139 (year 2011)NoStop
[Benetti and Marcos(2017)]benetti+marcos_2017
  author author F. P. C.  Benetti and author B. Marcos, https://doi.org/10.1103/PhysRevE.95.022111 journal journal Phys. Rev. E volume 95, eid
  022111 (year 2017)NoStop
[Chavanis(2022)]chavanis_2022PhyA
  author author P.-H. Chavanis, https://doi.org/10.1016/j.physa.2022.128089
  journal journal Physica A volume 606, eid 128089 (year
  2022)NoStop
[Chavanis(2008)]chavanis_serie_3
  author author P. H. Chavanis, @noop  journal journal
  Physica A volume 387, pages 787
  (year 2008)NoStop
[Janin(1971)]janin_1971
  author author G. Janin, @noop  journal journal Astron.
  Astrophys. volume 11, pages 188
  (year 1971)NoStop
[de Buyl and Gaspard(2011)]de_buly+gaspard_2011
  author author P. de
  Buyl and author P. Gaspard, https://doi.org/10.1103/PhysRevE.84.061139 journal journal  volume 84,  eid 061139 (year 2011), https://arxiv.org/abs/1108.0728 arXiv:1108.0728 [cond-mat.stat-mech]
  NoStop
[Colombi and Touma(2014)]colombi+touma_2014
  author author S. Colombi and author J. Touma, @noop  journal journal Mon.
  Not. R. Astr. Soc. volume 441, pages
  2414 (year 2014)NoStop
[Colombi(2015)]colombi_2015
  author author S. Colombi, https://doi.org/10.1093/mnras/stu2308 journal journal Mon. Not. R. Astr. Soc. volume 446, pages 2902 (year
  2015)NoStop
[Barnes and Ragan(2014)]barnes+ragan_2014
  author author E. I. Barnes and author R. J. Ragan, https://doi.org/10.1093/mnras/stt2051 journal journal Mon. Not. R. Astr. Soc. volume 437, pages 2340 (year
  2014)NoStop
[Ragan and Barnes(2019)]ragan+barnes_2019
  author author R. J. Ragan and author E. I. Barnes, https://doi.org/10.1093/mnras/stz1523 journal journal Mon. Not. R. Astr. Soc. volume 487, pages 3618 (year
  2019)NoStop
[Fanelli et al.(2003)Fanelli, Noullez, and Aurell]noullez_fanelli_aurell_2003_heap
  author author D. Fanelli, author A. Noullez,  and author E. Aurell,  @noop  journal journal J. Comp. Phys.  volume 186, pages 697 (year
  2003)NoStop
[Miller et al.(2023)Miller, Manfredi, Pirjol, and  Rouet]miller_et_al_2023
  author author B. Miller, author G. Manfredi, author D. Pirjol, and author J.-L. Rouet, https://doi.org/10.1088/1361-6382/acb8fb journal journal Class. Quantum Grav. volume 40, eid 073001 (year
  2023)NoStop
[Luwel et al.(1984)Luwel, Severne, and Rousseeuw]luwel_et_al_1984
  author author M. Luwel, author G. Severne,  and author P. J. Rousseeuw,  https://doi.org/10.1007/BF00651601 journal journal Astrophys. Space Sci. volume 100,  pages 261 (year 1984)NoStop
[Mineau et al.(1990)Mineau, Feix, and Rouet]mineau_et_al_1990
  author author P. Mineau, author M. R. Feix, and author J. L. Rouet, @noop  journal journal
  Astron. Astrophys. volume 228, pages
  344 (year 1990)NoStop
[Rouet and Feix(1999)]rouet+feix_1999
  author author J. L. Rouet and author M. R. Feix, @noop  journal journal Phys.
  Rev. E volume 59, pages 73 (year 1999)NoStop
[Levin et al.(2014)Levin, Pakter, Rizzato, Teles, and Benetti]levin_et_al_2014_pr
  author author Y. Levin, author R. Pakter,
  author F. B. Rizzato, author T. N. Teles, and author F. P. C. Benetti, @noop
   journal journal Phys. Rep. volume 535, pages 1 (year
  2014)NoStop
[Gluckstern(1994)]gluckstern_1994
  author author R. L. Gluckstern, https://doi.org/10.1103/PhysRevLett.73.1247
  journal journal  volume 73, pages 1247 (year
  1994)NoStop
[Levin et al.(2008a)Levin, Pakter,  and Teles]levin+pakter+teles_2008
  author author Y. Levin, author R. Pakter, and  author T. N. Teles, @noop
   journal journal Phys. Rev. Lett. volume 100, pages 040604 (year
  2008a)NoStop
[Levin et al.(2008b)Levin, Pakter,  and Rizzato]levin+pakter+rizzato_2008
  author author Y. Levin, author R. Pakter, and  author F. B. Rizzato,  @noop  journal journal Phys. Rev. E  volume 78, pages 021130 (year 2008b)NoStop
[Teles et al.(2011)Teles, Levin, and Pakter]teles_levin_pakter_2011
  author author T. N. Teles, author Y. Levin,  and author R. Pakter,  @noop  journal journal Mon. Not. R.
  Astron. Soc. volume 417, pages L21
  (year 2011)NoStop
[Pakter and Levin(2011)]pakter+levin_2011
  author author R. Pakter and author Y. Levin,  @noop  journal journal Phys. Rev. Lett.  volume 106 (year 2011)NoStop
[Tashiro(2019)]tashiro_2019
  author author T. Tashiro, https://doi.org/10.1016/j.physleta.2019.07.017
  journal journal Physics Letters A volume 383, eid 125829 (year
  2019), https://arxiv.org/abs/1903.03307 arXiv:1903.03307
  [cond-mat.stat-mech] NoStop
[Lecar and Cohen(1971)]lecar+cohen_1971
  author author M. Lecar and author L. Cohen,  @noop  journal journal Astrophys. Space
  Sci. volume 13, pages 397 (year 1971)NoStop
[Teles et al.(2010)Teles,
  Levin, Pakter, and Rizzato]teles+levin+pakter+rizzato_2010
  author author T. N. Teles, author Y. Levin,
  author R. Pakter, and author F. B. Rizzato, @noop 
  journal journal J. Stat. Mech. Theory Exp.  volume 05, pages 05007 (year
  2010)NoStop
[Reidl and Miller(1995)]reidl+miller_1995
  author author C. Reidl and author B. N. Miller, @noop  journal journal Phys.
  Rev. E volume 51, pages 884 (year 1995)NoStop
[Hohl and Feix(1967)]hohl+feix_1967
  author author F. Hohl and author M. Feix,  @noop  journal journal Astrophys. J.  volume 147, pages 147 (year
  1967)NoStop
[Yamashiro et al.(1992)Yamashiro, Gouda, and Sakagami]yamashiro+gouda+sakagami_1992
  author author T. Yamashiro, author N. Gouda,  and author M. Sakagami,  @noop  journal journal Progress of
  Theoretical Physics volume 88, pages
  269 (year 1992)NoStop
[Chavanis(2006)]chavanis_2006PhyA
  author author P.-H. Chavanis, https://doi.org/10.1016/j.physa.2006.01.006
  journal journal Physica A volume 365, pages 102 (year
  2006)NoStop







§ BOOTSTRAPPING TEST FOR THE QSS



In this appendix, we report the results of the bootstrapping
test for the robustness of the QSS. We choose the cases with b_0=0.1
and N=1000, which is furthest away from
the Vlasov limit, and b_0=0.1 and N=32000, which 
has the lowest number of realizations. We compute 
the phase-space density distribution f(ϵ)
(Fig. <ref>) and the
density profile of the core ρ (x) (Fig. <ref>)
for 5 different subsamples.
The horizontal line denoting f_0 is also provided in
the f(ϵ) plots. 
In each subsample, 20% of the realizations are randomly
pulled off. Although the variation of the QSS among realizations
is expected, especially for cases with low N, both figures
ascertain that the number of realizations for our study
is sufficient to yield the dependable results as
the computed profiles almost coincide.










