# Supervised Contrastive ResNet
and Transfer Learning for
the In-vehicle Intrusion Detection System

Thien-Nu Hoang, and Daehee Kim Authors are with the Department of Future
Convergence Technology, Soonchunhyang University, Asan, Republic of KoreaE-
mail<EMAIL_ADDRESS><EMAIL_ADDRESS>created October, 2020.

###### Abstract

High-end vehicles have been furnished with a number of electronic control
units (ECUs), which provide upgrading functions to enhance the driving
experience. The controller area network (CAN) is a well-known protocol that
connects these ECUs because of its modesty and efficiency. However, the CAN
bus is vulnerable to various types of attacks. Although the intrusion
detection system (IDS) is proposed to address the security problem of the CAN
bus, most previous studies only provide alerts when attacks occur without
knowing the specific type of attack. Moreover, an IDS is designed for a
specific car model due to diverse car manufacturers. In this study, we
proposed a novel deep learning model called supervised contrastive (SupCon)
ResNet, which can handle multiple attack identification on the CAN bus.
Furthermore, the model can be used to improve the performance of a limited-
size dataset using a transfer learning technique. The capability of the
proposed model is evaluated on two real car datasets. When tested with the car
hacking dataset, the experiment results show that the SupCon ResNet model
improves the overall false-negative rates of four types of attack by four
times on average, compared to other models. In addition, the model achieves
the highest F1 score at 0.9994 on the survival dataset by utilizing transfer
learning. Finally, the model can adapt to hardware constraints in terms of
memory size and running time.

###### Index Terms:

controller area network, intrusion detection, supervised contrastive learning,
transfer learning

##  I Introduction

The intelligent vehicle industry has gained considerable attention and
interest from companies, researchers, and consumers. Many electronic control
units (ECUs) are installed inside a smart vehicle to assist with more advanced
functions. These ECUs are interconnected through an in-vehicle network (IVN)
in which the controller area network (CAN) protocol is widely used, although
other technologies, such as the local interconnected network (LIN), FlexRay,
and Ethernet, are also available. The CAN protocol lacks encryption and
authentication mechanisms despite its fast speed and simplicity. Hence, there
is a trade-off between efficiency and security. Many studies demonstrate that
a CAN bus network can be attacked in various ways [1, 2, 3]. Thus, different
mechanisms have been introduced to fill the gap in the security of the CAN
bus. One of them is developing a system monitoring and detecting attacks in
the CAN bus network, which is called an intrusion detection system (IDS). With
the rapid development of the machine learning field, various studies applied
simple machine learning models, including K-means [4], K nearest neighbors
(KNN) [5], and one-class support vector machine (OSVM) [6], or used complex
deep learning models, such as convolutional neural network (CNN) [7, 8] and
long-short term memory (LSTM) [9, 10, 11] to build an efficient IDS. Although
the previous studies have achieved noticeable results with high detection
accuracy and low error rate, there are still two main issues. First, most of
them solve the problem of binary classification, which consists of two
classes: normal and abnormal. In this study, we aim to solve the multiclass
classification that distinguishes normal traffic and different types of
attacks, which will then notify the user of the specific type of attack. Thus,
we propose the concept of contrastive learning [12], which learns from
dissimilarity and similarity between training samples because multiclass
classification is more challenging than binary classification. In addition,
contrastive learning-based models can solve the challenge of class imbalance
in the CAN bus dataset. Second, previous machine learning-based IDS learns the
behavior of CAN message transmissions of a specific vehicle model. Further,
CAN ID meanings are distinct for different vehicle models although the CAN
messages follow the same structure. Therefore, it cannot be applied to other
vehicle models that require the collection of a new large dataset and training
of a new model if we want to develop a new efficient IDS for a newly launched
vehicle model. Hence, we apply the transfer learning technique to address the
problem in this study. The main contributions of this study can be summarized
as follows:

  * •

We introduce a novel model called supervised contrastive (SupCon) ResNet,
which combines supervised contrastive learning with the traditional ResNet and
minimizes the supervised contrastive loss, rather than the traditional cross-
entropy loss. To the best of our knowledge, our study is the first to apply
contrastive learning to the in-vehicle IDS.

  * •

We build a domain-adaptation transfer learning system, which facilitates the
generalization of a deep learning model across different CAN bus data from
different car models. We used the pre-trained SupCon ResNet model to transfer
learned knowledge from a rich source dataset to another target dataset. The
experiment results show that the proposed SupCon-based transfer learning
system lowers the false-negative rate of the training from scratch approach
and CE-based transfer learning by 0.4% and 0.2%, respectively.

  * •

We prove the efficiency of our proposed system through comprehensive
experiments and analyses on two popular real car datasets: the car hacking [7]
and survival [13] datasets. The results show that our proposed model is
superior to other baselines in terms of detection performance (obtained the
highest F1 score of 0.9998), memory size, and running time.

After providing an overview of the problem and stating the contributions of
the study in Section I, the relevant background knowledge and security issues
in the CAN bus protocol are introduced in Section II. In Section III, we
summarize various studies related to the in-vehicle IDS. Our proposed system
is described in detail in Section IV. Next, the experimental setup and results
are presented in Section V, where we showed the efficiency of the proposed
method empirically. Finally, Section VI presents the conclusions and
directions for future works.

##  II Can Bus Background

Figure 1: CAN2.0B data frame format.

###  II-A CAN Bus System

The controller area network (CAN) is one of the protocols used in the in-
vehicle network comprising many connected ECUs. The CAN bus, which is invented
by Robert Bosh GmbH in the early 1980s, has become common in automotive
systems because of its advantages, including high bit rates (up to 1 Mbit/s),
cost-effectiveness, and system flexibility [14]. The concept of CAN protocol
comprises message prioritization and multi-master. Every CAN message contains
an ID that indicates the content of the message. Based on this, each ECU
obtains its relevant messages from a message filterer. There is no source and
destination information in CAN messages. In addition, the ID implies the
priority of a message: the lower ID, the higher priority. When the bus is
idle, any unit can start to transmit a message. The message whose higher
priority wins the bus and takes the right of sending the message. During
transmission, any number of nodes can receive and simultaneously act upon the
same message according to the concept of message filtering.

A CAN message is called a frame belonging to four types, such as data, remote,
error, and overload frames. While the data frame contains the main information
in ECU communication, the remote frame is used when a unit requests peculiar
information from other ECUs. Meanwhile, the error frame supports error
detection and fault confinement. An overload frame will be sent whenever a
node requires a delay of the next frame. Every CAN frame is made up of a
sequence of dominant bits (0) and recessive bits (1). The data and remote
frames have the same structure, as shown in Fig. 1. The data and remote frames
start with a bit called the start of frame (SOF), which is always a dominant
bit (0). According to the length of the ID, the CAN data are classified into
CAN 2.0A and CAN 2.0B. CAN 2.0A is a shortened version of CAN 2.0B, without
the extended identifier. Hence, the CAN 2.0B format structure is described in
Fig. 1. The 11-bit length base identifier is followed by a 1-bit substitute
remote request (SRR) and 1-bit identifier extension (IDE). The successive 18
and 6 bits represent the extended identifier and control field, respectively.
The data field containing 64 bits is the content of the message. The following
16 bits for cyclic redundancy check (CRC) and 2 bits for acknowledgment are
utilized in error detection and correction. The end of the frame (EOF) is
recognized by seven consecutive bits. For more information regarding the other
frame types, we refer to the CAN specification reported by Bosh [14].

###  II-B Security in CAN Bus Protocol

Despite the advantages mentioned above, the most severe problem in CAN
protocol is message transmission without any authentication and encryption.
Particularly, any node connected to the bus can obtain the message during
message transmission. Consequently, an attacker who compromises a node on the
bus can easily sniff the information on the bus. There are two levels of an
adversary: weak and strong [15]. A weak attacker can prevent the ECU from
transmitting messages or keep the ECU in a listen-only mode. Meanwhile, a
strong attacker can completely control an ECU and have access to the memory
data. Consequently, malicious messages can be injected to launch the attack in
addition to the abilities of the weak attacker. These two levels of attacks
are proved in several studies [1, 2, 3]. In this study, our proposed IDS
handles the strong attacker including three different types of injected
messages:

  * •

DoS: The attacker floods the bus by injecting messages containing the highest
priority ID 0x0000 and arbitrary data fields. The legitimate messages are
prevented from being transmitted, resulting in unusual effects, such as
flashing dash indicators, intermittent accelerator/steering control, and even
full vehicle shutdown, since the message with the ID 0x0000 always wins the
bus [15].

  * •

Fuzzy: The attacker injects messages with an arbitrary ID and data fields at a
high frequency. Aside from being randomly generated, the injected IDs can be
chosen from those appearing in the normal traffic, which is supposed to be
difficultly detected. The effect of this attack is similar to the DoS attack.

  * •

Target ID: The attacker injects messages with a specific ID and manipulated
data fields as his intent. To achieve this, the attack requires the attacker
to have knowledge of the meaning of the specific ID by reverse engineering
techniques.

##  III Related Works

###  III-A Deep Learning-based In-vehicle IDS

Two main approaches to counter in-vehicle cyber-attacks are message
authentication and intrusion detection [15]. Message authentication consumes
more resources and time, although it provides a higher security level.
Consequently, these approaches decrease the CAN bus performance. Therefore,
intrusion detection is preferred. Different types of IDS have been proposed to
monitor and analyze transferred messages on the bus, detect vicious behaviors,
and make an alert if detected.

