# Cores of Games via Total Dual Integrality,
with Applications to Perfect Graphs and Polymatroids

Vijay V. Vazirani111Supported in part by NSF grants CCF-1815901 and
CCF-2230414. University of California, Irvine

###### Abstract

LP-duality theory has played a central role in the study of cores of games,
right from the early days of this notion to the present time. The classic
paper of Shapley and Shubik [SS71] introduced the “right” way of exploiting
the power of this theory, namely picking problems whose LP-relaxations support
polyhedra having integral vertices. So far, the latter fact was established by
showing that the constraint matrix of the underlying linear system is totally
unimodular.

We attempt to take this methodology to its logical next step – using total
dual integrality – thereby addressing new classes of games which have their
origins in two major theories within combinatorial optimization, namely
perfect graphs and polymatroids. In the former, we address the stable set and
clique games and in the latter, we address the matroid independent set game.
For each of these games, we prove that the set of core imputations is
precisely the set of optimal solutions to the dual LPs.

Another novelty is the way the worth of the game is allocated among sub-
coalitions. Previous works follow the bottom-up process of allocating to
individual agents; the allocation to a sub-coalition is simply the sum of the
allocations to its agents. The natural process for our games is top-down. The
optimal dual allocates to “objects” in the grand coalition; a sub-coalition
inherits the allocation of each object with which it has non-empty
intersection.

##  1 Introduction

The core is a quintessential solution concept in cooperative game theory. It
captures all possible ways of distributing the total worth of a game among
individual agents and sub-coalitions in such a way that the grand coalition
remains intact, i.e., a sub-coalition will not be able to generate more
profits by itself and therefore has no incentive to secede from the grand
coalition.

LP-duality theory has played a central role in the study of the core, right
from its early days all the way to the present time. The early works of
Bondareva [Bon63] and Shapley [Sha65] gave a necessary and sufficient
condition for non-emptiness of the core of a game, namely that it is
“balanced”. The proof of this fact is based on LP-duality theory and uses the
Farkas Lemma. Although this is an insightful fact, the LP employed for proving
it involves exponentially many variables, one for each subset of players.
Consequently, to the best of our knowledge, it did not lead to an efficient
way of computing core imputations for any natural game.

The classic paper of Shapley and Shubik [SS71] introduced the “right” way of
exploiting the power of LP-duality theory for characterizing core imputations.
They did this in the context of the assignment game222The assignment game can
also be viewed as a matching market in which utilities of the agents are
stated in monetary terms and side payments are allowed, i.e., it is a
transferable utility (TU) market. For an extensive coverage of these notions,
see the book by Moulin [Mou14]., for which they showed that the set of core
imputations is precisely the set of optimal solutions to the dual of the LP-
relaxation of the maximum weight matching problem in the underlying bipartite
graph. See Section 3 for details.

It is important to note that this dual LP is small — it has one variable for
each agent and one constraint for each “team”, consisting of one agent from
each side of the bipartition. In general, an optimal solution to such an LP
will be fractional, i.e., it will pick agents fractionally. Clearly, that is
not very meaningful. In retrospect, a key ingenuity of Shapley and Shubik lay
in choosing a game whose LP-relaxation has the following important property:
the polytope defined by its constraints has integral vertices; in this case,
they are matchings in the graph.

Deng et al. [DIN99] distilled the ideas underlying the Shapley-Shubik Theorem
to obtain a general framework, see Section 3.1, which helps characterize the
cores of several games that are based on fundamental combinatorial
optimization problems, including maximum flow in unit capacity networks both
directed and undirected, maximum number of edge-disjoint $s$-$t$ paths,
maximum number of vertex-disjoint $s$-$t$ paths, maximum number of disjoint
arborescences rooted at a vertex $r$ and concurrent games (i.e., general graph
matching games having a non-empty core). This was followed by other works in
the same vein, see Section 2. Similar to the Shapley-Shubik result, a central
fact underlying these successes is the integrality of the vertices of their
polyhedra, which was so far established by showing that the underlying linear
system is totally unimodular, see Definition 1.

We attempt to take this methodology to its logical next step – using total
dual integrality (TDI)333See Definition 2. – thereby addressing new classes of
games which have their origins in two major theories within combinatorial
optimization, namely perfect graphs and polymatroids. In the former, we
address the stable set and clique games and in the latter, we address the
matroid independent set game, given via a rank oracle. For each of these
games, we prove that the set of core imputations is precisely the set of
optimal solutions to the dual LPs. The TDI of the underlying linear systems,
and the consequent integrality of their polytopes, plays a key role in the
proofs. Our results for the perfect graph and polymatroid games are presented
in Sections 4 and 5, respectively.

Perhaps the most intriguing aspect of our work is the novel manner in which we
need to allocate the worth of the game among sub-coalitions. As stated in
Remark 2, whereas in the assignment game this is done via a bottom-up process,
we need to follow a top-down process for both games we study. To the best of
our knowledge, all previous works characterizing the cores of games follow a
bottom-up process, similar to the assignment game.

######  Definition 1.

Let $Ax\leq b$ be a linear system where $A$ is an $m\times n$ matrix and $b$
an $m$-dimensional vector, both with integral entries. $A$ is said to be
totally unimodular (TUM) if every submatrix of $A$ has determinant $0,1$ or
$-1$. If so, the vertices of the polytope of this linear system are integral.

######  Definition 2.

Let $Ax\leq b$ be a linear system where $A$ is an $m\times n$ matrix and $b$
an $m$-dimensional vector, both with rational entries. We will say that this
linear system is totally dual integral (TDI) if for any integer-valued vector
$c^{T}$ such that the linear program

$\max\\{cx:Ax\leq b\\}$

has an optimum solution, the corresponding dual linear program has an integer
optimal solution. If so, the vertices of the polytope of this linear system
are integral.

Note that TDI is a more general condition than TUM for integrality of
polyhedra. If $A$ is TUM then the polyhedron of the linear system $Ax\leq b$
has integral vertices for every integral vector $b$. However, even if $A$ is
not TUM, for specific choices of an integral vector $b$, the vertices of the
polytope of the linear system $Ax\leq b$ may be integral, and TDI may apply in
this situation. It is important to remark that TDI is not a property of the
polytope but of the particular linear system chosen to define it. See [GLS88,
Sch86] for further details.

##  2 Related Works

An imputation in the core has to ensure that each of the exponentially many
sub-coalitions is “happy” — clearly, that is a lot of constraints. As a
result, the core is non-empty only for a handful of games, some of which are
mentioned in the Introduction. A different kind of game, in which preferences
are ordinal, is based on the stable matching problem defined by Gale and
Shapley [GS62]. The only coalitions that matter in this game are ones formed
by one agent from each side of the bipartition. A stable matching ensures that
no such coalition has the incentive to secede and the set of such matchings
constitute the core of this game. Vande Vate [Vat89] and Rothblum [Rot92] gave
linear programming formulations for stable matchings; the vertices of their
underlying polytopes are integral and are stable matchings. More recently,
Kiraly and Pap [KP08] showed that the linear system of Rothblum is in fact
TDI.

To deal with games having an empty core, e.g., the general graph matching
game, the following two notions have been given in the past. The first is that
of least core, defined by Mascher et al. [MPS79]. If the core is empty, there
will necessarily be sets $S\subseteq V$ such that $v(S)<p(S)$ for any
imputation $v$. The least core maximizes the minimum of $v(S)-p(S)$ over all
sets $S\subseteq V$, subject to $v(\emptyset)=0$ and $v(V)=p(V)$. This
involves solving an LP with exponentially many constraints, though, if a
separation oracle can be implemented in polynomial time, then the ellipsoid
algorithm will accomplish this in polynomial time [GLS88]; see below for a
resolution for the case of the matching game.