With the rapid development of deep learning models, many deep learning-based
in-vehicle IDS have been introduced. In 2016, a previous study in [16] was one
of the first research that applied deep learning for CAN bus intrusion
detection. The authors developed a deep neural network where each neuron takes
each bit of data payload in the CAN message as the input. The network was
trained in a supervised manner with two classes: normal and abnormal. However,
the proposed model was investigated using a self-synthesized dataset with
simple attack models. Although the result was not impressive, this study is
the first to apply advanced deep learning models to the CAN bus network.
Several studies applied recurrent neural networks (RNN) to identify injected
messages since CAN messages possess sequential patterns. In [9], the authors
used the long short term memory (LSTM), which is the improvement of RNN, to
predict the next bit in the data payload of normal CAN message sequences. The
invasion is detected by checking the difference between the predicted and
receiving values on the bus. Meanwhile, the LSTM-based autoencoder, which was
adopted by [10], was trained with special statistic features of the CAN
message sequences. These models were designed for a message sequence of a
specific CAN ID. This implies that we need to train many models separately,
corresponding to the number of IDs in a CAN network of a vehicle. To reduce
the number of models, the authors in [11] employed the bi-directional
generative pretrained transformer (GPT) network to predict the next CAN ID in
the CAN message sequence, which was used for abnormal detection to identify
attacks.

Aside from RNN, the convolutional neural network (CNN) was exploited in many
CAN bus IDS studies. CNN is commonly used for various computer vision tasks,
such as object detection, image colorization, and image segmentation. Hence, a
CNN-based IDS requires an image representation of CAN messages. In [17],
multiple one hot vectors representing consecutive CAN IDs were stacked to form
a CAN image, which is fed into the CNN generative adversarial networks (GAN).
The proposed method can perform binary classification (normal/abnormal) with
an overall accuracy of 97% and detect unknown attacks. Similarly, a previous
study [7] manipulated CAN ID sequences in binary form, which is used to train
a simplified Inception Resnet. The proposed model achieves a low false-
negative rate in the case of DoS and spoofing attacks, but the result of the
fuzzy attack is not good. Meanwhile, the authors in [8] proposed a concept of
recurrence plots, which is the matrix of subtractions value of multiple CAN
IDs within a specific window size. In addition, they suggested using a
lightweight CNN model with the proposed recurrence plots to foster the speed
of IDS. However, the accuracy of multiclass classification is not high.

###  III-B Transfer Learning in the In-vehicle IDS

All previous models can only be applied to a specific car model. This
indicates retraining a new model on a new dataset is required if we build an
IDS for a newly released car. In reality, it is time-consuming and requires
lots of effort to collect data for every newly released car model. Therefore,
transfer learning is proposed to address this problem. As far as we know, the
study [18] is the first study that employed transfer learning for CAN bus data
and presented an LSTM-based model to solve the binary classification. The
authors tested their proposed scheme on the survival dataset [13], which is a
small size dataset. The data of Kia Soul and Chevrolet Spark were used for
testing the transfer learning, while the data of Hyundai Sonata car was
employed for training the proposed model. However, transfer learning should be
used when there is a large source dataset, and the source and target tasks
should have the same domain. Therefore, the proposed scheme did not achieve
good results.

###  III-C Contrastive Learning in IDS

To the best of our knowledge, there are no studies applying contrastive
learning in the in-vehicle IDS design. However, some recent works have
proposed the use of contrastive learning to the network IDS [19, 20, 21]. For
example, [19] combined the autoencoder and triplet loss to demonstrate an IDS
on various network datasets, such as
KDDCUP99111http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html,
AAGM17222https://www.unb.ca/cic/datasets/android-adware.html, and
CICIDS17333https://www.unb.ca/cic/datasets/ids-2017.html. Meanwhile,
autoencoder-based contrastive learning was introduced as a part of a multi-
task model in [20]. In addition, [21] suggested a novel concept of contrastive
learning, in which they projected the labels and features into the same
representation space. In the classification phase, the predicted label is a
class having the closest distance to the input features in that representation
space. All studies mentioned suggested that contrastive learning is suitable
to solve the problem of class imbalance, which occurs in IDS research.

##  IV Methodology

Figure 2: Overview of the proposed system.

###  IV-A Problem Formulation

The intrusion detection system is installed in the CAN bus to monitor the CAN
network and provide an alert if there is any malicious message. By using the
fluctuation in the CAN ID sequence of consecutive CAN messages, this study
aims to build an intrusion detection system $f$ that classifies a sequence of
CAN IDs into a set of $C$ attack types, denoted as
$\mathcal{A}=\\{a_{0},a_{1},a_{2},..,a_{C}\\}$, where $a_{0}$ is a normal
sequence. Suppose that several CAN log messages are collected from a car model
under both normal and attack circumstances. After processing, the data
$D=\\{(X_{i},y_{i})\\}$ are ready to train a machine learning model for the
intrusion detection task, where $X_{i}$ is a sequence of CAN IDs and
$y_{i}\in\mathcal{A}$ is the corresponding label. The objective is to build a
multiple class classification IDS $f(\cdot|X)$, which determines whether a CAN
IDs sequence $X_{i}$ belongs to a normal class or one of $C$ predefined attack
types as accurately as possible.

The most important challenge for a supervised-based IDS is how to collect
enough attack samples, especially for a newly released car model. We cannot
apply a model trained on a CAN dataset A to another CAN dataset B, if A and B
are from different manufacturers because each car model owns a specific CAN
IDs pattern. Transfer learning can address this problem by extracting
meaningful knowledge from a pre-trained model. We assume that an elegant
source dataset $D_{s}$ exists collected from a car model $s$. From this, we
can build an efficient IDS $f_{s}(\cdot|X_{s})$ for the source car. The final
objective is to build an IDS model for the target model
$f_{t}(\cdot|X_{t},f_{s})$, with a limited target dataset $D_{t}$ where
$|D_{t}|\ll|D_{s}|$.

###  IV-B Proposed System

The proposed system presented in Fig. 2 comprises two stages: supervised
contrastive learning from the source data and transfer learning for the target
data. While supervised contrastive learning reduces the effect of class
imbalance in the training data, transfer learning addresses the problem of
limited data on a newly released car model.

In the first stage, the source data are fed into the SupCon ResNet model after
being preprocessed. The proposed SupCon ResNet model includes three parts: the
encoder, projection, and classification networks.

  * •

The encoder network (Encoder) maps the input $X$ into the representation
feature space. In this study, we designed a compact ResNet-18 architecture
used as the encoder network.

  * •

The projection network (Projector) projects the outputs from the encoder
network into another space, which is normalized and used to measure the
distance for the supervised contrastive loss. From the experimental results,
the model with a projector produces a higher accuracy, compared to the model
without it [22]. In addition, a non-linear projection network is assumed to be
a better choice than the linear one. Hence, a multiple layer perception (MLP)
with one hidden layer is used in our model.

  * •

The classification network (Classifier) is a linear layer performing the final
classification task that maps the outputs from the encoder network to $C$
classes.

The classifier is trained with traditional CE-entropy loss thereafter, while
the encoder and the project networks are trained together with SupCon loss.

In the second stage, we performed data processing on the target data the same
as the source data. Training the target model from scratch may lead to the
overfitting problem because of the small size dataset. Hence, we utilized the
pre-trained model $f(s)$ for fine-tuning the target model. After completing
the training of the source model, we transferred the weights of the encoder
from the source model to the target model. After training the classifier on
top of the frozen encoder, the entire model is trained with a small learning
rate and the unfrozen encoder. The outcome of the system is the IDS model
composed of the encoder and classification networks, which will be deployed in
a real environment.

###  IV-C Data Preprocessing

In this study, only CAN IDs in the CAN messages were used, which was motivated
by the message transmission mechanism in the CAN bus protocol [7]. Each CAN ID
in the CAN bus has its message transmission cycle. Meanwhile, the message is
broadcasted on the CAN bus based on the priority mechanism. Therefore, we
believe that a certain sequential pattern exists in the CAN ID sequence. If an
attacker injects some messages into the CAN bus, the pattern will be broken.
We stacked $N$ sequential CAN IDs representing 29-bits, to create a CAN ID
sequence frame. We choose $N=29$ because of the convenience of processing
square matrix in the CNN architecture network. Consequently, the input fed
into the proposed model is a matrix of $29\times 29$, which is visualized in
Fig. 5.

###  IV-D ResNet Architecture

Figure 3: The residual block in ResNet

The ResNet (Residual Network) [23], which is introduced first in 2015 by He et
al., has become a popular architecture in deep learning. This allows us to
train a much deeper network without performance degradation. Theoretically,
the more layer, the better result a deep learning model provides. However, in
reality, training a deep network is challenging because of vanishing gradient
and saturated training error. In a previous study regarding ResNet, the
authors introduced the concept of skip connection to solve the problems. At
that time, the 152-layer ResNet architecture outperformed the other well-known
models, such as AlexNet and VGG, to win first place in the 2015-ILSVRC
competition.

The Fig. 3 displays a residual block, the core element in the ResNet. The
block takes the input $\mathbf{x}$ and outputs $\mathbf{y}$ as follows:

$\mathbf{y}=\mathcal{F}(\mathbf{x},\\{W_{i}\\})+\mathbf{x}$ (1)

where the function $\mathcal{F}(\mathbf{x},\\{W_{i}\\})$ represents the
residual mapping to be learned with the weight of layers $W_{i}$. The
$\mathbf{x}$ identity mapping is called a “shortcut connection” that is
element-wise addition to the function $\mathcal{F}$. Consequently, the simple
concept of adding the identity $\mathbf{x}$ to the outcome does not increase
the total number of parameters, depth, width, and computational cost (except
for the negligible element-wise addition), compared to the plain network. The
authors also showed that ResNet, which was 20 and 8 times deeper than AlexNet
and VGG respectively, still has lower complexity.

The CAN bus IDS requires not only high accuracy but also low running time.
Therefore, we design a compact ResNet-18 to meet these requirements. The Fig.
4 shows the detail of our architecture used for the encoder network.
Specifically, we reduced the number of channels of each layer in the original
ResNet-18 architecture by four times to produce the compact version. After
passing through the convolutional layer with a kernel size of three, the
sample goes through eight residual blocks to produce a hidden features
representation with the size of 128.

Figure 4: Detail of the compact ResNet encoder network

###  IV-E Supervised Contrastive Model

Figure 5: Difference between the SupCon ResNet and traditional CE ResNet
models

Compared to the traditional cross entropy that learns from the label,
contrastive learning learns from the dissimilarity between samples. In
general, the goal of contrastive learning is to produce embedding features
wherein similar samples stay close to each other, while dissimilar ones are
distant from each other. The original contrastive loss attempts to minimize
the embedding distance when they are from the same class by taking a pair of
$(\mathbf{x}_{i},\mathbf{x}_{j})$ as the input; otherwise, it maximizes the
distance:

$\mathcal{L}_{\text{cont}}(\mathbf{x}_{i},\mathbf{x}_{j},\theta)=\mathbf{1}[y_{i}=y_{j}]||f_{\theta}(\mathbf{x}_{i})-f_{\theta}(\mathbf{x}_{j})||^{2}_{2}\\\
+\mathbf{1}[y_{i}\neq
y_{j}]\text{max}(0,\epsilon-||f_{\theta}(\mathbf{x}_{i})-f_{\theta}(\mathbf{x}_{j})||^{2}_{2})$
(2)

where $y_{i}$, $y_{j}$ are the labels of $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$
respectively; $f_{\theta}(.)$ is the encoder function; and $\epsilon$ is a
hyperparameter defining the lower bound distance between samples of different
classes. Contrastive loss can be better than the cross-entropy if it is
trained with hard positive/negative pairs. For example, a hard positive pair
are two samples that belong to the same class but appear different. Meanwhile,
a negative pair are those from different classes but appear quite similar. The
challenge of how to design an efficient hard mining technique for training
contrastive loss has still been being developed. In 2020, the SupCon loss [24]
was proposed to solve the problem of the original contrastive loss. It
outperformed the traditional cross-entropy on the ImageNet datasets in terms
of the accuracy. The SupCon loss is extended from the contrastive learning in
a self-supervised manner [22] by leveraging the label information. Fig. 5
illustrates the structure of the supervised contrastive model, compared to
conventional neural networks trained with cross-entropy loss. The SupCon
architecture has three main parts: the encoder network, projector network, and
classifier. The encoder network maps the input $\mathbf{x}$ to a
representation vector, $\mathbf{r}=Enc(\mathbf{x})$. After being normalized,
$\mathbf{r}$ is fed into the projector network and produces the
$\mathbf{z}=Proj(\mathbf{r})$. The encoder and projector network are trained
with the supervised contrastive loss as below:

$\mathcal{L}_{\text{supcon}}=-\sum_{i=1}^{N}\frac{1}{N_{y_{i}}-1}\times\\\
\Bigg{[}\sum_{j=1}^{N}\mathbf{1}_{i\neq
j}\mathbf{1}_{y_{i}=y_{j}}\text{log}\frac{\text{exp}(\mathbf{z}_{i}\cdot\mathbf{z}_{j}/\tau)}{\sum_{k=1}^{N}\mathbf{1}_{k\neq
i}\text{exp}(\mathbf{z}_{i}\cdot\mathbf{z}_{k}/\tau)}\Bigg{]}$ (3)

where $N$ is the total number of training samples, $N_{y_{i}}$ is the number
of positive samples that have the same label $y_{i}$ with the sample $i$, $j$
and $k$ are the indices of positive samples and all samples in the training
set, respectively. In addition, temperature parameter $\tau$ controls the
smoothness of probability distribution. A small $\tau$ will be good for
training, but too small $\tau$ can lead to unstable training because of
numerical instability. In the original paper, the authors state that
$\mathcal{L}_{supcon}$ possesses the implicit property that encourages the
hard positive/negative mining without performing it explicitly. As a result,
the learning process for the supervised contrastive model does not need to
perform hard mining, which slows down the training time. For more detail, we
refer to the proof of the original paper [24]. By learning from
dissimilarity/similarity, the model produces good representations
$\mathbf{r}$, which are then used to train the classifier with cross-entropy
loss as usual.

###  IV-F Transfer Learning

Transfer learning is a technique for improving the performance of a target
model on a limited target data by utilizing knowledge from a related source
model, trained on an abundant source data. According to [25], the formal
definition of transfer learning is described as follows:

  * •

Domain: A domain $\mathcal{D}=\\{\mathcal{X},P(X)\\}$ includes two components:
a feature space $\mathcal{X}$ and a marginal distribution $P(X)$ of data $X$.

  * •

Task: A task $\mathcal{T}=\\{\mathcal{Y},f\\}$ comprises a label space
$\mathcal{Y}$ and a predictive function $f=P(Y|X)$, which is trained on a
dataset consists of pairs $\\{x_{i},y_{i}\\}$, where $x_{i}\in X$ and
$y_{i}\in Y$.

  * •

Transfer learning: Given a source domain $\mathcal{D}_{s}$, a corresponding
source task $\mathcal{T}_{s}$, and a target domain $\mathcal{D}_{T}$ and
target task $\mathcal{T}_{T}$, the objective of transfer learning is to learn
the target conditional probability distribution $P(Y_{T}|X_{T})$ in
$\mathcal{D}_{T}$ with the information gained from $\mathcal{D}_{S}$ and
$\mathcal{T}_{S}$. Here, $\mathcal{D}_{S}\neq\mathcal{D}_{T}$ or
$\mathcal{T}_{S}\neq\mathcal{T}_{T}$. In most cases, a sufficient number of
labeled source examples, which is extremely larger than the number of labeled
target examples, are assumed to be available.

Based on the relationship between the source and target domains, source and
target tasks, and availability of data, transfer learning techniques have
various types: inductive, transductive, and unsupervised. For more details on
these categories, we refer to the survey of transfer learning [25].

We have the $\mathcal{D}_{S}$ and $\mathcal{D}_{T}$ domains from different car
models to apply the transfer learning definition to the IDS for the CAN bus.
Although they have the same data format that follows the predefined CAN
messages structure, each dataset owns a different message transmission
behavior. This implies the $\mathcal{X}_{S}=\mathcal{X}_{T}$, whereas
$P(X_{S})\neq P(X_{T})$. These properties prevent the ability of a generalized
IDS model for all car models. In practice, it is difficult to collect a large
amount of training data for each car model. Based on the framework of [25],
the domain-adaptation inductive transfer learning technique is suitable for
our problem. We utilized the encoder in the supervised contrastive learning
model, which is trained on a source data, as a pre-trained model. The
predictive target model includes the encoder that has the same structure as
the source model and top classifier, which is attached to solve the target
task. First, we trained the target model with the frozen encoder. This step
can be considered as the initialization of the weights of the classifiers.
Finally, we trained the entire model after unfreezing the encoder network with
a small learning rate to avoid overfitting.

##  V Experimental Results

###  V-A Experiment Setup

TABLE I: Number of messages in car hacking and survival datasets.

Attack type |  #Messages (Sonata) |  #Messages (Soul) |  #Messages (Spark)
---|---|---|---
DoS Attack |  3,078,250 587,521 |  181,901 33,141 |  120,570 22,587
Fuzzy Attack |  3,347,013 491,847 |  249,990 39,812 |  65,665 5,812
Malfunctioning |  - |  173,436 7,401 |  79,787 8,047
Gear Spoofing |  4,443,142 597,252 |  - |  -
RPM Spoofing |  4,621,702 654,897 |  - |  -

  * •

Note: The Sonata values come from the car hacking dataset, while the other
values are from the survival dataset. The red color indicates the injected
messages, while the black indicates the normal messages.

We used two popular datasets produced by the Hacking and Countermeasure
Research Lab (HCRL) of Korea University: the car hacking [7] and survival [13]
datasets. The car hacking dataset was collected from the Hyundai Sonata model,
and the survival dataset was collected from three specific car models (i.e.,
Hyundai Sonata, KIA Soul, and Chevrolet Spark). We used the data from the
Sonata Hyundai car model in the car hacking dataset as the source and treated
the others as the target ones since the data of the Hyundai Sonata car model
is much larger than the others. The label space of source and target tasks are
different: $\mathcal{Y}_{S}\in\mathbb{R}^{5}$ including normal, DoS, fuzzy,
spoofing RPM, and spoofing gear information. Meanwhile,
$\mathcal{Y}_{T}\in\mathbb{R}^{4}$, including normal, DoS, fuzzy, and
malfunction. Table I lists the details of the data size.

Each message in both two datasets follows the same format containing a
timestamp, CAN ID in HEX, the number of data bytes, 8 bytes of data, and the
label. We extracted CAN IDs and transformed them from hexadecimal to 29-bit
representation. Then, the CAN ID sequence frame was built by stacking 29
sequential samples together. A stride value $s$ is the number of messages
between two continuous frames. More frames are created if a stride value is a
small but less variant between these frames. We choose $s=15$ for the car
hacking dataset and $s=10$ for the survival dataset because the size of the
survival dataset is quite smaller than that of the car hacking dataset.