A more well known notion is that of nucleolus which is contained in the least
core. After maximizing the minimum of $v(S)-p(S)$ over all sets $S\subseteq
V$, it does the same for all remaining sets and so on. A formal definition is
given below.

######  Definition 3.

For an imputation $v:{V}\rightarrow\mathcal{R}_{+}$, let $\theta(v)$ be the
vector obtained by sorting the $2^{|V|}-2$ values $v(S)-p(S)$ for each
$\emptyset\subset S\subset V$ in non-decreasing order. Then the unique
imputation, $v$, that lexicographically maximizes $\theta(v)$ is called the
nucleolus and is denoted $\nu(G)$.

In 1998, [FKFH98] stated the problem of computing the nucleolus of the
matching game in polynomial time. For the unit weight case, this was done by
Kern and Paulusma [KP03], and for arbitrary weights by Konemann et al.
[KPT20]. However, their algorithm makes extensive use of the ellipsoid
algorithm and is therefore neither efficient nor does it give deep insights
into the underlying combinatorial structure. They leave the open problem of
finding a combinatorial polynomial time algorithm. We note that the difference
$v(S)-p(S)$ appearing in the least core and nucleolus has not been upper-
bounded for any standard family of games, including the general graph matching
game.

A different notion was recently proposed in [Vaz22a], namely approximate core,
and was used to obtain an imputation in the $2/3$-approximate core for the
general graph matching game. This imputation can be computed in polynomial
time, again using the power of LP-duality theory. This result was extended to
$b$-matching games in general graphs by Xiaoet al. [XLF21].

Konemann et al. [KTZ21] showed that computing the nucleolus of the constrained
bipartite $b$-matching game, in which an edge can be matched at most once, is
NP-hard even for the case $b=3$ for all vertices. Biro et al. [BKP12] showed
that the core non-emptiness and core membership problems for the $b$-matching
game are solvable in polynomial time if $b\leq 2$ and are co-NP-hard even for
$b=3$.

Deng and Papadimitriou [DP94] define several computational questions regarding
the core, e.g., testing non-emptyness, membership etc. and study them for a
specific game on an undirected graph with positive and negative weight edges.
Granot and Huberman [GH81, GH84] showed that the core of the minimum cost
spanning tree game is non-empty and gave an algorithm for finding an
imputation in it. Megiddo [Meg78] shows how to compute the nucleolus and
Shapley value of this game. Koh and Sanita [KS20] settle the question of
efficiently determining if a spanning tree game is submodular; the core of
such games is always non-empty. Nagamochi et al. [NZKI97] characterize non-
emptyness of core for the minimum base game in a matroid; the minimum spanning
tree game is a special case of this game. For the $b$-matching game, non-
emptyness of the core was established by showing that every optimal dual
solution is a core imputation [Vaz22b]; this involved exploiting total
unimodularity of the constraint matrix.

##  3 The Core of the Assignment Game

In this section, we will give a short overview of the Shapley-Shubik work as
well as the framework of Deng et al. [DIN99], since they provide a valuable
background for the rest of the paper.

The following setting, taken from [EK01] and [BKP12], vividly captures the
issues underlying profit-sharing in an assignment game. Suppose a coed tennis
club has sets $U$ and $V$ of women and men players, respectively, who can
participate in an upcoming mixed doubles tournament. Assume $|U|=m$ and
$|V|=n$, where $m,n$ are arbitrary. Let $G=(U,V,E)$ be a bipartite graph whose
vertices are the women and men players and an edge $(i,j)$ represents the fact
that agents $i\in U$ and $j\in V$ are eligible to participate as a mixed
doubles team in the tournament. Let $w$ be an edge-weight function for $G$,
where $w_{ij}>0$ represents the expected earnings if $i$ and $j$ do
participate as a team in the tournament. The total worth of the game is the
weight of a maximum weight matching in $G$.

Assume that the club picks such a matching for the tournament. The question is
how to distribute the total profit among the agents — strong players, weak
players and unmatched players — so that no subset of players feel they will be
better off seceding and forming their own tennis club.

######  Definition 4.

The assignment game, as defined above, will be denoted by $G=(U,V,E),\
w:E\rightarrow\mathcal{R}_{+}$. The set of all players, $U\cup V$, is called
the grand coalition. A subset of the players, $(S_{u}\cup S_{v})$, with
$S_{u}\subseteq U$ and $S_{v}\subseteq V$, is called a coalition or a sub-
coalition.

######  Definition 5.

The worth of a coalition $(S_{u}\cup S_{v})$ is defined to be the maximum
profit that can be generated by teams within $(S_{u}\cup S_{v})$ and is
denoted by $p(S_{u}\cup S_{v})$. Formally, $p(S_{u}\cup S_{v})$ is the weight
of a maximum weight matching in the graph $G$ restricted to vertices in
$(S_{u}\cup S_{v})$ only. $p(U\cup V)$ is called the worth of the game. The
characteristic function of the game is defined to be $p:2^{U\cup
V}\rightarrow\mathcal{R}_{+}$.

######  Definition 6.

An imputation444Some authors prefer to call this a pre-imputation, while using
the term imputation when individual rationality is also satisfied. gives a way
of dividing the worth of the game, $p(U\cup V)$, among the agents. It consists
of two functions $u:{U}\rightarrow\mathcal{R}_{+}$ and
$v:{V}\rightarrow\mathcal{R}_{+}$ such that $\sum_{i\in U}{u(i)}+\sum_{j\in
V}{v(j)}=p(U\cup V)$.

######  Definition 7.

An imputation $(u,v)$ is said to be in the core of the assignment game if for
any coalition $(S_{u}\cup S_{v})$, the total worth allocated to agents in the
coalition is at least as large as the worth that they can generate by
themselves, i.e., $\sum_{i\in S_{u}}{u(i)}+\sum_{j\in S_{v}}{v(j)}\geq p(S)$.

We next describe the characterization of the core of the assignment game given
by Shapley and Shubik [SS71]. As stated in Definition 5, the worth of the
game, $G=(U,V,E),\ w:E\rightarrow\mathcal{R}_{+}$, is the weight of a maximum
weight matching in $G$. Linear program (7) gives the LP-relaxation of the
problem of finding such a matching. In this program, variable $x_{ij}$
indicates the extent to which edge $(i,j)$ is picked in the solution. Matching
theory tells us that this LP always has an integral optimal solution [LP86];
the latter is a maximum weight matching in $G$.

∑_(i, j) ∈E w_ij x_ij ∑_(i, j) ∈E x_ij≤1 ∀i ∈U ∑_(i, j) ∈E x_ij≤1 ∀j ∈V
x_ij≥0∀(i, j) ∈E

Taking $u_{i}$ and $v_{j}$ to be the dual variables for the first and second
constraints of (7), we obtain the dual LP:

∑_i ∈U u_i + ∑_j ∈V v_j u_i + v_j ≥w_ij ∀(i, j) ∈E u_i≥0∀i ∈U v_j≥0∀j ∈V

######  Theorem 1.

(Shapley and Shubik [SS71]) The imputation $(u,v)$ is in the core of the
assignment game if and only if it is an optimal solution to the dual LP, (7).

By Theorem 1, the core of the assignment game is a convex polyhedron. Finally,
we state a fundamental fact about LP (7); its proof follows from the fact that
the constraint matrix of this LP is totally unimodular.