The dataset comprises multiple files corresponding to the attack types. For
each file, the frame was labeled as normal if there is no injected message,
whereas the label of the frame is a non-zero integer number considering the
type of attack. We combined all of the processed frames to obtain the final
dataset. The train/test splitting rate is 7:3. We trained the model on the
training set and evaluated it on the test set. For each model, we performed
the entire process from dataset splitting to evaluation five times and
reported the average results.

The experiments were conducted on a server provided with 32 Intel(R) Xeon(R)
Silver 4108 CPUs @ 1.80 GHz, a memory of 128 GB, and an Nvidia Titan RTX 24GB
GPU. All of the code written in Python 3.7.10 and Pytorch 1.9.0 is published
on the Github444Source code is available at https://github.com/htn274/CAN-
SupCon-IDS.

###  V-B Evaluation Metrics

As our problem is multiclass classification, we evaluated the proposed model
using the false-negative rate (FNR), recall (Rec), precision (Prec), and F1
score (F1). These metrics were calculated from the confusion matrix using
true/false positive/negative samples. When we calculated a metric of a
specific class, this class is considered positive, whereas the others are
considered negative. Concretely, the FNR is the fraction between the false-
negative and the total number of positive.

$\text{FNR}=\frac{FN}{TP+FN}$ (4)

In the case of the IDS, the FNR must be as small as possible. In addition, the
FNRs of the attack classes are expected to be smaller than of the normal class
because the consequence of a missed detected attack sample is more dangerous
than that of a false alarm alert. Meanwhile, the recall is calculated as the
fraction between true positive and the total sum of true positive and false
negative.

$\text{Rec}=\frac{TP}{TP+FN}$ (5)

The precision is calculated as the fraction between the true positive and the
total sum of true and false positives.

$\text{Pre}=\frac{TP}{TP+FP}$ (6)

For our problem, the proposed model should achieve both high recall and high
precision. Therefore, the F1 score is used to balance these two metrics. The
F1 score is defined as the harmonic mean of precision and recall.

$\text{F1}=2\times\frac{Prec\times Rec}{Prec+Recall}.$ (7)

###  V-C Choosing Hyperparameters for the SupCon ResNet Model

Choosing the right hyperparameters is important in deep learning because it
can boost the performance of a model. Three hyperparameters should be
considered in the SupCon ResNet model: the learning rate, batch size, and
$\tau$ value. For simplicity, we set the $\tau$ value as 0.07 as advised in
the original paper [24] and tuned only the learning rate and batch size. The
authors of the SupCon study stated that the SupCon loss takes advantage of
large batch sizes. Hence, we tried different batch sizes of 512, 1024, 2048,
and 4096. We set the learning rate as 0.05 for 512 and 1024 batch sizes.
Meanwhile, the learning rate of 0.1 was used for the model with the batch size
of 2048 and 4096. Finally, we trained the top classifier for four models with
the same configurations with 256 batch sizes and a learning rate of 0.01 after
150 epochs. We reported the false negative errors grouped by attack types of
four models in Fig. 6. The results showed that the classifier trained with
SupCon loss classifies normal and attack samples well on the Car Hacking
dataset, as all models have false-negative errors lower than 0.1%. The Fuzzy
attack has the highest false-negative errors in almost settings, while the DoS
attack is easily detected. This finding is aligned with the study in [7].
Moreover, we compared these models in terms of the average false-negative
errors. The results suggest that the larger batch size does not show better
results. Overall, the model with a batch size of 512, which is the best model
with the lowest mean of false-negative errors, is chosen for the following
experiments.

Figure 6: False negative errors grouped by attack types of each set of
hyperparameters

###  V-D Comparison with other Supervised Methods

We compared the proposed SupCon ResNet model to two other models: the
Inception ResNet [7], which is the first study that demonstrates the concept
of using CNN for CAN bus messages data, the baseline with the same compact
ResNet encoder but using the cross-entropy loss. Regarding the Inception
ResNet, we reimplemented the model with the hyperparameters reported in the
study. Meanwhile, the ResNet trained with the cross-entropy loss was trained
with a batch size of 256 and a learning rate of 0.001. Table II summarizes the
results of these models, which are the average taken from five times running.

Although the CE ResNet obtains the highest false-negative error of the DoS
attack, at 0.21%, those of the other attacks are slightly lower than those of
the Inception ResNet, by 0.06%, 0.08%, 0.09% for the fuzzy attack, gear
spoofing attack, and RPM spoofing attack, respectively. Surprisingly, the
results of the SupCon ResNet outperform the others. Based on Table II, the
SupCon ResNet reduces the false-negative errors of all attacks almost six
times on average. In addition, the proposed model achieves the highest F1
score for all of the attacks, roughly 0.9998 on average.

Fig. 7 illustrates the confusion matrices of the two models to show the
superior of the SupCon ResNet compared to the baseline model. By learning from
(dis)similarity, the SupCon ResNet decreases more than half of the number of
miss-detected attack samples. There are 51 DoS samples classified falsely as a
normal class in the case of the CE one, whereas there is only 1 sample in the
case of the SupCon ResNet. Moreover, the SupCon ResNet also slightly decreases
the false-negative samples within the attack classes and the number of false
alarm cases.

In summary, these results indicate that the SupCon ResNet model separates the
classes well on the car hacking dataset. This is because of the motivation of
contrastive learning: the similar samples are pulled together while the
dissimilar samples are pushed far apart. Due to the richness of the car
hacking dataset, the SupCon ResNet can learn good representations, which are
useful for classification, not only on the car hacking dataset but also on
other CAN datasets, which will be proved empirically in the next section.

(a)

(b)

Figure 7: Comparison of confusion matrix between CE and SupCon models TABLE
II: Comparison of the proposed SupCon model
and other supervised models Attack Type | Model | FNR | Rec | Prec | F1
---|---|---|---|---|---
DoS | Inception ResNet | 0.12% | 0.9988 | 0.9998 | 0.9993
| CE ResNet | 0.21% | 0.9979 | 0.9997 | 0.9987
| SupCon ResNet | 0.02% | 0.9998 | 1.0 | 0.9999
Fuzzy | Inception ResNet | 0.30% | 0.9970 | 0.9982 | 0.9976
| CE ResNet | 0.24% | 0.9976 | 0.9997 | 0.9987
| SupCon ResNet | 0.06% | 0.9994 | 1.0 | 0.9997
Gear Spoofing | Inception ResNet | 0.19% | 0.9981 | 0.9996 | 0.9988
| CE ResNet | 0.11% | 0.9989 | 0.9997 | 0.9993
| SupCon ResNet | 0.06% | 0.9994 | 0.9999 | 0.9997
RPM Spoofing | Inception ResNet | 0.21% | 0.9979 | 0.9998 | 0.9988
| CE ResNet | 0.12% | 0.9988 | 0.9998 | 0.9993
| SupCon ResNet | 0.06% | 0.9994 | 1.0 | 0.9997

###  V-E Transfer Learning Results

The SupCon ResNet produces good representations that contribute learned
knowledge to train an efficient classification model using a small dataset. In
this section, we will analyze the results of transfer learning from the
proposed model, trained with the car hacking dataset, to other smaller-size
and different car model datasets. We set up the experiment with three
different configurations:

  1. 1.

Random: We trained the model on the new dataset from scratch. This means that
the weight of the model is randomly initialized.

  2. 2.

CE ResNet: We used the CE ResNet trained on the car hacking dataset as the
pretrained model.

  3. 3.

SupCon ResNet: We used the SupCon ResNet trained on the car hacking dataset as
the pretrained model.

The results of KIA Soul and Chevrolet Spark models are presented in Table III
and Table IV, respectively.

Overall, transfer learning reduces false-negative errors and increases the F1
score, especially for the fuzzy attack. For example, on both two target
datasets, the false-negative errors of the fuzzy attack reduce by two times
when using the CE ResNet as a pretrained model, compared to the random model.
The number decreases significantly to 0.06% and 0.47% in the case of KIA Soul
and Chevrolet Spark models, respectively, when the CE ResNet is replaced by
the SupCon ResNet model. From the experiment results, we observe two important
conclusions as follows:

  * •

For the IDS CAN domain, in the case of a small size dataset, transfer learning
from a different car model can increase the detection abilities in terms of
false-negative errors and F1 scores.

  * •

Knowledge obtained from the SupCon ResNet model is better than that of the CE
ResNet model. In fact, it increases the overall F1 score.

TABLE III: Transfer learning results for Kia Soul data Attack Type | Model | FNR | Rec | Prec | F1
---|---|---|---|---|---
Normal | Random | 0.04% | 0.9996 | 0.9980 | 0.9988
CE ResNet | 0.04% | 0.9996 | 0.9990 | 0.9993
SupCon ResNet | 0.0% | 1.0 | 0.9999 | 0.9999
DoS | Random | 0.12% | 0.9988 | 1.0 | 0.9994
CE ResNet | 0.05% | 0.9995 | 1.0 | 0.9997
SupCon ResNet | 0.01% | 0.9999 | 1.0 | 1.0
Fuzzy | Random | 0.49% | 0.9951 | 0.9994 | 0.9973
CE ResNet | 0.20% | 0.9980 | 0.9992 | 0.9986
SupCon ResNet | 0.06% | 0.9995 | 0.9999 | 0.9997
Malfunction | Random | 0.20% | 0.9980 | 0.9973 | 0.9977
CE ResNet | 0.21% | 0.9979 | 0.9984 | 0.9981
SupCon ResNet | 0.06% | 0.9995 | 0.9999 | 0.9997
Overall | Random | 0.21% | 0.9979 | 0.9987 | 0.9983
CE ResNet | 0.13% | 0.9988 | 0.9992 | 0.9989
SupCon ResNet | 0.03% | 0.9997 | 1.0 | 0.9998
TABLE IV: Transfer learning results for Chevrolet Spark data Attack Type | Model | FNR | Rec | Prec | F1
---|---|---|---|---|---
Normal | Random | 0.25% | 0.9975 | 0.9960 | 0.9967
CE ResNet | 0.11% | 0.9989 | 0.9982 | 0.9985
SupCon ResNet | 0.11% | 0.9989 | 0.9992 | 0.9991
DoS | Random | 0.18% | 0.9982 | 1.0 | 0.9991
CE ResNet | 0.05% | 0.9995 | 1.0 | 0.9998
SupCon ResNet | 0.04% | 0.9996 | 1.0 | 0.9998
Fuzzy | Random | 2.99% | 0.9701 | 0.9737 | 0.9719
CE ResNet | 1.46% | 0.9854 | 0.9901 | 0.9878
SupCon ResNet | 0.47% | 0.9953 | 0.9902 | 0.9927
Malfunction | Random | 0.28% | 0.9972 | 0.9994 | 0.9983
CE ResNet | 0.13% | 0.9987 | 0.9989 | 0.9988
SupCon ResNet | 0.04% | 0.9996 | 1.0 | 0.9998
Overall | Random | 0.93% | 0.9908 | 0.9923 | 0.9915
CE ResNet | 0.44% | 0.9956 | 0.9968 | 0.9962
SupCon ResNet | 0.17% | 0.9984 | 0.9974 | 0.9979

###  V-F Model Complexity Analysis

Along with high detection accuracy, a CAN-bus IDS must follow hardware
constraints, such as low computing power, small memory size, and small
response time, such that it can be installed and deployed in an ECU. Table V
shows the complexity of the proposed model in terms of memory and running
time, compared to the Inception ResNet proposed by [7]. From the results, our
model is more lightweight and runs faster than the state of the art of CAN bus
IDS. The model spends approximately 3.65 ms when inferencing one frame
including 29 messages on a 64-bit Intel (R) Core(TM) i7-7700 CPU @ 3.6 GHz.
This implies that the model can process up to approximately 8000 messages
within a second. Simultaneously, there are approximately 2000 CAN messages on
the bus. Therefore, our proposed model has enough capability to be deployed in
a real ECU.

TABLE V: Comparison of model complexity Model |  Total parameters (M) |  Total size (MB) |  Inference time (ms)
---|---|---|---
Inception ResNet [7] |  1.69 |  9.33 |  4.49
SupCon ResNet |  0.70 |  5.13 |  3.65

##  VI Conclusion

This study aims to develop a lightweight and efficient CAN bus intrusion
detection system capable of detecting and identifying specific types of
attacks. We introduced a SupCon ResNet, which was trained with the supervised
contrastive loss, by manipulating the CAN ID sequences in binary form. To the
best of our knowledge, this is the first study that combined supervised
contrastive learning and transfer learning to the in-vehicle IDS. The results
of experiments on real datasets indicate that the supervised contrastive loss
significantly reduced the false-negative rate, especially for the fuzzy
attack. Consequently, a lightweight deep learning model, which occupies low
memory size and runs faster, can be deployed in an ECU to serve the intrusion
detection task without detection accuracy degradation. Moreover, the proposed
model is useful for transfer learning, solving the challenge of data lacking
in the new release car model. We showed that the transfer learning from the
SupCon ResNet model is superior to other baselines. Concretely, the proposed
model decreases the false-negative rate by 4 and 2.5 times when tested with
KIA Soul and Chevrolet Spark models, respectively. From these results, we
concluded that the SupCon ResNet can deal with the multiple attack
classification on the CAN bus. Furthermore, the model saves time and labor
costs for data collection for the new car model. However, the proposed system
contains a limitation of no unknown attack detection. The model should be
frequently updated to detect new types of attacks since the multiclass
classification cannot be solved using unsupervised models. The challenge of
adding a new attack type without training the model from the beginning can be
a promising direction for future studies. In addition, the idea of
personalized federating learning can be adopted to CAN IDS such that each
vehicle can contribute its own data to build a competent global model for all
vehicles.

## References

  * [1] Karl Koscher et al.  “Experimental Security Analysis of a Modern Automobile”  In _2010 IEEE Symposium on Security and Privacy_ , 2010, pp. 447–462  URL: http://www.autosec.org/
  * [2] Tobias Hoppe, Stefan Kiltz and Jana Dittmann  “Security threats to automotive CAN networks—Practical examples and selected short-term countermeasures”  In _Reliability Engineering & System Safety_ 96.1 Elsevier, 2011, pp. 11–25  DOI: 10.1016/J.RESS.2010.06.026
  * [3] Hyo Jin Jo and Wonsuk Choi  “A Survey of Attacks on Controller Area Networks and Corresponding Countermeasures”  In _IEEE Transactions on Intelligent Transportation Systems_ Institute of ElectricalElectronics Engineers Inc., 2021  DOI: 10.1109/TITS.2021.3078740
  * [4] Gianni D’Angelo, Arcangelo Castiglione and Francesco Palmieri  “A cluster-based multidimensional approach for detecting attacks on connected vehicles”  In _IEEE Internet of Things Journal_ Institute of ElectricalElectronics Engineers Inc., 2020  DOI: 10.1109/JIOT.2020.3032935
  * [5] Abdelouahid Derhab et al.  “Histogram-Based Intrusion Detection and Filtering Framework for Secure and Safe In-Vehicle Networks”  In _IEEE Transactions on Intelligent Transportation Systems_ , 2021, pp. 1–14  DOI: 10.1109/TITS.2021.3088998
  * [6] Omid Avatefipour et al.  “An intelligent secured framework for cyberattack detection in electric vehicles’ can bus using machine learning”  In _IEEE Access_ 7 Institute of ElectricalElectronics Engineers Inc., 2019, pp. 127580–127592  DOI: 10.1109/ACCESS.2019.2937576
  * [7] Hyun Min Song, Jiyoung Woo and Huy Kang Kim  “In-vehicle network intrusion detection using deep convolutional neural network”  In _Vehicular Communications_ 21 Elsevier Inc., 2020, pp. 100198  DOI: 10.1016/j.vehcom.2019.100198
  * [8] Araya Kibrom Desta, Shuji Ohira, Ismail Arai and Kazutoshi Fujikawa  “Rec-CNN: In-vehicle networks intrusion detection using convolutional neural networks trained on recurrence plots”  In _Vehicular Communications_ 35 Elsevier BV, 2022, pp. 100470  DOI: 10.1016/J.VEHCOM.2022.100470
  * [9] Adrian Taylor, Sylvain Leblanc and Nathalie Japkowicz  “Anomaly detection in automobile control network data with long short-term memory networks”  In _Proceedings - 3rd IEEE International Conference on Data Science and Advanced Analytics, DSAA 2016_ Institute of ElectricalElectronics Engineers Inc., 2016, pp. 130–139  DOI: 10.1109/DSAA.2016.20
  * [10] Javed Ashraf et al.  “Novel Deep Learning-Enabled LSTM Autoencoder Architecture for Discovering Anomalous Events From Intelligent Transportation Systems”  In _IEEE Transactions on Intelligent Transportation Systems_ Institute of ElectricalElectronics Engineers (IEEE), 2020, pp. 1–12  DOI: 10.1109/tits.2020.3017882
  * [11] Minki Nam, Seungyoung Park and Duk Soo Kim  “Intrusion Detection Method Using Bi-Directional GPT for in-Vehicle Controller Area Networks”  In _IEEE Access_ 9 Institute of ElectricalElectronics Engineers Inc., 2021, pp. 124931–124944  DOI: 10.1109/ACCESS.2021.3110524
  * [12] Sumit Chopra, Raia Hadsell and Yann Lecun  “Learning a Similarity Metric Discriminatively, with Application to Face Verification”
  * [13] Mee Lan Han, Byung Il Kwak and Huy Kang Kim  “Anomaly intrusion detection method for vehicular networks based on survival analysis”  In _Vehicular Communications_ 14 Elsevier, 2018, pp. 52–63  DOI: 10.1016/J.VEHCOM.2018.09.004
  * [14] “BOSCH CAN Specification Version 2.0”, 1991
  * [15] Kyong-Tak Cho and Kang G Shin  “Fingerprinting Electronic Control Units for Vehicle Intrusion Detection”, 2016, pp. 911  URL: https://www.usenix.org/conference/usenixsecurity16/technical-sessions/presentation/cho
  * [16] Min Ju Kang and Je Won Kang  “A novel intrusion detection method using deep neural network for in-vehicle network security”  In _IEEE Vehicular Technology Conference_ 2016-July Institute of ElectricalElectronics Engineers Inc., 2016  DOI: 10.1109/VTCSPRING.2016.7504089
  * [17] Eunbi Seo, Hyun Min Song and Huy Kang Kim  “GIDS: GAN based Intrusion Detection System for In-Vehicle Network”  In _2018 16th Annual Conference on Privacy, Security and Trust, PST 2018_ Institute of ElectricalElectronics Engineers Inc., 2018  DOI: 10.1109/PST.2018.8514157
  * [18] Liuwang Kang and Haiying Shen  “A Transfer Learning based Abnormal CAN Bus Message Detection System”  In _Proceedings - 2021 IEEE 18th International Conference on Mobile Ad Hoc and Smart Systems, MASS 2021_ Institute of ElectricalElectronics Engineers Inc., 2021, pp. 545–553  DOI: 10.1109/MASS52906.2021.00073
  * [19] Giuseppina Andresini, Annalisa Appice and Donato Malerba  “Autoencoder-based deep metric learning for network intrusion detection”  In _Information Sciences_ 569 Elsevier, 2021, pp. 706–727  DOI: 10.1016/J.INS.2021.05.016
  * [20] Qigang Liu et al.  “A multi-task based deep learning approach for intrusion detection”  In _Knowledge-Based Systems_ 238 Elsevier, 2022, pp. 107852  DOI: 10.1016/J.KNOSYS.2021.107852
  * [21] Manuel Lopez-Martin, Antonio Sanchez-Esguevillas, Juan Ignacio Arribas and Belen Carro  “Supervised contrastive learning over prototype-label embeddings for network intrusion detection”  In _Information Fusion_ 79 Elsevier, 2022, pp. 200–228  DOI: 10.1016/J.INFFUS.2021.09.014
  * [22] Ting Chen, Simon Kornblith, Mohammad Norouzi and Geoffrey Hinton  “A Simple Framework for Contrastive Learning of Visual Representations”  In _37th International Conference on Machine Learning, ICML 2020_ PartF168147-3 International Machine Learning Society (IMLS), 2020, pp. 1575–1585  DOI: 10.48550/arxiv.2002.05709
  * [23] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun  “Deep Residual Learning for Image Recognition”  In _Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ 2016-December IEEE Computer Society, 2015, pp. 770–778  DOI: 10.48550/arxiv.1512.03385
  * [24] Prannay Khosla et al.  “Supervised Contrastive Learning”, 2020  DOI: 10.48550/arxiv.2004.11362
  * [25] Sinno Jialin Pan and Qiang Yang  “A survey on transfer learning”  In _IEEE Transactions on Knowledge and Data Engineering_ 22.10, 2010, pp. 1345–1359  DOI: 10.1109/TKDE.2009.191