######  Theorem 2.

(Birkhoff [Bir46]) The vertices of the polytope defined by the constraints of
LP (7) are $0/1$ vectors, i.e., they are matchings in $G$.

###  3.1 The Framework of Deng et al. [DIN99]

In this section, we present the framework of Deng et al. [DIN99], which was
mentioned in the Introduction. Let $T=\\{1,\cdots,n\\}$ be the set of $n$
agents of the game. Let $w\in\mathbb{R}^{m}_{+}$ be an $m$-dimensional non-
negative real vector specifying the weights of certain objects; in the
assignment game, the objects are edges of the underlying graph. Let $A$ be an
$n\times m$ matrix with $0/1$ entries whose $i^{th}$ row corresponds to agent
$i\in T$. Let $x$ be an $m$-dimensional vector of variables and $\mathbb{1}$
be the $n$-dimensional vector of all 1s. Assume that the worth of the game is
given by the objective function value of following integer program.

w ⋅x Ax ≤1 x ∈{0, 1}

Moreover, for a sub-coalition, $T^{\prime}\subseteq T$ assume that its worth
is given by the integer program obtained by replacing $A$ by $A^{\prime}$ in
(3.1), where $A^{\prime}$ picks the set of rows corresponding to agents in
$T^{\prime}$. The LP-relaxation of (3.1) is:

w ⋅x Ax ≤1 x ≥0

Deng et al. proved that if LP (3.1) always has an integral optimal solution,
then the set of core imputations of this game is exactly the set of optimal
solutions to the dual of LP (3.1). So far, for specific games, integrality of
the vertices of the polytope defined by this linear system was established by
showing that its constraint matrix is totally unimodular, as in Theorem 2.

##  4 The Stable Set Game in Perfect Graphs

In Section 4.1 we give the required definitions and facts from the (extensive)
theory of perfect graphs. We will not credit individual papers for these
facts; instead, we refer the reader to Chapter 9 of the book [GLS88] as well
as the remarkably clear and concise exposition of this theory, presented as an
“appetizer” by Groetschel [Grö99].

###  4.1 Definitions and Preliminaries

######  Definition 8.

Given a graph $G=(V,E)$, $\omega(G)$ denotes its clique number, i.e., the size
of the largest clique in it and $\chi(G)$ denotes its chromatic number, i.e.,
the minimum number of colors needed for its vertices so that the two endpoints
of any edge get different colors.

######  Definition 9.

A graph $G=(V,E)$ is said to be perfect if and only if for each vertex-induced
subgraph $G^{\prime}\subseteq G$, $\omega(G^{\prime})=\chi(G^{\prime})$.

Let $\overline{G}$ denote the complement of $G$, i.e.,
$\overline{G}=(V,\overline{E})$, where $\overline{E}$ is the complement of
$E$, with $\forall\ u,v\in V,(u,v)\in E$ if and only if
$(u,v)\notin\overline{E}$. A central fact about perfect graphs is that $G$ is
perfect if and only if $\overline{G}$ is perfect.

######  Definition 10.

$V^{\prime}\subseteq V$ is said to be a stable set in $G$, also sometimes
called an independent set, if no two vertices of $V^{\prime}$ are connected by
an edge, i.e., $\forall\ u,v\in V^{\prime},\ (u,v)\notin E$. Let
$w:V\rightarrow\boldsymbol{\mathit{Q}}_{+}$ be a weight function on the
vertices of $G$. Then the worth of the stable set game is defined to be the
weight of a maximum weight stable set in $G$ and is denoted by $\mbox{\rm
worth}(G)$.

The game defined above is well-motivated: $V$ represents a set of agents from
which a subset needs to be chosen for performing a certain task. For agent
$v\in V$, $w(v)$ specifies the contribution of agent $v$ to performing the
task. An edge $(u,v)$ indicates that agents $u$ and $v$ do not get along, to
the extent that a sub-coalition involving both of them is bound to fail.
Therefore, the most effective sub-coalition is a maximum weight stable set in
$G$. However, as is well known, the maximum stable set problem is NP-hard in
arbitrary graphs, even with unit weights. For this reason, we restrict to
perfect graphs where, as we will see below, it can be solved in polynomial
time.

######  Remark 1.

Since the complement of a perfect graph is also perfect, our result for the
stable set game in perfect graphs yields an analogous result for the clique
game in perfect graphs. The latter game is also well-motivated.

Let $G$ be an arbitrary graph. Clearly any clique in $G$ can intersect a
stable set in at most one vertex, and therefore the constraint in LP (1) is
satisfied by every stable set; note that variable $x_{v}$ indicates the extent
to which $v$ is picked in a fractional stable set. LP (1) contains such a
constraint for each clique in $G$ and is an LP-relaxation of the maximum
weight stable set problem in $G$. However, LP (1) has exponentially many
constraints, one corresponding to each clique in $G$; moreover, it is NP-hard
to solve in general [GLS88].

∑_v ∈V w_v x_v x(Q)≤1 ∀ clique $Q$ in G x_v≥0∀v ∈V

The situation is salvaged in case $G$ is a perfect graph: via a rather long
route, touching profound ideas, e.g., a foray into the Lovasz theta function
and Shannon entropy, one can can show that LP (1) can be solved in polynomial
time, see [GLS88, Grö99]. This requires the use of the ellipsoid algorithm;
unfortunately, at present, an efficient combinatorial algorithm is not known.

Another useful fact is that the constraint system of LP (1) is TDI for perfect
graphs, see Definition 2. As a consequence, the polytope defined by this
system has integral vertices, i.e., they are stable sets. This together with
the previous fact implies that in perfect graphs, the worth of the stable set
game can be computed in polynomial time. Furthermore, since for each vertex-
induced subgraph $G^{\prime}\subseteq G$, $G^{\prime}$ is also perfect, this
holds for the game restricted to $G^{\prime}$ as well.

In order to characterize core imputations of this game, we will need the dual
LP. This is obtain by taking $y_{Q}$ to be the dual variable for the
constraint of LP (1), and is given in LP (1) below.

∑_clique $Q$ in G y_Q ∑_Q ∋v y_Q ≥w_v ∀v ∈V y_Q≥0∀ clique $Q$ in G

The dual LP is solving a clique covering problem. In particular, if $w$ is 1
for every vertex, the minimum number of cliques which cover each vertex at
least once is an optimal solution to the dual – this follows from the fact
that its constraint system is TDI.

###  4.2 Characterizing the Core

An imputation describes how the worth of the game leads to allocation of
“profits” to sub-coalitions. In the assignment game, an optimal dual
distributes the worth of the game among the agents, i.e., the vertices of the
graph. The allocation to a sub-coalition is simply the sum of the allocations
to the agents in it.

In our setting, an optimal dual distributes the worth of the game among the
cliques of $G$. However, dividing the worth $y_{Q}$, given to clique $Q$,
among the agents in the clique in not very meaningful. Instead, for a clique
$Q$, we will view the worth $y_{Q}$ as residing in every sub-clique of $Q$.
However, this is done not by simply duplicating the worth onto each sub-
clique, an operation that is clearly not legitimate. It is done in a more
interesting manner, so as to define allocations to sub-coalitions in a
consistent manner.