# Formation of Multiple Counter-propagating Clusters in the Attractive
Hamiltonian Mean-field Model

Danilo M. Rivera<EMAIL_ADDRESS>Roberto E. Navarro
<EMAIL_ADDRESS>Departamento de Física, Facultad de Ciencias Físicas y
Matemáticas, Universidad de Concepción. Concepción, Chile.

###### Abstract

Many-body long-range interacting systems can remain approximately in a quasi-
stationary state far-from-thermodynamic equilibrium. These states are
typically characterized by a pair of counter-propagating density clusters, or
by a single non-homogeneous core-halo in the phase-space of the particles. By
using particle simulations based on the Hamiltonian mean-field model, we show
that this model supports stationary states with multiple cluster or particle
holes in phase-space density. We also propose a mechanism based on wave-wave
and wave-particle interactions that lead to the formation of these clusters,
and characterize these new quasi-stationary states in terms of the initial
parameters of the simulations.

††preprint: APS/123-QED

##  I Introduction

Systems with long-range interactions can show particular characteristics both
in their evolution towards the equilibrium state and in the equilibrium state
itself [1, 2]. In particular, macroscopic quantities of the system can remain
approximately in a steady state but far from the expected Boltzmann-Gibbs
equilibrium state for long periods of time [3]. This metastable or quasi-
stationary state (QSS) is characterized by non-Gaussian distribution functions
[4], with a lifetime that diverges with the number of particles $N$ composing
the system [5, 6].

The Hamiltonian Mean-field (HMF) model is one of the most studied examples of
long-range interacting systems [7]. It has been shown that QSSs in the HMF
model can be classified according to the initial total energy $u_{0}$ and the
initial degree of homogeneity $M_{0}$[8]. An homogeneous QSS is typically
characterized by two counter-propagating clusters in phase-space. The
formation of these two counter-propagating clusters have been explained as the
result of a wave-particle interaction, where a certain group of particles
resonate with counter-propagating density waves [9].

Below a certain threshold for $u_{0}$, a non-homogeneous QSS characterized by
a core-halo structure is observed [10, 11, 12]. This threshold corresponds to
a first-order phase transition between homogeneous and non-homogeneous QSSs
for $M_{0}<0.6$. For $M_{0}>0.6$, the order of this transition is unclear
because the time-asymptotic state of the HMF model is highly dependent of the
initial conditions [13]. On the other hand, the formation of the core-halo
structure in non-homogeneous QSSs has been linked to a Landau-like damping
mechanism. Under this mechanism, particles resonate with density waves, thus
gaining a large amount of energy and forming a diffuse halo in phase-space.
This also results in damping of the density waves, making the remaining
particles to be condensed in a dense low-energy core [12].

In this paper, long-term particle simulations based on the HMF model are
performed. It is shown that the HMF model also supports QSSs characterized by
more than two counter-propagating clusters or by core-halo distributions with
an irregular halo. Depending on the initial energy $u_{0}$ and homogeneity
degree of the distribution function, an stationary cluster can also be formed
in presence of counter-propagating waves. We classify the different QSS
according to their degree of homogeneity, number, and type of structures
present in phase-space, and organized depending on the properties of the
initial distribution function. The classification considers both known and new
cases obtained as a result of molecular dynamics simulations performed for
this study.

##  II  The Hamiltonian Mean-field model

The HMF model corresponds to a system of $N$ particles of unitary mass moving
in a circle of unit radius, interacting through a potential of infinite range
[14]. The Hamiltonian describing this system is

$\displaystyle
H=\sum^{N}_{i=1}\frac{p_{i}^{2}}{2}+\frac{1}{2N}\sum^{N}_{i,j=1}[1-\cos(\theta_{i}-\theta_{j})]\,,$
(1)

where $|\theta_{i}|<\pi$ is the position of the $i$-th particle, and $p_{i}$
is its conjugate momentum. The potential energy —the second term in Eq. (1)—,
is considered to be attractive. Thus, the Hamilton equations of motion are:

$\displaystyle\dot{\theta}$ $\displaystyle=p\,,$ $\displaystyle\dot{p}$
$\displaystyle=-M_{x}\sin\theta_{i}+M_{y}\cos\theta_{i}\,,$ (2)

where $M_{x}=\sum_{i}\cos\theta_{i}/N$ and $M_{y}=\sum_{i}\sin\theta_{i}/N$.
We also define $M=\sqrt{M_{x}^{2}+M_{y}^{2}}$, usually termed as the
magnetization, which provides a measure of the degree of homogeneity of the
system. Note that the magnetization only takes values in the domain $0\leq
M\leq 1$, where $M=0$ corresponds to particles distributed homogeneously in
space, while $M=1$ is an extreme case where all particles are concentrated in
a single point in space.

The positions and momentums of the particles are randomly initialized
following a water-bag phase-space distribution, defined as

$\displaystyle
f_{0}(\theta,p)=\frac{1}{4\theta_{0}p_{0}}\Theta(\theta_{0}-|\theta|)\Theta(p_{0}-|p|)\,,$
(3)

where $\Theta$ is the Heaviside function, and $\theta_{0}\geq 0$ and
$p_{0}\geq 0$ are the boundary values of the initial domain of particle
positions and velocities, respectively. The initial conditions can be
characterized by either $\theta_{0}$ and $p_{0}$, or by the initial energy
$u_{0}=p_{0}^{2}/6+(1-M_{0}^{2})/2$ and the initial magnetization
$M_{0}=M_{x}(0)=\sin(\theta_{0})/\theta_{0}$. Here, $M_{y}(0)=0$ since $f_{0}$
is symmetric with respect to $\theta$. In the thermodynamic limit
$N\rightarrow\infty$, it can be shown that $M_{y}(t)=0$ at all times [1].

To study the QSS of the HMF model, $816$ particle simulations performed. We
used $N=10^{6}$ particles, whose equations of motion Eq. (2) are solved using
a second-order symplectic integrator with a time-step $\Delta t=0.01$ and
$2^{17}$ integration steps (corresponding to a maximum time of
$t_{\text{max}}=1310.72$). The initial conditions were chosen in the range
$0<M_{0}<1$ and $\frac{1}{2}(1-M_{0}^{2})<u_{0}<2$, where the lower bound in
$u_{0}$ corresponds to the minimum value of the total energy allowed for the
HMF model when initialized with a water-bag particle distribution. The upper
bound was chosen arbitrarily, however for values of $u_{0}>2$ no different
conclusions were obtained.

###  II.1 Magnetization characteristics

Figure 1 shows the time evolution of the magnetization $M_{x}(t)$ for six
representative cases of this study, all initialized with the same
magnetization $M_{0}=M_{x}(0)=0.8$. Figures 1(a)–(d) correspond to simulations
initialized with energies between $0.78\leq u_{0}\leq 1.45$, showing that the
magnetization violently decreases and start oscillating around the time-
averaged $\langle M_{x}\rangle=0$. Although this is expected value of the
magnetization at equilibrium [2], we will see later that the phase-space
distribution function supports non-Gaussian characteristics for a long time,
but eventually the distribution will relax to a Gaussian [15]. Although the
figures show results up to $t<410$, this state remains until the end of the
simulation at $t_{\text{max}}=1310.72$. Due to the magnetization fluctuating
around $M_{x}=0$, these cases are usually associated to homogeneous or
paramagnetic QQSs [1]. However, the term is misleading because the
magnetization also exhibit strong fluctuations with amplitudes $|\delta
M_{x}(t)|$ of the order of 0.1. These fluctuations stabilize after $t>200$,
with oscillating frequency increasing with the initial energy $u_{0}$. Figures
1(a) and (b) show clearly that the magnetization fluctuates at a single
frequency. Meanwhile, Figs. 1(c) and (d) show that the fluctuations in
$M_{x}(t)$ display irregular periodicity, suggesting a bifurcation of the
periodic motion [16].

Figure 1:  Time evolution of the magnetization $M_{x}(t)$ for $M_{0}=0.80$,
$N=10^{6}$ particles, and (a) $u_{0}=0.78$, (b) $u_{0}=0.90$, (c)
$u_{0}=1.10$, (d) $u_{0}=1.45$, (e) $u_{0}=0.3297$, and (f) $u_{0}=0.45$.
Figures in the same row share axes. Notice that figures (e) and (f) use a
different vertical range than figures (a)–(d).

Figure 1(e) shows the magnetization for $u_{0}=0.3297$ and $M_{0}=0.8$. In
this case, $M_{x}(t)$ oscillates around the initial value, and with a
fluctuating amplitude much lower than the discussed in Figs. 1(a)–(d). In
fact, this case is consistent with a generalized virial condition defined by
Benetti  _et al._ [13]. They show that a water-bag distribution initialized
with values satisfying

$\displaystyle(1-2u_{0})\theta_{0}=\sin\theta_{0}\cos\theta_{0}\,,$ (4)

will evolve in such a way that its spatial width remains invariant,
$\theta_{\text{max}}(t)=\theta_{0}$. This means that the magnetization also
remains approximately invariant.

Figure 1(d) shows the magnetization for $u_{0}=0.45$ and $M_{0}=0.8$. Here,
$M_{x}(t)$ also violently decreases as in Figs. 1(a)–(d), but then oscillates
around $\langle M_{x}\rangle\simeq 0.675$. Both cases in Figs 1(e) and (d) are
representatives of non-homogeneous or ferromagnetic QSSs [13, 12], i.e. whose
time-averaged magnetization is $\langle M_{x}\rangle\neq 0.675$. The amplitude
of the fluctuations is clearly modulated in the case shown in Fig. 1(d),
suggesting that at least two frequencies are at work in the evolution of the
system.

###  II.2 Quasi-homogeneous states

Since the magnetization fluctuates in all cases shown in Fig. 1, even when
$\langle M_{x}\rangle=0$, it means that the particle density distribution is
not exactly homogeneous at all times and that some structures (like clusters
of particles or holes) are propagating in phase-space. This can be seen in
Fig. 2, where a snapshot of the phase-space density at $t=350.08$ for the same
parameters as in Fig. 1 is shown.

Figure 2: Phase-space density at $t=350.08$ for the same parameters as in Fig.
1. Figures in the same row share axes. The colorbar represents the number of
particles per bin. White areas correspond to empty bins.

Figure 2(a) shows that the homogeneous QQS case presented in Fig. 1(a) is
characterized by two counter-propagating clusters. The observed finite-
amplitude fluctuation in the magnetization is due to the particles trying to
organize themselves collectively in phase-space but failing to do so [10]. By
using phenomenological arguments, Antoniazzi  _et al._ [9] suggested that —for
short times—, the one-particle Hamiltonian can be understood as the
Hamiltonian of one particle interacting with two counter-propagating waves
with phase-speed $\Delta_{p}$ (although they do not show how this speed
depends on the initial values of $u_{0}$ and $p_{0}$). This suggests that,
depending on the initial conditions of the system, certain particles can be
trapped in one of the two resonances produced by the propagation of these
waves in the system, leading to the concentration of particles around the
resonance regions and to the formation of these two counter-propagating
clusters. Our simulations indicate that $\Delta_{p}$ grows mostly with the
initial energy $u_{0}$. Due to this, as the initial energy is higher, the
distance between the two counter-propagating clusters also increases. On the
other hand, our simulations also indicate that the initial magnetization
mostly affects how many particles compose each traveling cluster.

Since $M_{y}(t)=0$, then Eq. (2) is consistent with the reflection symmetry
$f(\theta,p,t)=f(-\theta,-p,t)$. This means that the counter-propagating waves
must propagate at the same phase-speed and with the same amplitude. For
sufficiently low energies $u_{0}$, then their linear superposition must form a
standing wave that traps particles in the center of phase-space. This is shown
in Fig. 2(b) for $u_{0}=0.9$, where a third cluster is centered at $\theta=0$
and $p=0$. A similar case has been reported by Yamaguchi [17], Martelloni  _et
al._ [4], without a proper phenomenological explanation. This third cluster
appears early in the simulation, around the same time as the counter-
propagating clusters. Although it deforms in time due to the traveling of the
other two clusters, it conserves the area and thus the number of particles it
has trapped. In the previous case discussed in Fig. 2(a), a standing wave is
also present, but the counter-propagating waves travel a low enough speeds
$\Delta_{p}$. This means that most particles are trapped by the traveling
waves, with only few particles being trapped at the center of the phase-space.
When $\Delta_{p}$ increases, the two resonance zones in phase-space are quite
far apart from each other so that in the center of the phase-space there
remains a non-negligible population of particles that do not interact with the
two clusters, but form a third smaller stationary cluster. Thus, this case
corresponds to a homogeneous QSS, but with different characteristics to the
ones discussed in Fig. 2(a).

Figure 2(c) shows that, for $u_{0}=1.1$, a new secondary pair of counter-
propagating clusters have formed around $p=\pm 0.5$. For higher energies
$u_{0}$, non-linear wave-wave coupling is expected between the traveling and
stationary waves discussed in Figs. 2(b). When the system enters the quasi-
stationary stage, there is a population of particles forming a streamline with
speed between 0 and $\Delta_{p}$. Then, this streamline condenses, thus
forming these secondary counter-propagating clusters. As $u_{0}$ increases,
the density ratio between the secondary clusters and the stationary cluster
also increases, meaning that the non-linear coupling becomes strong enough to
carry away particles from the low-energy central cluster. This kind of QSS has
been reported in [4] and are explained through a maximum entropy scheme of the
Lynden-Bell theory [18], without linking it to a physical mechanism.

As the energy $u_{0}$ increases further, and within the simulations performed,
we also obtained homogeneous QSSs characterized by multiple small clusters
forming spontaneously between the two counter-propagating clusters [see Fig.
2(d)]. These structures are not evident during the whole quasi-stationary
stage of the system: sometimes, the particles composing these small structures
tend to condense making the clusters more evident, while at another stages
these clusters show a dilation making their observation more complex. This may
explain the irregular periodicity of the magnetization shown in Fig. 1(d).
Moreover, these small clusters do not show to be recurrent structures, since
for repeated simulations with the parameters the small structures may not
emerge in phase-space density.

###  II.3 Far-from-homogeneity states

The non-homogeneous QSS characterized by $\langle M_{x}\rangle\neq 0$. Unlike
the quasi-homogeneous cases, there are no counter-propagating clusters in
phase-space, but a single central structure composed of a dense core of
particles surrounded by a small, sparse halo of particles, as seen in the
phase-space density of Fig. 2(e). Benetti  _et al._ [13] argued that cases
initialized with initial conditions belonging to Eq. (4), are subject to a
mean-field potential that acts on each particle. This generates only
microscopic oscillations, so that parametric resonances are suppressed. This
implies that the macroscopic oscillations in the system are suppressed, and
magnetization fluctuations are limited as shown in Fig. 1(e). So, there is no
physical mechanism that can cause certain particles to move to highly
energetic regions, or to energetically low regions as a result of oscillation
damping. The particle distribution of these QSS can be well described by the
Lyndell-Bell theory [18], which has been shown to maximize the coarse-grained
entropy characteristic of this theory [13, 1].

For states slightly above or below the virial line Eq. (4), core-halo
distributions can still be observed. They have been explained as a product of
the Landau damping mechanism present in the system [12]. This implies that a
density wave propagates in the system with which certain particles interact by
gaining energies at the expense of the collective motion and reaching high
energy states forming a halo. As a result of this energy transfer, the wave is
damped by carrying a large number of particles to low energy states, forming
the dense core.

However, as shown in Fig. 2(f), the halo exhibits holes or regions free from
particles in phase-space. These particle holes orbit clock wise around dense
core. After averaging in phase-space, this results in the modulation of the
magnetization, as shown in Fig. 1. The holes are not formed until the system
has been in the quasi-stationary stage for some time. The presence of holes in
the halo have already been reported by Benetti  _et al._ [13], where they
studied regular particle trajectories in the phase-space. They argue that,
when the initial conditions do not coincide with the curve (4), parametric
resonances arise in the system, turning the dynamics chaotic and producing
macroscopic oscillations in the system. These oscillations lead to Landau
damping, which is consistent to explain the formation of the halo, but not the
partial depletion of particles in the halo. The relationship between the
modulation of the magnetization and these halo structures merits a detailed
study.

###  II.4 Classification of the phase-space distribution