Before providing further details, let us establish terminology that is more
appropriate for the setting at hand. We will think of $y_{Q}$ as the
“satisfaction” of clique $Q$. Consider a sub-coalition $T\subseteq V$ such
that $Q\cap T\neq\emptyset$, and let $Q^{\prime}=Q\cap T$. As stated above,
the satisfaction of $y_{Q}$ resides on every sub-clique of $Q$, including
$Q^{\prime}$. The crucial point is that the latter viewpoint will manifest
itself only when we “think” of the sub-coalition $T$. Moreover, $Q$ will
allocate only to its largest sub-clique in $T$. The details given below will
clarify these ideas.

Let $y:(\mbox{cliques in}\ G)\rightarrow\boldsymbol{\mathit{R}}_{+}$ be a
function assigning satisfaction to the cliques of $G$. Define

$\mbox{\rm satisfaction}(G)=\sum_{\mbox{clique $Q$ in}\ G}{y_{Q}}.$

Let $T\subseteq V$ and let $G^{\prime}$ be the subgraph of $G$ induced on $T$.
Each clique $Q$ in $G$ will allocate satisfaction to at most one clique in
$G^{\prime}$ as follows. If $Q\cap T=\emptyset$, then $Q$ will make no
allocation to any clique in $G^{\prime}$. Otherwise, it will allocate $y_{Q}$
to $Q^{\prime}=Q\cap T$. In general, there may be several cliques in $G$ which
may allocate to $Q^{\prime}$. Hence the total allocation to $Q^{\prime}$ is
defined to be

$z_{Q^{\prime}}=\sum_{Q\ {\mbox{in}\ G}:Q\cap T=Q^{\prime}}{y_{Q}},$

where $z:(\mbox{cliques in}\
G^{\prime})\rightarrow\boldsymbol{\mathit{R}}_{+}$ is the function which
assigns satisfaction to the cliques of $G^{\prime}$.

Finally, define

$\mbox{\rm satisfaction}(G^{\prime})=\sum_{\mbox{clique $Q^{\prime}$ in}\
G^{\prime}}{z_{Q^{\prime}}}.$

######  Definition 11.

We will say that imputation $y$ is in the core of the stable set game if

  1. 1.

$\mbox{\rm worth}(G)=\mbox{\rm satisfaction}(G)$

  2. 2.

$\forall G^{\prime}\subseteq G,\ \mbox{\rm worth}(G^{\prime})\leq\mbox{\rm
satisfaction}(G^{\prime})$

######  Remark 2.

The contrast between the way allocations are made to sub-coalitions in the
assignment game and in our game can be summarized as follows: Whereas the
former follows a bottom-up process, we follow a top-down process. To the best
of our knowledge, all previous works characterizing the cores of games follow
a bottom-up process, similar to the assignment game.

######  Theorem 3.

An imputation $y$ is in the core of the stable set game over a perfect graph
if and only if it is an optimal solution to the dual, LP(1).

###### Proof.

($\Leftarrow$) First assume that $y$ is an optimal solution to the dual,
LP(1). By the LP-duality theorem and the TDI of its linear system, $\mbox{\rm
worth}(G)=\mbox{\rm satisfaction}(G)$.

Let $T\subseteq V$ and let $G^{\prime}$ be the graph induced on $T$. By the
allocation procedure described above, for $v\in T$, each clique $Q$ in $G$
containing $v$ will contribute $y_{Q}$ to its largest sub-clique $Q^{\prime}$
in $G^{\prime}$, which will contain $v$. Therefore,

$\forall v\in T:\ \sum_{Q^{\prime}\ni v}{z_{Q^{\prime}}}=\sum_{Q\ni
v}{y_{Q}}\geq w_{v}.$

Therefore $z$ is feasible for the dual LP (1), restricted to $G^{\prime}$.
Hence by the weak duality theorem, we get that the objective value function of
the primal LP restricted to $G^{\prime}$ is bounded by $\mbox{\rm
satisfaction}(G^{\prime})$. Finally, since this restricted system is also TDI,
we get $\mbox{\rm worth}(G^{\prime})\leq\mbox{\rm satisfaction}(G^{\prime})$.

($\Rightarrow$) Next, assume that imputation $y$ is in the core of the stable
set game over a perfect graph. By definition, $\mbox{\rm worth}(G)=\mbox{\rm
satisfaction}(G)=\sum_{\mbox{clique $Q$ in}\ G}{y_{Q}}$.

Let $v\in V$ and let $G^{\prime}$ be $G$ restricted to $v$. By the definition
of core and the way $y$ yields $z$,

$w_{v}=\mbox{\rm worth}(G^{\prime})\leq\mbox{\rm
satisfaction}(G^{\prime})=\sum_{Q^{\prime}\ni v}{z_{Q^{\prime}}}=\sum_{Q\ni
v}{y_{Q}}.$

Therefore, $y$ satisfies the constraint for each $v\in V$ in the dual, LP(1).
Hence $y$ is dual feasible. This together with the previous statement gives
that $y$ is an optimal solution to the dual, LP(1). ∎

######  Corollary 1.

The core of the stable set game over a perfect graph is non-empty.

##  5 The Matroid Independent Set Game, Given via a Rank Oracle

In Section 5.1, we present the necessary background information. Once again,
we will not credit individual papers for these (mostly well-known) facts;
instead we refer the reader to Chapter 10 in [GLS88].

###  5.1 Definitions and Preliminaries

Let $\mathcal{M}=(U,\mathcal{I})$ be a matroid, where $U$ is the ground set
with $n$ elements and $\mathcal{I}$ is its set of independent sets. Let
$w:U\rightarrow\boldsymbol{\mathit{R}}_{+}$ be a function that assigns weights
to the elements. Assume further that we are given a rank oracle $r$: for any
set $S\subseteq U$, $r(S)$ is the rank of $S$, i.e., the size of a maximum
independent set in $S$. The rank function of a matroid is submodular, i.e.,

$r(S)+r(T)\geq r(S\cap T)+r(S\cup T).$

It is well known that a maximum weight independent set can be found via the
following greedy algorithm. Sort the elements of $U$ by decreasing weight, say
$e_{1},e_{2},\ldots,e_{n}$. For $1\leq i\leq n$ do: Pick $e_{i}$ if
$r(e_{1},\ldots e_{i})>r(e_{1},\ldots,e_{i-1})$.

######  Definition 12.

Given matroid $\mathcal{M}=(U,\mathcal{I})$ and weight function $w$ on its
elements, the worth of the matroid independent set game is defined to be the
weight of a maximum weight independent set in $\mathcal{M}$ and is denoted by
$\mbox{\rm worth}(\mathcal{M})$.

We wish to characterize the core imputations of the matroid independent set
game. LP (12) is the LP-relaxation of the maximum weight independent set
problem; here $x_{e}$ is the extent to which element $e$ is picked.

∑_e ∈U w_e x_e x(S)≤r(S) ∀S ⊆U x_e≥0∀e ∈U

The linear system of LP (12) is TDI and therefore its polytope has integral
vertices, i.e., they are independent sets of the matroid. This polytope is
called a polymatroid. Taking $u_{i}$ and $v_{j}$ to be the dual variables for
the first and second constraints of (12), we obtain the dual LP:

∑_S ⊆U r(S) ⋅y_S ∑_S ∋e y_S ≥w_e ∀e ∈U y_S≥0∀S ⊆U

Given an oracle for the rank function $r$, the dual LP can be solved in
polynomial time using the ellipsoid algorithm [GLS88].

###  5.2 Characterizing the Core

The process of allocating the worth of the game to sub-coalitions is similar
to that followed in the stable set game over perfect graphs, i.e., it is top-
down, as described in Remark 2. An optimal dual distributes the worth of the
game among subsets $S\subseteq U$. Once again, dividing the worth $y_{S}$,
given to set $S$, among the agents in the set in not very meaningful. Instead,
again we will view the worth $y_{S}$ as residing in every subset of $S$, and
this viewpoint manifests itself only when we “think” of the sub-coalition $T$.
As before, we will think of $y_{S}$ as the “satisfaction” of set $S$.

######  Remark 3.

An added source of confusion in this setting is that the worth of the game is
distributed among subsets of $U$, and sub-coalitions are also subsets of $U$.
We will try to mitigate this confusion by naming the former $S,S^{\prime}$
etc. and the latter as $T,T^{\prime}$ etc.

Let $y:2^{U}\rightarrow\boldsymbol{\mathit{R}}_{+}$ be a function assigning
satisfaction to the subsets of $U$. Define

$\mbox{\rm satisfaction}(U)=\sum_{S\subseteq U}{y_{S}}.$

Consider sub-coalition $T\subseteq U$. Each subset $S$ of $U$ will allocate
satisfaction to at most one subset of $T$ as follows. If $S\cap T=\emptyset$,
then $S$ will make no allocation to any subset of $T$. Otherwise, it will
allocate $y_{S}$ to the set $S^{\prime}=S\cap T$, i.e., the largest subset of
$S$ in sub-coalition $T$. Once again, there may be several subsets of $U$
which may allocate to $S^{\prime}$. The total allocation to $S^{\prime}$ is
defined to be

$z_{S^{\prime}}=\sum_{S\subseteq U:S\cap T=S^{\prime}}{y_{S}},$

where $z:2^{T}\rightarrow\boldsymbol{\mathit{R}}_{+}$ is the function which
assigns satisfaction to subsets of $T$.

Finally, define

$\mbox{\rm satisfaction}(T)=\sum_{S^{\prime}\subseteq T}{z_{S^{\prime}}}.$

######  Definition 13.

We will say that imputation $y$ is in the core of the matroid independent set
game if

  1. 1.

$\mbox{\rm worth}(U)=\mbox{\rm satisfaction}(U)$

  2. 2.

$\forall T\subseteq U,\ \mbox{\rm worth}(T)\leq\mbox{\rm satisfaction}(T)$

######  Theorem 4.

An imputation $y$ is in the core of the matroid independent set game if and
only if it is an optimal solution to the dual, LP(12).

###### Proof.

($\Leftarrow$) First assume that $y$ is an optimal solution to the dual,
LP(12). By the LP-duality theorem and the TDI of this linear system,
$\mbox{\rm worth}(U)=\mbox{\rm satisfaction}(U)$.

Consider sub-coalition $T\subseteq U$. By the allocation procedure described
above,

$\forall e\in T:\ \sum_{S^{\prime}\subseteq T:\ S^{\prime}\ni
e}{z_{S^{\prime}}}=\sum_{S\subseteq U:\ S\ni e}{y_{S}}\geq w_{e}.$

Therefore $z$ is feasible for the dual LP (12), restricted to $T$. Therefore,
by the weak duality theorem, we get that the objective value function of the
primal LP restricted to $T$ is bounded by $\mbox{\rm satisfaction}(T)$.
Finally, since this restricted system is also TDI, we get $\mbox{\rm
worth}(T)\leq\mbox{\rm satisfaction}(T)$.

($\Rightarrow$) Next, assume that imputation $y$ is in the core of the stable
set game over a perfect graph. By definition, $\mbox{\rm worth}(U)=\mbox{\rm
satisfaction}(U)=\sum_{S\subseteq U}{y_{S}}$.

Let $e\in U$ and let $T^{\prime}$ be the sub-coalition $\\{e\\}$. By the
definition of core and the way $y$ yields $z$,

$w_{e}=\mbox{\rm worth}(T^{\prime})\leq\mbox{\rm
satisfaction}(T^{\prime})=\sum_{S^{\prime}\subseteq T^{\prime}:\ S^{\prime}\ni
e}{z_{S^{\prime}}}=\sum_{S\subseteq U:\ S\ni e}{y_{S}}.$

Therefore, $y$ satisfies the constraint for each $e\in U$ in the dual LP (12).
Hence $y$ is dual feasible. This together with the previous statement gives
that $y$ is an optimal solution to the dual LP (12). ∎

######  Corollary 2.

The core of the matroid independent set game is non-empty.

##  6 Acknowledgements

I wish to thank Gerard Cornuejols, Federico Echenique, Martin Groetschel,
Hervé Moulin and Joseph Root for valuable discussions.

## References

  * [Bir46] Garrett Birkhoff.  Three observations on linear algebra.  Univ. Nac. Tacuman, Rev. Ser. A, 5:147–151, 1946.
  * [BKP12] Péter Biró, Walter Kern, and Daniël Paulusma.  Computing solutions for matching games.  International journal of game theory, 41(1):75–90, 2012.
  * [Bon63] Olga N Bondareva.  Some applications of linear programming methods to the theory of cooperative games.  Problemy Kibernetiki, 10(119):139, 1963.
  * [DIN99] Xiaotie Deng, Toshihide Ibaraki, and Hiroshi Nagamochi.  Algorithmic aspects of the core of combinatorial optimization games.  Mathematics of Operations Research, 24(3):751–766, 1999.
  * [DP94] Xiaotie Deng and Christos H Papadimitriou.  On the complexity of cooperative solution concepts.  Mathematics of operations research, 19(2):257–266, 1994.
  * [EK01] Kimmo Eriksson and Johan Karlander.  Stable outcomes of the roommate game with transferable utility.  International Journal of Game Theory, 29(4):555–569, 2001.
  * [FKFH98] Ulrich Faigle, Walter Kern, Sándor P Fekete, and Winfried Hochstättler.  The nucleon of cooperative games and an algorithm for matching games.  Mathematical Programming, 83(1):195–211, 1998.
  * [GH81] Daniel Granot and Gur Huberman.  Minimum cost spanning tree games.  Mathematical programming, 21(1), 1981.
  * [GH84] Daniel Granot and Gur Huberman.  On the core and nucleolus of minimum cost spanning tree games.  Mathematical programming, 29(3):323–347, 1984.
  * [GLS88] M. Grotschel, L. Lovasz, and A. Schirjver.  Geometric Algorithms and Combinatorial Optimization.  Springer-Verlag, 1988.
  * [Grö99] Martin Grötschel.  My favorite theorem: characterizations of perfect graphs.  1999\.
  * [GS62] David Gale and Lloyd S Shapley.  College admissions and the stability of marriage.  The American Mathematical Monthly, 69(1):9–15, 1962.
  * [KP03] Walter Kern and Daniël Paulusma.  Matching games: the least core and the nucleolus.  Mathematics of operations research, 28(2):294–308, 2003.
  * [KP08] Tamás Király and Júlia Pap.  Total dual integrality of rothblum’s description of the stable-marriage polyhedron.  Mathematics of Operations Research, 33(2):283–290, 2008.
  * [KPT20] Jochen Könemann, Kanstantsin Pashkovich, and Justin Toth.  Computing the nucleolus of weighted cooperative matching games in polynomial time.  Mathematical Programming, 183(1):555–581, 2020.
  * [KS20] Zhuan Khye Koh and Laura Sanità.  An efficient characterization of submodular spanning tree games.  Mathematical Programming, 183(1):359–377, 2020.
  * [KTZ21] Jochen Könemann, Justin Toth, and Felix Zhou.  On the complexity of nucleolus computation for bipartite b-matching games.  In International Symposium on Algorithmic Game Theory, pages 171–185. Springer, 2021.
  * [LP86] L. Lovász and M.D. Plummer.  Matching Theory.  North-Holland, Amsterdam–New York, 1986.
  * [Meg78] Nimrod Megiddo.  Computational complexity of the game theory approach to cost allocation for a tree.  Mathematics of Operations Research, 3(3):189–196, 1978.
  * [Mou14] Hervé Moulin.  Cooperative microeconomics: a game-theoretic introduction, volume 313.  Princeton University Press, 2014.
  * [MPS79] Michael Maschler, Bezalel Peleg, and Lloyd S Shapley.  Geometric properties of the kernel, nucleolus, and related solution concepts.  Mathematics of operations research, 4(4):303–338, 1979.
  * [NZKI97] Hiroshi Nagamochi, Dao-Zhi Zeng, Naohiśa Kabutoya, and Toshihide Ibaraki.  Complexity of the minimum base game on matroids.  Mathematics of Operations Research, 22(1):146–164, 1997.
  * [Rot92] Uriel G Rothblum.  Characterization of stable matchings as extreme points of a polytope.  Mathematical Programming, 54(1):57–67, 1992.
  * [Sch86] A. Schrijver.  Theory of Linear and Integer Programming.  John Wiley & Sons, New York, NY, 1986.
  * [Sha65] Lloyd S Shapley.  On balanced sets and cores.  Technical report, RAND CORP SANTA MONICA CALIF, 1965.
  * [SS71] Lloyd S Shapley and Martin Shubik.  The assignment game I: The core.  International Journal of Game Theory, 1(1):111–130, 1971.
  * [Vat89] John H Vande Vate.  Linear programming brings marital bliss.  Operations Research Letters, 8(3):147–153, 1989.
  * [Vaz22a] Vijay V Vazirani.  The general graph matching game: Approximate core.  Games and Economic Behavior, 132, 2022.
  * [Vaz22b] Vijay V Vazirani.  New characterizations of core imputations of matching and $b$-matching games.  In Foundations of Software Technology and Theoretical Computer Science, 2022.
  * [XLF21] Han Xiao, Tianhang Lu, and Qizhi Fang.  Approximate core allocations for multiple partners matching games.  arXiv preprint arXiv:2107.01442, 2021.


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