Figure 3 shows the phase diagram as a function of the initial parameters
$M_{0}$ and $u_{0}$, where we have classified the different structures
observed in the QSS phase-space distribution. Following the characteristics
discussed in Figs. 1 and 2, we identify: Homogenous QSS with two clusters
[H2C, as in Fig. 2(a)], three clusters [H3C, as in Fig. 2(b)], five clusters
[H5C, as in Fig. 2(c)], and two clusters with multiple spontaneous low-density
clusters [H2C*, as in Fig. 2(d)]. Also, non-homogeneous states (NHS)
characterized by core-halo distributions with particle holes in the halo.

Figure 3:  Initial conditions space. Classification diagram of the QSS
associated with different initial conditions of the HMF model. Homogeneous QSS
with: (yellow–H2C) two counter-propagating clusters; (yellow–H2C*) two
counter-propagating clusters and small non-recurrent clusters; (red–H3C) three
clusters; (blue–H5C) five clusters. (Green–NHS) Non-homogeneous QSS with a
core-halo with particle holes. The gray area containing the dashed line
corresponds to non-homogeneous QSS with a core-halo. The dashed line Eq. (4)
corresponds to non-homogeneous QSS where the magnetization of the system
remains constant at all times. Points (a)-(f) are the corresponding initial
states labeled in Figs. 1 and 2.

The area marked as H2C* in Fig. 3 occur for high initial energies
($u_{0}>1.2$) and high initial magnetizations ($M_{0}>0.5$), and we interpret
it as a sub-zone of the H2C zone. There is no unique classification category
for this sub-region, since it was not possible to distinguish unequivocally
all the characteristics of phase-spaces densities. However, due to the high
energies, the formation of spontaneous clusters may be linked to some physical
mechanism different from the one responsible for the two main clusters.

On the other hand, Eq. (4) determines the values (dashed lines in Fig. 3) for
which the magnetization does not change from the initial condition $M_{0}$,
and where the phase-space distribution function is perfectly bounded. Due to
numerical uncertainties, this curve was experimentally verified (to a small
error) with the gray area in Fig. 3. Note that independent of the previously
commented curve (4) is in the center of the gray zone NH which shows a
concordance between the theory and the simulations.

##  III Conclusion

Many-body long-range interacting systems can remain approximately in a quasi-
stationary state far-from-thermodynamic equilibrium. These states are
typically characterized by a pair of counter-propagating density clusters, or
by a single non-homogeneous core-halo in the phase-space of the particles.

In this work, we analyze a set of $816$ particle simulations based on the
Hamiltonian mean-field model. The initial conditions are characterized by a
water-bag distribution of $N=10^{6}$ particles, energy $u_{0}$ and
magnetization $M_{0}$. The magnetization is a measure of the homogeneity
degree of the space density of particles.

We observe that the quasi-stationary states (QSS) of the HMF model are
characterized by strong fluctuations in the magnetization around its time
averaged value $\langle M_{x}(t)\rangle$, with amplitude fluctuations of the
order of $0.1$. Thus, we will refer as quasi-homogeneous QSS to the cases when
$\langle M_{x}(t)\rangle=0$ and with finite-amplitude fluctuations. As the
initial energy $u_{0}$ increases, the frequency of these fluctuations also
increases. At some high-enough value of $u_{0}$, the fluctuations display
irregular periodicity, suggesting a bifurcation of the periodic motion [16].

For quasi-homogeneous QSS, i.e. with fluctuating magnetization around $\langle
M_{x}(t)\rangle=0$, the phase-space particle density distribution is
characterized by multiple pairs of counter-propagating clusters, whose number
tends to increase with the initial energy $u_{0}$. Since the distribution
function must follow a reflection symmetry $f(x,p)=f(-x,-p)$ at all times,
then each pair of counter-propagating clusters must propagate at the same
speed and with the same density. This suggests the presence of counter-
propagating wave densities that superpose to form a standing wave, thus
trapping particles near the center $x=0$ and $p=0$ of the phase-space
distribution. As the energy $u_{0}$ increases, the pair of counter-propagating
clusters separate from each other in momentum space. Thus, through non-linear
wave-wave interactions, secondary clusters are formed in phase-space. For
high-enough energies, the distribution is characterized by a flat homogeneous
momentum distribution for velocities bounded by the two principal counter-
propagating waves, although spontaneous bumps can form through the simulation.

For non-homogeneous QSS and for distributions far from the virial condition
Eq. (4), the magnetization exhibits modulated finite-amplitude fluctuations
around $\langle M_{x}(t)\rangle\neq 0$. This suggests that at least two
frequencies are at work through the simulation. The phase-space density in
this case is characterized by a relatively dense core of particles, and a
relatively tenuous halo probably formed by a Landau-like damping mechanism. In
some portions of the halo, holes or regions free-from-particles can be seen,
that orbit in a clock-wise direction around the core. The holes are not formed
until the system has been in the quasi-stationary stage for some time, much
after the core and halo have formed.

Finally, the above description was used to characterize the QSS in terms of
the initial energy $u_{0}$, initial magnetization $M_{0}$, and the number of
cluster and their characteristics in phase-space density.

In summary, we reinterpret the notion of homogeneous QSS; provide examples
where standing clusters or multiple counter-propagating clusters exist in the
phase-space of particles; and provide possible mechanisms that may be
responsible for the formation of these clusters. Thus, we suggest that other
phase-transitions may occur, not only from a non-homogeneous to a quasi-
homogeneous QSS, but also in the formation of multiple clusters that may be
related to the periodicity and modulation of the density waves, as suggested
by [16].

###### Acknowledgements.

This work has been financially supported by ANID, Chile through the projects
FONDECyT No. 11180947 and 1191351 (R.E.N.). We thank Nicol Gutiérrez for
useful comments.

## References

  * Levin  _et al._ [2014] Y. Levin, R. Pakter, F. B. Rizzato, T. N. Teles, and F. P. Benetti, Nonequilibrium statistical mechanics of systems with long-range interactions, Physics Reports 535, 1 (2014).
  * Campa  _et al._ [2009] A. Campa, T. Dauxois, and S. Ruffo, Statistical mechanics and dynamics of solvable models with long-range interactions, Physics Reports 480, 57 (2009).
  * Pluchino  _et al._ [2007] A. Pluchino, A. Rapisarda, and C. Tsallis, Nonergodicity and central-limit behavior for long-range hamiltonians, Europhysics Letters (EPL) 80, 26002 (2007).
  * Martelloni  _et al._ [2016] G. Martelloni, G. Martelloni, P. de Buyl, and D. Fanelli, Generalized maximum entropy approach to quasistationary states in long-range systems, Phys. Rev. E 93, 022107 (2016).
  * Roupas [2020] Z. Roupas, Statistical mechanics of gravitational systems with regular orbits: rigid body model of vector resonant relaxation, Journal of Physics A: Mathematical and Theoretical 53, 045002 (2020).
  * Schütz and Morigi [2014] S. Schütz and G. Morigi, Prethermalization of atoms due to photon-mediated long-range interactions, Phys. Rev. Lett. 113, 203002 (2014).
  * Santini  _et al._ [2022] A. Santini, G. Giachetti, and L. Casetti, Violent relaxation in the hamiltonian mean field model: II. non-equilibrium phase diagrams, Journal of Statistical Mechanics: Theory and Experiment 2022, 013210 (2022).
  * Bachelard  _et al._ [2008] R. Bachelard, C. Chandre, D. Fanelli, X. Leoncini, and S. Ruffo, Abundance of regular orbits and nonequilibrium phase transitions in the thermodynamic limit for long-range systems, Phys. Rev. Lett. 101, 260603 (2008).
  * Antoniazzi  _et al._ [2007] A. Antoniazzi, D. Fanelli, J. Barré, P.-H. Chavanis, T. Dauxois, and S. Ruffo, Maximum entropy principle explains quasistationary states in systems with long-range interactions: The example of the hamiltonian mean-field model, Phys. Rev. E 75, 011112 (2007).
  * Bachelard  _et al._ [2009] R. Bachelard, C. Chandre, A. Ciani, D. Fanelli, and Y. Yamaguchi, Analytical results on the magnetization of the hamiltonian mean-field model, Physics Letters A 373, 4239 (2009).
  * de Buyl and Gaspard [2011] P. de Buyl and P. Gaspard, Effectiveness of mixing in violent relaxation, Phys. Rev. E 84, 061139 (2011).
  * Pakter and Levin [2011] R. Pakter and Y. Levin, Core-halo distribution in the hamiltonian mean-field model, Phys. Rev. Lett. 106, 200603 (2011).
  * Benetti  _et al._ [2012] F. P. d. C. Benetti, T. N. Teles, R. Pakter, and Y. Levin, Ergodicity breaking and parametric resonances in systems with long-range interactions, Phys. Rev. Lett. 108, 140601 (2012).
  * Antoni and Ruffo [1995] M. Antoni and S. Ruffo, Clustering and relaxation in hamiltonian long-range dynamics, Phys. Rev. E 52, 2361 (1995).
  * Pluchino  _et al._ [2004] A. Pluchino, V. Latora, and A. Rapisarda, Metastable states, anomalous distributions and correlations in the HMF model, Phys. D 193, 315 (2004).
  * Morita and Kaneko [2006] H. Morita and K. Kaneko, Collective oscillation in a hamiltonian system, Phys. Rev. Lett. 96, 050602 (2006).
  * Yamaguchi [2011] Y. Y. Yamaguchi, Construction of traveling clusters in the hamiltonian mean-field model by nonequilibrium statistical mechanics and bernstein-greene-kruskal waves, Phys. Rev. E 84, 016211 (2011).
  * Lynden-Bell [1967] D. Lynden-Bell, Monthly Notices of the Royal Astronomical Society 136, 101 (1967).