Tirawut Worrakitpoonpon<EMAIL_ADDRESS>Institute of Science, Suranaree
University of Technology, Nakhon Ratchasima 30000, Thailand

###### Abstract

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.

##  I Introduction

Long-range interacting systems constitute a major branch in statistical
mechanics and attract much interest due to their complex nature (see [1] 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) [2], 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 [3, 4, 5]. 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 [6, 7, 8] and the Hamiltonian Mean-Field (HMF) model [9, 10]. 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 [11, 12, 13,
14]. 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
[15]. 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 [16, 17, 18, 19, 20].

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. II, 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. III
details of the initial condition for the simulations, the units, and the
simulation method. In Sec. IV, the numerical results and the comparison with
the theoretical framework are presented. Finally, Sec. V concludes this study.

##  II Analysis of the Vlasov equation with local fluctuations

We recall the Vlasov equation that governs the collisionless evolution of the
distribution function $f(\mathbf{x},\mathbf{v})$

$\frac{df}{dt}=\frac{\partial f}{\partial t}+\mathbf{v}\cdot\frac{\partial
f}{\partial\mathbf{x}}+\mathbf{a}\cdot\frac{\partial f}{\partial\mathbf{v}}=0$
(1)

where $\mathbf{a}$ 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 $\Phi$ as

$\mathbf{a}(\mathbf{x})=-\nabla_{\mathbf{x}}\Phi.$ (2)

By definition, $\Phi$ can be written in terms of the integral of $f$ as

$\Phi(\mathbf{x})=\int\int\varphi(\mathbf{x},\mathbf{x^{\prime}})f(\mathbf{x^{\prime}},\mathbf{v^{\prime}})d\mathbf{x^{\prime}}d\mathbf{v^{\prime}}$
(3)

where $\varphi(\mathbf{x},\mathbf{x^{\prime}})$ is the pair potential at
$\mathbf{x}$ from the mass element at $\mathbf{x^{\prime}}$. According to the
LB analysis [2], 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(\mathbf{x},\mathbf{v})$ denotes the point-wise one-particle distribution
function, Eq. (1) is known as the Klimontovich equation. In the case of
Newtonian gravity, Eq. (1) 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=\bar{f}+\delta f,$ (4)

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 [1] for review). The variation term vanishes by ensemble
average, i.e., $<\delta f>=0$, so $<f>$ converges to $\bar{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 $\delta\mathbf{a}$
accordingly, i.e., $\mathbf{a}=\mathbf{\bar{a}}+\delta\mathbf{a}$ where
$\mathbf{\bar{a}}$ is the mean-field acceleration from $\bar{f}$. From Eqs.
(3) and (4), we can deduce that $\delta\mathbf{a}$ is an integral of $\delta
f$ which vanishes after averaging in the same way as $\delta f$. All variation
terms are scaled by $1/\sqrt{N}$ as they originate from the Poissonian noise.
By substituting the perturbed quantities into Eq. (1) and performing the
ensemble average, the perturbed Vlasov equation becomes

$\bigg{<}\frac{d\bar{f}}{dt}\bigg{>}=-\bigg{<}\delta\mathbf{a}\cdot\frac{\partial\delta
f}{\partial\mathbf{v}}\bigg{>}\sim\frac{1}{N},$ (5)

where

$\frac{d\bar{f}}{dt}=\frac{\partial\bar{f}}{\partial
t}+\mathbf{v}\cdot\frac{\partial\bar{f}}{\partial\mathbf{x}}+\mathbf{\bar{a}}\cdot\frac{\partial\bar{f}}{\partial\mathbf{v}}.$
(6)

Although the ensemble average eliminates all terms involving first-order
fluctuations in Eq. (5), 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 [21, 22]. 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 $\delta 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 [23]. We have an alternative view on that term starting from
a simple dynamical perspective that $\delta\mathbf{a}$ which involves the
integral of $\delta f$ over the entire phase space, designates how large the
local fluctuation of mean-field acceleration is at any location. The product
$\delta f\delta\mathbf{a}$ in that term can then be regarded as the local
fluctuations of accelerations amplified by local variations $\delta f$, which
can also be interpreted as the strength of the deviation from mean-field
approximation. As the seed of this term, namely $\delta\mathbf{a}$, 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/\sqrt{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 $\bar{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\rightarrow\infty$. In that limit, the evolution follows
strictly the Vlasov equation, i.e., $df/dt=d\bar{f}/dt=0$. This happens when
one directly solves the Vlasov equation for the evolution of $f$, either
numerically or semi-analytically [24, 25, 26, 27, 28, 29]. 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. (5). We do not provide a full description of
how the perturbation (or $\delta 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 $\delta f$ can be assumed
to take a form of the standing wave. A recent work by [22] 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 $\delta 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.

##  III 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}=gm\frac{x_{j}-x_{i}}{|x_{j}-x_{i}|}$ (7)

where $g$ is the gravitational constant. The force expressed in Eq. (7) 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 [30]. 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}=\frac{M}{4x_{0}v_{0}}$ (8)

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}\equiv\frac{2T_{0}}{|U_{0}|}$, where $T_{0}$ and $U_{0}$ are the initial
kinetic and potential energies of the system, respectively, as

$f_{0}=\sqrt{\frac{M}{2gL^{3}b_{0}}}$ (9)

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}=\sqrt{\frac{L}{gM}}$ (10)

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.

##  IV Numerical results

###  IV.1 Waterbag evolution and phase-space density of QSS

Figure 1: Configurations in phase space for $b_{0}=0.1$ and $N=32000$ at
different times.

In this section, we consider first of all the evolution of the waterbag. Shown
in Fig. 1 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 [31] 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.

Figure 2: Ensemble-averaged phase-space density as a function of the energy
per unit mass $f(\epsilon)$ for $b_{0}=0.3$ (top panel) and $0.8$ (bottom
panel) and for different $N$. Plots are taken at $300\ t_{d}$. Horizontal
dashed line indicates the corresponding $f_{0}$. Size of error bars
corresponds to the standard deviation of the mean.

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. 2
is the ensemble-averaged phase-space density as a function of the energy per
unit mass $f(\epsilon)$ at $300\ t_{d}$ for $b_{0}=0.3$ and $0.8$, and for
various $N$. To compute $f(\epsilon)$, we first determine the energy
distribution of particles $F(\epsilon)$ such that $\int
F(\epsilon)d\epsilon=M$. Then, the conversion of the distribution from the
energy domain (or $F(\epsilon)$) to the phase-space domain (or $f(\epsilon)$)
is carried out by the relation

$F(\epsilon)=g(\epsilon)f(\epsilon)$ (11)

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

$g(\epsilon)=\int_{\Phi_{0}}^{\epsilon}\frac{2\sqrt{2}}{a(\Phi)\sqrt{\epsilon-\Phi}}d\Phi.$
(12)

In Eq. (12), $a(\Phi)=\frac{d\Phi}{dx}$ is the magnitude of the acceleration.
The lower integral limit $\Phi_{0}$ corresponds to the minimum potential
energy which is located at $x=0$. For simplicity, we can set $\Phi_{0}=0$
without loss of generality. We determine the conversion factor $g(\epsilon)$
by the numerical integration of Eq. (12). The ensemble-averaged $f(\epsilon)$
is obtained by averaging $f(\epsilon)$ 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(\epsilon)$ is tested by the bootstrapping method and the results
are reported in Appendix A.

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}$
[32, 33, 34]. 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 [8], 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. 1), $f(\epsilon)$ is well stationary in
the last $100\ t_{d}$, which is the time frame of consideration for further
analysis.

From the $f(\epsilon)$ plot, the QSS exhibits two distinct components: the
inner part has $f(\epsilon)$ close to $f_{0}$ while the outer part at higher
$\epsilon$ has much lower density. This form corresponds to the core-halo
structure as shown in Fig. 1 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 [35] for review). That theory was initially developed for explaining the
halo formation in an ion system [36], but the adaptation to self-gravitating
systems was also considered [37, 38, 39, 40]. 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 [39, 41]. The fact that the core ends up in the LB
degenerate state has been verified in past works on the self-gravitating sheet
model [42, 39, 20] and on other long-range interacting systems [38, 43], and
our frame of work is based on that fact. Because the filaments are still
intact at $\sim 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(\epsilon)=\begin{cases}f_{0},&\text{if $\Phi_{0}<\epsilon<\epsilon_{f}$}\\\
0,&\text{otherwise}\end{cases}$ (13)

where $\Phi_{0}$ is the minimum energy and $\epsilon_{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(\epsilon)$ can reach the value above $f_{0}$
and this excess is more noticeable for lower $N$. If $N$ is higher, the core
$f(\epsilon)$ is more uniform and closer to $f_{0}$, albeit with fluctuations,
and $f(\epsilon)$ falls more sharply to the halo level beyond the core energy
limit. This implies that the core is closer to the LB degenerate
$f(\epsilon)$. In fluid limit, the Vlasov equation suggests that the core
distribution function takes the degenerate form (13). 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. 2 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. (5). 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 [44].

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. 1, 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 [45, 46, 19] but this component was found to be meta-stable as it
disappeared when the system reached the thermal equilibrium [8].

###  IV.2 $N$-dependence of core phase-space density

Figure 3: Time evolution of the ensemble-averaged core phase-space density
$f_{c}$ for $b_{0}=0.3$ and for various $N$. Size of error bars is the average
of the mean errors in the core region.

In continuity with the observed excess of the phase-space density of the core
in Sec. IV.1, 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(\epsilon)$ 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. 3.
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. II.

The deviation of $f_{c}$ from $f_{0}$ can be investigated in a more
quantitative way by defining

$\Delta f\equiv\frac{f_{c}-f_{0}}{f_{0}}$ (14)

to be the phase-space density excess of the core and the plot of the magnitude
of the ensemble-averaged $\Delta f$ (or $|<\Delta f>|$) for different $b_{0}$
and $N$ is shown in Fig. 4. 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$, $|<\Delta 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 $|<\Delta 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^{\gamma}$ decay, where
$\gamma$ is a positive real number, for all $b_{0}$ and the results are
depicted in Fig. 4. 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 $\gamma$ 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. II. 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., [11, 12]). 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.

Figure 4: Magnitude of the averaged $\Delta f$ (or $|<\Delta f>|$) as a
function of $N$ for different $b_{0}$. The solid line is the best-fitting
$1/N^{\gamma}$ function with the best-fitting $\gamma$ reported for each case.
The dashed line is the $1/N$ decay for comparison. Size of error bars
corresponds to the averaged mean error in that time window.

###  IV.3 Core density profile and its agreement with LB theory

Figure 5: Density profile of the core $\rho(x)$ for $b_{0}=0.3$ (two top
panels) and $0.8$ (two bottom panels) in histogram. For each $b_{0}$, we
depict $\rho(x)$ for $N=1000$ and $32000$. Solid line is the averaged density
profile of the LB degenerate limit.

To inspect the core in more details, the ensemble-averaged density profile of
the core $\rho(x)$ is plotted in Fig. 5 for some selected $b_{0}$ and $N$. We
define the core to be the component below $\epsilon_{c}$ at which
$f(\epsilon)$ 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 (13) via the Poisson equation. The
plot of the LB density profile in Fig. 5 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 A. From the plot, the deviation from the LB
degenerate limit is remarkable when one inspects the density profile although
$|<\Delta 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(\epsilon)$
plots in Fig. 2. The plots in Fig. 5 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$ [19, 20]. 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 [47]. 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 [16,
17].

Figure 6: Density profile of the central component $\rho(x)$ below the energy
$\epsilon_{c,10}$ at which $f(\epsilon)$ falls to $0.1f_{0}$ for $b_{0}=0.3$
(left panel) and $0.5$ (right panel), both of which have $8000$ particles.
Solid line represents the density profile of the LB degenerate limit for the
corresponding $f_{0}$ and the mass of this component.

Finally, we revisit the ansatz of Teles et al (2011) [39] 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(\epsilon)$ is more complex than the three-level configuration (see Fig. 2).
We note the intermediate region where $f(\epsilon)$ decreases from the core
density to the halo density that spans a considerable range of $\epsilon$, 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(\epsilon)$ falls to $0.1f_{0}$, i.e.,
$\epsilon_{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 $\epsilon_{c,10}$ and the
corresponding LB degenerate profile are plotted in Fig. 6 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) [39]. 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 [39] 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.

##  V 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.

###### Acknowledgements.

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.

## References

  * Campa  _et al._ [2009] A. Campa, T. Dauxois, and S. Ruffo, Phys. Rep. 480, 57 (2009).
  * Lynden-Bell [1967] D. Lynden-Bell, Mon. Not. R. Astr. Soc. 136, 101 (1967).
  * Braun and Hepp [1977] W. Braun and K. Hepp, Commun. Math. Phys. 56, 101 (1977).
  * Chavanis [2010] P. H. Chavanis, J. Stat. Mech. Theory Exp. 5, 19 (2010).
  * Fouvry and Bar-Or [2018] J.-B. Fouvry and B. Bar-Or, Mon. Not. R. Astr. Soc. 481, 4566 (2018).
  * Miller [1996] B. N. Miller, Phys. Rev. E 53, R4279 (1996).
  * Tsuchiya  _et al._ [1996] T. Tsuchiya, N. Gouda, and T. Konishi, Phys. Rev. E 53, 2210 (1996).
  * Joyce and Worrakitpoonpon [2010] M. Joyce and T. Worrakitpoonpon, J. Stat. Mech. Theory Exp. 10, 10012 (2010).
  * Yamaguchi [2003] Y. Y. Yamaguchi, Phys. Rev. E 68, 066210 (2003).
  * Jain  _et al._ [2007] K. Jain, F. Bouchet, and D. Mukamel, J. Stat. Mech. Theory Exp. 2007, 11008 (2007).
  * Yamaguchi  _et al._ [2004] Y. Y. Yamaguchi, J. Barré, F. Bouchet, T. Dauxois, and S. Ruffo, Physica A 337, 36 (2004).
  * Antoniazzi  _et al._ [2007] A. Antoniazzi, F. Califano, D. Fanelli, and S. Ruffo, Phys. Rev. Lett. 98, 150602 (2007).
  * Campa and Chavanis [2017] A. Campa and P.-H. Chavanis, Journal of Statistical Mechanics: Theory and Experiment 5, 053202 (2017).
  * Santini  _et al._ [2022] A. Santini, G. Giachetti, and L. Casetti, J. Stat. Mech. Theory Exp. 2022, 013210 (2022).
  * Gabrielli  _et al._ [2014] A. Gabrielli, M. Joyce, and J. Morand, Phys. Rev. E 90, 062910 (2014), arXiv:1408.0999 [cond-mat.stat-mech] .
  * Yawn and Miller [1997] K. R. Yawn and B. N. Miller, Phys. Rev. Lett. 79, 3561 (1997).
  * Yawn and Miller [2003] K. R. Yawn and B. N. Miller, Phys. Rev. E 68, 056120 (2003).
  * Arad and Lynden-Bell [2005] I. Arad and D. Lynden-Bell, Mon. Not. R. Astr. Soc. 361, 385 (2005).
  * Yamaguchi [2008] Y. Y. Yamaguchi, Phys. Rev. E 78, 041114 (2008).
  * Joyce and Worrakitpoonpon [2011] M. Joyce and T. Worrakitpoonpon, Phys. Rev. E 84, 011139 (2011).
  * Benetti and Marcos [2017] F. P. C. Benetti and B. Marcos, Phys. Rev. E 95, 022111 (2017).
  * Chavanis [2022] P.-H. Chavanis, Physica A 606, 128089 (2022).
  * Chavanis [2008] P. H. Chavanis, Physica A 387, 787 (2008).
  * Janin [1971] G. Janin, Astron. Astrophys. 11, 188 (1971).
  * de Buyl and Gaspard [2011] P. de Buyl and P. Gaspard, Phys. Rev. E 84, 061139 (2011), arXiv:1108.0728 [cond-mat.stat-mech] .
  * Colombi and Touma [2014] S. Colombi and J. Touma, Mon. Not. R. Astr. Soc. 441, 2414 (2014).
  * Colombi [2015] S. Colombi, Mon. Not. R. Astr. Soc. 446, 2902 (2015).
  * Barnes and Ragan [2014] E. I. Barnes and R. J. Ragan, Mon. Not. R. Astr. Soc. 437, 2340 (2014).
  * Ragan and Barnes [2019] R. J. Ragan and E. I. Barnes, Mon. Not. R. Astr. Soc. 487, 3618 (2019).
  * Fanelli  _et al._ [2003] D. Fanelli, A. Noullez, and E. Aurell, J. Comp. Phys. 186, 697 (2003).
  * Miller  _et al._ [2023] B. Miller, G. Manfredi, D. Pirjol, and J.-L. Rouet, Class. Quantum Grav. 40, 073001 (2023).
  * Luwel  _et al._ [1984] M. Luwel, G. Severne, and P. J. Rousseeuw, Astrophys. Space Sci. 100, 261 (1984).
  * Mineau  _et al._ [1990] P. Mineau, M. R. Feix, and J. L. Rouet, Astron. Astrophys. 228, 344 (1990).
  * Rouet and Feix [1999] J. L. Rouet and M. R. Feix, Phys. Rev. E 59, 73 (1999).
  * Levin  _et al._ [2014] Y. Levin, R. Pakter, F. B. Rizzato, T. N. Teles, and F. P. C. Benetti, Phys. Rep. 535, 1 (2014).
  * Gluckstern [1994] R. L. Gluckstern, Phys. Rev. Lett. 73, 1247 (1994).
  * Levin  _et al._ [2008a] Y. Levin, R. Pakter, and T. N. Teles, Phys. Rev. Lett. 100, 040604 (2008a).
  * Levin  _et al._ [2008b] Y. Levin, R. Pakter, and F. B. Rizzato, Phys. Rev. E 78, 021130 (2008b).
  * Teles  _et al._ [2011] T. N. Teles, Y. Levin, and R. Pakter, Mon. Not. R. Astron. Soc. 417, L21 (2011).
  * Pakter and Levin [2011] R. Pakter and Y. Levin, Phys. Rev. Lett. 106 (2011).
  * Tashiro [2019] T. Tashiro, Physics Letters A 383, 125829 (2019), arXiv:1903.03307 [cond-mat.stat-mech] .
  * Lecar and Cohen [1971] M. Lecar and L. Cohen, Astrophys. Space Sci. 13, 397 (1971).
  * Teles  _et al._ [2010] T. N. Teles, Y. Levin, R. Pakter, and F. B. Rizzato, J. Stat. Mech. Theory Exp. 05, 05007 (2010).
  * Reidl and Miller [1995] C. Reidl and B. N. Miller, Phys. Rev. E 51, 884 (1995).
  * Hohl and Feix [1967] F. Hohl and M. Feix, Astrophys. J. 147, 147 (1967).
  * Yamashiro  _et al._ [1992] T. Yamashiro, N. Gouda, and M. Sakagami, Progress of Theoretical Physics 88, 269 (1992).
  * Chavanis [2006] P.-H. Chavanis, Physica A 365, 102 (2006).

##  Appendix A 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(\epsilon)$ (Fig. 7) and the density profile of the
core $\rho(x)$ (Fig. 8) for 5 different subsamples. The horizontal line
denoting $f_{0}$ is also provided in the $f(\epsilon)$ 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.

Figure 7: Phase-space density distribution $f(\epsilon)$ for $b_{0}=0.1$ and
$N=1000$ (left panel) and $b_{0}=0.1$ and $N=32000$ (right panel) for 5
different subsamples, each of which consists of randomly chosen $80\%$ of
realizations. The horizontal dashed line indicates the value of $f_{0}$.
Figure 8: Density profile of the core $\rho(x)$ for $b_{0}=0.1$ and $N=1000$
(left panel) and $b_{0}=0.1$ and $N=32000$ (right panel) for 5 different
subsamples, each of which consists of randomly chosen $80\%$ of realizations.