CAN-BERT do it? Controller Area Network Intrusion Detection System based on BERT Language Model
Natasha Alkhatib, Maria Mushtaq, Hadi Ghauch, Jean-Luc Danger

Télécom Paris, IP Paris, Palaiseau, France

qcr{natasha.alkhatib, maria.mushtaq, hadi.ghauch<EMAIL_ADDRESS>
Due to the rising number of sophisticated customer functionalities, electronic control units (ECUs) are increasingly integrated into modern automotive systems.
However, the high connectivity between the in-vehicle and the external networks paves the way for hackers who could exploit in-vehicle network protocols' vulnerabilities.
Among these protocols, the Controller Area Network (CAN), known as the most widely used in-vehicle networking technology, lacks encryption and authentication mechanisms, making the communications delivered by distributed ECUs insecure.
Inspired by the outstanding performance of bidirectional encoder representations from transformers (BERT) for improving  many natural language processing tasks, we propose in this paper “CAN-BERT", a deep learning based network intrusion detection system, to detect cyber attacks on CAN bus protocol.
We show that the BERT model can learn the sequence of arbitration identifiers (IDs) in the CAN bus for anomaly detection using the “masked language model" unsupervised training objective.
The experimental results on the “Car Hacking: Attack & Defense Challenge 2020" dataset show that “CAN-BERT" outperforms state-of-the-art approaches. In addition to being able to identify in-vehicle intrusions in real-time within 0.8 ms to 3 ms w.r.t CAN ID sequence length, it can also detect a wide variety of cyberattacks with an F1-score of between 0.81 and 0.99.

controller area network, CAN, Intrusion Detection, bidirectional encoder representations from transformers, BERT, in-vehicle network, cyberattacks.

§ INTRODUCTION

To fulfill automotive features, the Controller Area Network (CAN) bus is widely used as the de-facto standard for message communication between different electronic control units (ECUs) in today's vehicles. It is sometimes referred to as a "message-based" system, since it focuses on the transmission of diagnostic, informative and controlling data through messages, also known as CAN data frames. In fact, while developing a vehicle, all conceivable CAN bus messages and their respective priority, encoded into an identifier called “CAN ID", must be determined beforehand.
Due to the lack of authentication, any device can connect physically or wirelessly to the CAN bus, broadcast CAN data frames, and all the participants on the CAN bus can receive it without verifying its source. Consequently, since CAN security was not a design priority, many message injection attacks have become widely implemented. These attacks can interfere with the desired function of the system, shut down some connected nodes, and make the vehicle behave abnormally, putting at risk the safety of the driver and the passengers.

A CAN bus network exploited by attackers through different physical and wireless interfaces.

To address these security flaws, researchers have proposed intrusion detection as a supplementary layer of protection to specialized security solutions. By monitoring the communication between different ECUs within a CAN bus system, a network intrusion detection system (N-IDS) can detect deviations from the normal message exchange behavior and, thereby, identify both anticipated and novel cyberattacks.

CAN packet structure.

The adoption of deep neural networks for in-vehicle intrusion detection have lately proliferated, with impressive results. Since a message injection attack can alter the normal order of occurence of several CAN IDs, researchers have deployed deep learning based sequential models, to model the CAN ID sequences.
Some studies have proposed the use of autoregressive models that are trained to capture the patterns of regular CAN ID sequences by predicting the future CAN ID sequence based on the preceding one, such as recurrent neural network (RNN) models and its variants and the generative pretrained transformer (GPT).
However, these models identify malicious network intrusions on CAN ID traffic by focusing primarily on the exchange of CAN ID messages from earlier steps rather than integrating the left and right context of a CAN ID sequence, limiting the model's capacity to grasp the whole context information representation. Additionally, these algorithms focus largely on the correlation between CAN ID messages in normal sequences, which would result in false alarms for intrusion detection whenever the correlation is breached. Hence, due to these limitations, the autoregressive models do not adequately depict the natural communication behavior between the various ECUs.

To address these challenges, we propose CAN-BERT, an intrusion detection system based on a language representation model called Bidirectional Encoder Representations from Transformers (BERT). CAN-BERT, in contrast to autoregressive models, is a self-supervised model which can successfully depict deep bidirectional representations from CAN ID sequences by conditioning on both left and right context in its various layers. By using the “masked language model" unsupervised training objective, CAN-BERT model masks some of the CAN IDs in the input at random, with the goal of predicting the conventional ID of the masked word based on its left and right context.

We evaluate our approach using the recently published “Car Hacking: Attack & Defense Challenge 2020" collected from three different cars, Chevrolet Spark, Hyundai Sonata and Kia Soul and which contains diverse types of message injection attacks.

Our contributions are summarized below:

* Inspired by the outstanding performance of BERT model for improving  many natural language processing tasks, we propose “CAN-BERT", a novel BERT-based intrusion detection system architecture which can detect known and unprecedented cyberattacks in CAN ID sequences.

* We evaluate the performance of our approach with the recently published “Car Hacking: Attack & Defense Challenge 2020" collected from three different cars, Chevrolet Spark, Hyundai Sonata and KIA Soul and which contains diverse types of message injection attacks. We also compare our model with other baseline models.

Towards this end, our paper is organized into six sections. In Section <ref>, we present an overview of the Controller Area Network (CAN) and the Bidirectional Encoder Representations from Transformers model (BERT). In Section <ref>, we present an overview of our proposed framework “CAN-BERT". Section <ref> discusses the launched experiments with the corresponding dataset and the proposed metrics for IDS' evaluation. In Section <ref>, we discuss the obtained results showing the proposed model accuracy and complexity. Finally, we conclude our paper with future work direction.

§ RELATED WORK

Intrusion Detection systems (IDSs) have been widely used to detect intrusions on the Controller Area Network (CAN). Intrusions can be detected either by inspecting the content or the signals transmitted by CAN data frames [19], or by examining the order by which different CAN data frames' identifiers are exchanged between the ECUs [20].

§ PRELIMINARIES

§.§ Controller Area Network

The Controller Area Network (CAN), created by BOSCH in 1983, is a potent networking technology essential for the development of useful automotive features. Due to its robustness represented by its ability to allow various ECUs to be connected in almost all areas of a car, it still prevails in vehicles today. As seen in Figure, it is a bus system, meaning that all Electronic Control Units (ECUs) share the same wiring.

It is a "message-based" system in which messages, also known as CAN data frames, are transmitted between various ECUs.
As depicted in Figure. <ref>, each CAN data frame is composed of the following elements: Start of frame (SOF), identifier (ID), Remote frame transmission field (RTR), control field, data, cyclic redundancy check (CRC), delimiters (DEL) and acknowledgment fields (ACK), and end of frame fields. Each CAN bus message has a priority that is represented by the arbitration identifier field (ID) that can either be composed of 11 or 29 bits, depending on the car manufacturer and which will be mainly used in our work for detecting in-vehicle intrusions.

To avoid contention between multiple ECUs willing to transmit CAN messages in the medium, CAN employs a priority-based mechanism which allows the ECU with the highest priority/lowest value identifier to transmit before others. The procedure is termed “arbitration" because the message with the highest priority wins out over competing messages with a lower priority at the time of transmission.

ECU devices connected through CAN bus, source [1].

§.§.§ Security Weaknesses

CAN does not prohibit several ECUs from sharing the same IDs. Moreover, CAN messages are broadcast and do not contain any sender's address. Consequently, any device linked to the CAN bus can use any pre-defined ID, communicate its message without authentication or encryption, and all associated ECUs can receive it. The receiver defines whether or not a message identification causes the receiving ECU to retain and process the given data. Consequently, an attacker is able to broadcast spoofed CAN messages, eavesdrop on all the CAN traffic and collect detailed information about it, resulting in Fuzzing and Malfunction attacks.

Additionally, as previously mentionned, the CAN bus leverages the arbitration method which discerns between "dominant" (0) and "recessive" (1) bits in the message identifiers. Therefore, if several ECUs begin transmitting simultaneously, the ECU whose message begins with a greater number of dominant "0" bits will take over the CAN bus. As soon as a unit detects that the message on the bus is no longer the message it is transmitting, it halts its transmission, waits for the real transmission to conclude, then waits for the interframe gap to expire and retransmits its message. This phenomenon carries the risk that a message with a lower priority will never be delivered if the network is very congested and can be exploited by attackers to launch denial of service (DoS)/ flooding attacks.

§.§ BERT

CAN-BERT model architecture

Bidirectional Encoder Representations from Transformers (BERT), proposed by Devlin et al. [2], is a state-of-the-art language representation model which is designed to pretrain bidirectional representations from unlabeled text by jointly conditioning on both left and right context in all layers. Regarding its architecture, it is a multi-layer bidirectional Transformer encoder based on the original implementation proposed by Vaswani et al. [3].

Inspired by its outstanding performance in modeling sequential data, BERT is recently employed for sequence anomaly detection [6] [7] [8] [9] [10]. To the best of our knowledge, none of the previous works have evaluated the performance of the BERT model for in-vehicle intrusion detection on CAN protocol.

In order to detect anomalies in sequences, it's crucial to incorporate context from both left and right directions of the sequence. Sequential anomalies may be misdetected by traditional unidirectional models, such as OpenAI GPT and RNNs, where every token can only attend to context to its left. To solve this significant constraint, some researchers have proposed a shallow concatenation of both left-to-right and right-to-left architecture of the autoregressive models , such as Bi-RNN and Bi-GPT [5]. However, these approaches aren't as powerful as BERT which adopts a "masked language model" (MLM) training objective, in which input sequence tokens are randomly masked and the goal is to predict the original vocabulary id of the masked word based on its context. In contrast to denoising auto-encoders, BERT predicts the masked words instead of reconstructing the whole sequence [12].

§ PROPOSED FRAMEWORK: CAN-BERT

We propose "CAN-BERT", a pattern-based anomaly detection algorithmn, which leverages a BERT-based architecture to detect message injection intrusions in the CAN bus. As seen in Figure. <ref>, our model is composed of a multi-layer bidirectional Transformer encoder and is trained using the ”masked language model" self-supervised task to model normal CAN ID sequences. The following subsections elaborately describe the suggested framework.

§.§ Model description

Note that $\mathbf{S}=[ \mathbf{id}_{1},..., \mathbf{id}_{t},...,\mathbf{id}_{T}]$ is an observed sequence of $T$ CAN identifiers, arranged in their order of transmission in the CAN bus network, where an identifier  $\mathbf{id}_{t} \in \mathbb{ID}$ is an M-dimensional vector which denotes the CAN ID transmitted at time $t$ by an ECU, $\mathbb{ID}$ indicates the set of CAN IDs extracted from CAN messages, and M is the size of the $\mathbb{ID}$ set.

Since anomaly detection is an unsupervised learning-based technique in which only normal data are used for training, a collection of $N$ CAN ID sequences, represented as $\mathbb{D}_{training}=\{ \mathbf{S}_{1},..., \mathbf{S}_{j},...,\mathbf{S}_{N}\}$, is solely used as a training dataset.

Identifier Embeddings To feed the appropriate input to the BERT model, each identifier $\mathbf{id}^{j}_{t}$ with size (M,1) in a CAN sequence $\mathbf{S}^{j}$ is firstly projected into a $d$-dimensional space using a single linear layer, i.e.:

\begin{equation}
\mathbf{e^{j}_{t}}=\mathbf{W}^{e}\mathbf{id^{j}_{t}}+\mathbf{b}^{e}, \forall i \in \{1..T\},  \forall j \in \{1..N\}
\end{equation}

where $\mathbf{e^{j}_{t}}$ represents the identifier embedding with size ($d$,1), $\mathbf{W}^{e} \in \mathbb{R}^{d \times M}$ is the input embedding weight matrix, and $\mathbf{b}^{e} \in \mathbb{R}^{d}$ denotes the bias. Both $\mathbf{W}^{e}$ and $\mathbf{b}^{e}$ are trainable parameters.

Subsequently, the identifier's position is encoded into a $d$-dimensional positional embedding $\mathbf{p^{j}_{t}}$ using a sinusoidal function. To this end, the CAN ID, fed into the CAN-BERT model, is a summation of both the positional encoding and the input embedding :

\begin{equation}
\label{input}
\mathbf{x}^{j}_{t}=\mathbf{e^{j}_{t}}+\mathbf{p^{j}_{t}}
\end{equation}

where $\mathbf{x}^{j}_{t}$ is the total embedding $j$-th identifier in the $t$-th CAN ID sequence $\mathbf{id}^{j}_{t}$, thereby the convergence of the input sequence $\mathbf{S}^{j}$ into $\mathbf{X}^{j}=[\mathbf{x}^{j}_{1},..\mathbf{x}^{j}_{t}..,\mathbf{x}^{j}_{T}]^{T}$ with $\mathbf{X}^{j}$ a matrix with size $T \times d$.

Transformer Encoder The encoded input $\mathbf{X}^{j}$ is then delivered into a stack of $L$ transformer encoder layers, each of which has two sub-layers: a multi-head self-attention mechanism and a position-wise feed-forward network [3]. A residual connection is employed around each of these two sub-layers, followed by layer normalization [4], as follows:

\begin{equation}
\label{sublayers}
\begin{split}
\mathbf{E}^{(j,l)}=g(\mathbf{X}^{(j,l)}) + f(\mathbf{X}^{(j,l)}+g(\mathbf{X}^{(j,l)})) \\
\mathbf{H}^{(j,l)}=z(\mathbf{E}^{(j,l)}) + f(\mathbf{E}^{(j,l)}+z(\mathbf{E}^{(j,l)})) \\
\mathbf{X}^{(j,l+1)}=\mathbf{H}^{(j,l)},  \forall l < L
\end{split}
\end{equation}

where $\mathbf{E}^{(j,l)}$ represents the output of the first sublayer for the $l$-th transformer encoder layer with size $T \times d$, $\mathbf{H}^{(j,l)}$ represents the output of the second sublayer for the $l$-th transformer encoder layer with size $T \times d$, $g$ is the multi-head attention function, $z$ is the position wise feed forward function, and $f$ is the layer normalization function.

Attention We use the scaled dot-product attention proposed by [3], requiring query $\mathbf{Q}^{(j,l)}$, key $\mathbf{K}^{(j,l)}$, and value $\mathbf{V}^{(j,l)}$ representations, and which are projections of the embedded sequence $\mathbf{X}^{(j,l)} \in \mathbb{R}^{T \times d}$. In fact, we leverage the dot-product similarity to compare the query representation of a given CAN identifier to all other keys. Hence, if the query and key are comparable have a high attention weight, the matching value is deemed to be relevant. The output is therefore computed as a weighted sum
of the values $\mathbf{V}$:

\begin{equation}
\begin{split}
Attn(\mathbf{Q}^{(j,l)},\mathbf{K}^{(j,l)},\mathbf{V}^{(j,l)})=\sigma(\frac{\mathbf{Q}^{(j,l)}\mathbf{K}^{(j,l)T}}{\sqrt{d}})\mathbf{V^{(j,l)}} \\
\end{split}
\end{equation}

where $\sigma$ is the softmax function, $\mathbf{A} \in \mathbb{R}^{T \times T}$ denotes the attention weight matrix containing attention weights, and $d$ is the dimension of the $\mathbf{Q}^{(j,l)}$, $\mathbf{K}^{(j,l)}$,$\mathbf{V}^{(j,l)}$ vectors.

As described by [3], the multiple heads of attention allows the model to concurrently capture diverse aspect of data at distinct CAN IDs. Hence, we adopt a multi-head attention (MHA) mechanism in which the $d$-dimensional CAN identifers are projected into subspaces calculated by different attention heads $n \in \{1,..,H\}$:

\begin{equation}
\begin{split}
\mathbf{Q}^{(j,n,l)}= \mathbf{X}^{(j,l)}\mathbf{W}^{(Q,n)}, \mathbf{Q}^{(j,n,l)} \in \mathbb{R}^{T \times F} \\
\mathbf{K}^{(j,n,l)}= \mathbf{X}^{(j,l)}\mathbf{W}^{(K,n)}, \mathbf{K}^{(j,n,l)} \in \mathbb{R}^{T \times F} \\
\mathbf{V}^{(j,n,l)}= \mathbf{X}^{(j,l)}\mathbf{W}^{(V,n)}, \mathbf{V}^{(j,n,l)} \in \mathbb{R}^{T \times F}
\end{split}
\end{equation}

where $\mathbf{Q}^{(j,n,l)}$, $\mathbf{K}^{(j,n,l)}$ and $\mathbf{V}^{(j,n,l)}$ are the query, key and value vectors, respectively of the $j-th$ CAN ID sequence for the $l$-th transformer encoder layer and which are calculated using the $n$-th attention head. The $\mathbf{W}^{(Q,n)}$, $\mathbf{W}^{(K,n)}$ and $\mathbf{W}^{(V,n)}$ are their corresponding trainable weight matrices $\in \mathbb{R}^{d \times F}$, and $F$ is set to $D/H$. The results are then concatenated and projected back into representation space using the weight matrix $\mathbf{W}^{o} \in \mathbb{R}^{HF \times D}$ as follows:

\begin{equation}
\label{attention_head}
\mathbf{head}^{(j,l)}_{n}=Attn(\mathbf{Q}^{(j,n,l)},\mathbf{K}^{(j,n,l)},\mathbf{V}^{(j,n,l)})
\end{equation}

\begin{equation}
\label{concat_head}
\overline{\mathbf{X}}^{(j,l)}=[\mathbf{head}^{(j,l)}_{1},..\mathbf{head}^{(j,l)}_{n},..,\mathbf{head}^{(j,l)}_{H}]\textbf{W}^{O}
\end{equation}

where $\overline{\mathbf{X}}^{(j,l)}  \in \mathbb{R}^{(T,d)}$.

Position-wise feed-forward A position-wise feed-forward network with a ReLU activation is thereby applied to each representation in each of the layers of our encoder, in addition to attention sub-layers, using the following equation:

\begin{equation}
\label{ffn}
z(\mathbf{E}^{(j,l)})=[\mathbf{W_{1}}\mathbf{E}^{(j,l)} ]_{+} \circ \mathbf{W_{2}}
\end{equation}

where  $\mathbf{E}^{(j,l)}$ is previously defined in (<ref>), $\textbf{W}_{1}$ and $\textbf{W}_{2}$ are trainable projection matrices, where $\circ$ is the hadarmard product, and $[ ]_{+}$ is the element-wise maximum.

After passing through different transformer layers, the $L$-th contextual embedding vectors of the CAN IDs, denoted as $\mathbf{h}^{(j,L)}_{t}$ with size $(d,1) \in \mathbf{H}^{(j,L)}=[\mathbf{h}^{(j,L)}_{1},..,\mathbf{h}^{(j,L)}_{T}]^{T}$, are fed into a single linear layer which projects them back to the $M$-dimensional layer, as follows:

\begin{equation}
\mathbf{m^{j}_{t}}=\mathbf{W}^{m}\mathbf{h}^{(j,L)}_{t}+\mathbf{b}^{m}, \forall i \in \{1..T\},  \forall j \in \{1..N\}
\end{equation}

where $\mathbf{m^{j}_{t}}$ represents the projected output with size ($M$,1), $\mathbf{W}^{m} \in \mathbb{R}^{M \times d}$ is the input embedding weight matrix, and $\mathbf{b}^{e} \in \mathbb{R}^{M}$ denotes the bias. Both $\mathbf{W}^{m}$ and $\mathbf{b}^{m}$ are trainable parameters.

§.§ Training and Inference

We use the masked language model training method to train the CAN-BERT model on capturing the patterns of normal CAN ID sequences. Hence, CAN sequences with random mask as inputs, where we randomly replace a ratio of CAN IDs in a sequence with a specific MASK token, are fed into CAN-BERT. The purpose of training is to reliably anticipate the CAN IDs that have been randomly masked.

To achieve that, we feed the contextual embedding vector of the $u$-th MASK in the j-th CAN ID sequence $\textbf{m}^{j}_{MASK_{u}}$ to a softmax function, which will return a probability distribution for the whole set of CAN IDs $\mathbb{ID}$:

\begin{equation}
\label{output}
\hat{\mathbf{y}}^{j}_{[MASK_{u}]}=\sigma(\mathbf{W}^{c} \mathbf{m}^{j}_{[MASK_{u}]}+\mathbf{b}^{c})
\end{equation}

where $\hat{\mathbf{y}}^{j}_{[MASK_{u}]}$ is an $m$-dimensional vector, $\sigma$ is the softmax function, $\mathbf{m}^{j}_{[MASK_{u}]}$ and $\mathbf{b}^{c}$ are trainable parameters.

CAN-BERT is trained to minimize the cross entropy loss over a batch of $I$ sequences ( with $ I \leq N$), for masked CAN ID prediction, which is defined as:

\begin{equation}
\label{objective_function}
\mathcal{L}_{MASK}=-\frac{1}{IR}\sum^{N}_{j=1}\sum^{R}_{u=1}{\mathbf{y}}^{j}_{[MASK_{u}]}log\hat{\mathbf{y}}^{j}_{[MASK_{u}]}
\end{equation}

where the ground-truth $u$-th CAN ID in the $j$-th sequence is denoted as ${\mathbf{y}}^{j}_{[MASK_{u}]}$, $R$ is the total number of masked tokens in the $j$-th sequence, and $N$ is the number of training samples.

By modeling the normal exchange of messages through CAN bus using CAN-BERT, our model is expected, after training, to predict a candidate set with the normal CAN IDs having the highest likelihood for each masked token. Hence, for a randomly masked testing sequence, we calculate the probability distribution represented in (<ref>), for each masked CAN ID. Therefore, if the actual CAN ID is among the anticipated candidates, the corresponding CAN ID sequence is considered as normal. Otherwise, it is deemed abnormal.

§ EXPERIMENTAL SETTINGS

§.§ Dataset

To assess the proposed CAN-BERT, we leverage the “In-Vehicle Network Intrusion Detection Challenge" dataset [13] (presented in Table <ref>), which was used at the “In-vehicle Network Intrusion Detection track’ of ‘Information Security R&D Data Challenge 2019. It includes normal and abnormal in-vehicle network traffic data of HYUNDAI Sonata, KIA Soul, and CHEVROLET Spark vehicles collected during their stationary state. We have mainly used its preliminary dataset, which includes three types of attacks (Flooding, Fuzzy, and Malfunction). The dataset is labeled and each sample is represented by the following features: “Timestamp" representing the logging time, “CAN ID" representing the CAN Identifier, “DLC" indicating the Data length code, and the Payload indicating the “CAN data" field.

§.§.§ Attacks

The dataset contains the following attacks:

* Flooding Attack The flooding attack was carried out by injecting many messages with the CAN ID set to 0x000 into the CAN network. Consequently, an ECU that transmits CAN data frames with such CAN ID dominates the CAN bus, which could restrict the communications among the ECU nodes and impair normal in-vehicle functions.

* Fuzzy Attack To implement fuzzy attacks, the attacker injected every 0.0003 seconds random CAN packets, for both the ID field and the Data field. This process lead to abnormal automotive functionalities behavior such as short beeping sound repeatedly occurring, the heater turning on, etc.

* Malfunction Attack The malfunction attack was carried out by injecting messages with a specified CAN ID from among the extractable CAN IDs of a particular vehicle in order to disable relevant automotive functions, such as IDs 0×316, 0×153 and 0×18E for the HYUNDAI YF Sonata, KIA Soul, and CHEVROLET Spark vehicles, respectively.

[group style=group size=4 by 1,

[ymin=70,ymax=100,ylabel=F1-Score, xtick=1,...,4,ytick pos=left,title=$m= 0.15$,xticklabels=1,2,4,8,xlabel=$h$]
[ymin=70,ymax=100,title=$m= 0.3$,xtick=1,...,4,xticklabels=1,2,4,8,xlabel=$h$]

[ymin=70,ymax=100,title=$m= 0.45$,xtick=1,...,4,xticklabels=1,2,4,8,xlabel=$h$]

[ymin=70,ymax=100,title=$m= 0.6$,xtick=1,...,4, xlabel=$h$, xticklabels=1,2,4,8,legend style=at=($(-0.5,0)-((3cm,1.4cm)$),legend

Hyperparameter tuning the mask ratio $m$ and the number of attention heads $h$ on the “CHEVROLET Spark" dataset for $T$=32. We have obtained similar behavior pattern for the results w.r.t other sequence length and car types.

In-Vehicle Network Intrusion Detection Dataset

2*Vehicle     Dataset     # Normal     # Abnormal     Size

packets     packets     (MB)

Attack Free     136,933     N/A     6.2

CHEVROLET                         Flooding     70,001     14,999     4.2

Spark                       Fuzzy     37,957     3,043     2.0

Malfunction     47,005     3,995     2.5

Attack Free     117,172     N/A     5.8

HYUNDAI                      Flooding     78,907     17,093     4.9

Sonata                  Fuzzy     78,905     9,095     4.5

Malfunction     78,798     8,202     4.5

Attack Free     192,515     N/A     9.3

KIA                Flooding     103,928     16,072     6.2

Soul               Fuzzy     122,387     21,613     7.4

Malfunction     108,230     4,770     5.8

As mentioned in Section <ref>, we aim to detect if a sequence of ordered CAN ID contains injected messages. Hence, in order to represent CAN ID sequences, we use the Feature-based Sliding Window (FSW) to group CAN IDs which belong to the dataset into subsequences with fixed window size $\mathbf{T}$, where $\mathbf{T}$ $\in \{16,32,64,128,256\}$ and the slide size is 1. Moreover, each CAN ID sequence $\mathbf{S}=[ \mathbf{id}_{1},..., \mathbf{id}_{t},...,\mathbf{id}_{T}]$ has its corresponding labels represented by $\mathbf{Y}=[ \mathbf{y}_{1},..., \mathbf{y}_{t},...,\mathbf{y}_{T}]$ wherein each identifier $\mathbf{id}_{t} \in \mathbf{S} $ is labeled as $\mathbf{y}_{t}=1$ if $\mathbf{id}_{t} $ is an injected identifier in $\mathbf{S}$ or as $\mathbf{y}_{t}=0$ otherwise. However, to identify the state of each sequence, we have used the following criteria:

\begin{equation*}
z = \begin{cases} \mbox{1 (abnormal)} & \textup{if }  \exists \mathbf{y}_{t}=1 , \forall t \in \{1, .., T\} \nonumber \\
\mbox{0 (normal)} & \textup{otherwise}
\end{cases}
\end{equation*}

where $z$ is the CAN ID sequence's label.

§.§ Benchmark Models

The benchmark models for evaluating the performance of different CAN ID sequence anomaly detection algorithms with CAN-BERT on the chosen dataset, are detailed in this section.

* Principal Component Analysis (PCA): PCA [14] is a feature selection model which can be used to reduce data features from $m$ dimensions to $n$. Inverting the PCA transform does not retrieve the data lost during the application of the transform. The essence of PCA-based anomaly identification is that an anomalous sample should have more loss or reconstruction error than a normal sample. In other words, the loss sustained when an anomalous sample is processed by a PCA algorithm and projected back to its dimension using PCA also should be greater than when the same procedure is performed on a normal sample.

* Isolation Forest (iForest): Isolation forest (IF), proposed by Liu at al. [15], detects anomalies using isolation rather than modelling the normal points. In fact, this technique presents a novel approach for isolating anomalies using binary trees, providing a new prospect for a speedier anomaly detector that directly targets abnormalities rather than profiling all regular instances.

* Autoencoder (AE):
The autoencoder, introduced by Rumelhart et al. [16], is a deep learning based algorithm which seeks to learn a low-dimensional feature representation space suitable for reconstructing the provided data instances. During the encoding process, the encoder maps the original data onto low-dimensional feature space, while the decoder tries to retrieve the original data from the projected low-dimensional space. Reconstruction loss functions are used to learn the parameters of the encoder and decoder networks. Its reconstruction error value must be minimized during the training of normal instances and therefore used during testing as an anomaly score. In other words, compared to the typical data reconstruction error, anomalies that differ from the majority of the data have a large data reconstruction error. In our experiments, we have tested Long short-term based memory (LSTM) and Bidirectional LSTM (BiLSTM) models with different network hyperparameters: BiLSTM-AE-4 (with 4 layers), LSTM-AE-4 (with 4 layers), and LSTM-AE-8 (with 8 layers).

§.§ Evaluation metrics

For measuring the performance of different anomaly-based IDS, we use the F1-score metric, a weighted average result of both metrics precision and recall and which is specifically used when the dataset is imbalanced. The model has a large predictive power if the F1-score is near 1.0.

Precision is the ratio of correctly classified predicted abnormal observations of all the observations in the predicted class.

\begin{equation}
Precision = \frac{TP}{TP + FP}
\end{equation}

Recall is the ratio of correctly predicted abnormal observations of all observations in the actual class.

\begin{equation}
Recall = \frac{TP}{TP + FN}
\end{equation}

Hence, the F1-score is calculated using the following equation:

\begin{equation}
F1-score= 2 \cdot \frac{Precision \cdot Recall}{Precision + Recall}
\end{equation}

Where: TP= True Positive; FP=False Positive; TN= True Negative; FN=False Negative.

[group style=group size=3 by 3,

[xtick=1,...,5,title=CHEVROLET Spark,ylabel=Flooding, ytick pos=left,xticklabels=16,32,64,128,256]

[xtick=1,...,5,title=HYUNDAI Sonata,xticklabels=16,32,64,128,256]
[xtick=1,...,5,title=KIA Soul,xticklabels=16,32,64,128,256]

[xtick=1,...,5,ylabel=Fuzzing, ytick pos=left,xticklabels=16,32,64,128,256]

[xtick=1,...,5,ylabel=Malfunction, ytick pos=left,xticklabels=16,32,64,128,256]

[xtick=1,...,5,legend style=at=($(0,0)-((3cm,1.5cm)$),legend
mark=x,]  coordinates

Comparision of the CAN-BERT model with other anomaly detection baselines using the F1-score percentage metric, for different message injection attacks applied on different car models w.r.t sequence length $T$.

§ RESULTS

To evaluate our model, we leverage the Python deep learning framework Pytorch [17]. We train and evaluate them on NVIDIA® Tesla® V100S with 32 GB HBM2 memory.

§.§ Model Configuration & Hyperparameter tuning

As presented in Table <ref>, for CAN-BERT, we have chosen the total number of transformer encoder layers as 4. In each transformer layer, the position-wise feed forward network is composed of two dense layers where the first one projects $d$=256 dimensional CAN identifier embedding into $d_{ff}$=512 dimensional space, followed by a ReLU activation. The second dense layer maps back the 512-dimensional vector into the $d$-dimensional space.

CAN-BERT model configuration

Parameter    Value

N     4

$d_{model}$     256

$d_{ff}$     512

h     1

$P_{drop}$     0.1

$m$     0.45

$\#$ Candidates     5

Optimizer     Adam

Adam $\beta_{1}$     0.9

Adam $\beta_{2}$     0.999

Learning rate     0.001

Batch size     32

$\#$ Epochs     200

Patience     10

For training, we use a batch size of 32, a learning rate of 0.001 and the Adam optimizer [18] with its default parameters $\beta_{1}=0.9$ and $\beta_{2}=0.999$. To avoid overfitting, we employ the same dropout of $P_ {drop}$=0.1 for all dropout layers in our network. Moreover, we apply early stopping for a total number of 200 epochs and a patience of 10 epochs as a form of regularization used to avoid overfitting.

The hyperparameters, including the ratio of masked CAN IDs in a sequence $m$, and $h$ the number of attention heads are tuned based on a validation set for the three car types and the different message injection attacks. As seen in Figure <ref>, raising the ratios of masked CAN IDs in the sequences from 0.1 to 0.45 somewhat improves F1 scores, however increasing the ratios further degrades performance, as is the case for $m$=0.65. Furthermore, the model performance is relatively stable by setting different attention head $h \in \{1,2,4,8\}$ values for each mask ratio $m \in \{0.15,0.3,0.45,0.6\}$. Therefore, in our in-vehicle intrusion detection use case, a single attention head is sufficiant to detect different types of intrusions.
Note that, in our experiments, we use the same ratio of masked CAN IDS $m$=0.45 and $h$=1 for both training and inference phases.

§.§ Model Accuracy

As seen in Figure <ref>, we compare performance of the CAN-BERT model with other baselines approaches for different sequence length $T$ using the F1-score metric. In fact, we varied the sequence length among values of 16, 32, 64, 128 and 256 CAN IDs per sequence in the experiments. If our approaches could identify a message injection attack in a shorter sequence length, it would be more advantageous in a practical situation. The traditional machine learning algorithms such as Isolation Forest (iForest), and Principal Component Analysis (PCA) perform poorly and maintain the same F1-score metric w.r.t sequence length. Because these models presume small datasets with a limited number of features, they fail to discover abnormalities in  high-dimensional datasets. Because of this, a significant fraction of irrelevant features may effectively disguise the underlying abnormalities in the input data, resulting in poor anomaly identification performance when dealing with large input dimensions. Meanwhile, both deep learning based models autoencoder (AE) and CAN-BERT outperformed the traditional anomaly detection models over different window sizes. However, when the length of the CAN ID sequence is increased, both models performed  oppositely. In contrast to the baseline models, our suggested model, CAN-BERT, significantly outperforms them by huge margins and obtains respectable F1 scores $\in [0.85,0.99]$, demonstrating the usefulness of using BERT-based models to capture the patterns of CAN ID sequences when $T  \geq 32$. However, on short sequence length as is the case for $T=16$, the model performs modestly with F1-score $\in [0.6,0.9]$ for different kind of attacks. These experiments, therefore, reveal that by using self-supervised training tasks, CAN-BERT can effectively model medium to long normal CAN ID sequences and accurately detect anomalous sequences.

§.§ Model Complexity

From a practical point of view, we must assess not only the classification performance but also the model complexity to check if the model's ability for real-time in-vehicle intrusion detection in CAN networks. Therefore, we assessed the inference time per sample as well as the number of parameters for the CAN-BERT model w.r.t different car types.

CAN-BERT model complexity

Vehicle    Features     Values

3*CHEVROLET Spark   Number of Parameters     2,937,422

Model Size (MB)     $[20,70]$

Inference Time (ms)     $[0.8,3.1]$

3*HYUNDAI Sonata   Number of Parameters     3,149,291

Model Size (MB)     $[20,70]$

Inference Time (ms)    $[0.8,3.5]$

3*KIA Soul    Number of Parameters     3,163,142

Model Size (MB)     $[20,70]$

Inference Time (ms)    $[0.8,3.8]$

As depicted in Table <ref>, the intrusion detection inference time varies between 0.8 and 3.1 ms w.r.t CAN ID sequence length. Hence, when considering a sequence length of 32 CAN IDs, our model detects an intrusion in 0.9 to 1 ms, which is suitable for real-time detection. Furthermore, having a size between 20MB and 70 MB and a number of parameters ranging between 2 to 3 millions, our model can be either deployed in performant ECU or even on a cloud server wirelessly connected to the vehicle.

§ CONCLUSION

Identification of intrusions within the vehicle is critical for defending it against malicious cyberattacks. In this paper, we suggest CAN-BERT, a self-supervised intrusion detection system based on BERT model, for in-vehicle intrusion detection. Experimental results on benchmark datasets for different CAN ID sequence length have shown that CAN-BERT surpasses state-of-the-art techniques for CAN ID sequence anomaly detection with an F1-score ranging between 0.81 and 0.99 for different type of attacks and is appropriate for real-time detection with an inference time ranging between 0.8 and 3 ms w.r.t CAN ID sequence length. For future work, we aim to deploy our model on embedded electronic control units (ECU) and test the model efficiency in a real vehicle environment.

[1]Matheus, Kirsten, and Thomas Königseder. Automotive ethernet. Cambridge University Press, 2021.
[2]Devlin, Jacob, et al. "Bert: Pre-training of deep bidirectional transformers for language understanding." arXiv preprint arXiv:1810.04805 (2018).
[3]Vaswani, Ashish, et al. "Attention is all you need." Advances in neural information processing systems 30 (2017).
[4]Ba, Jimmy Lei, Jamie Ryan Kiros, and Geoffrey E. Hinton. "Layer normalization." arXiv preprint arXiv:1607.06450 (2016).
[5]Nam, Minki, Seungyoung Park, and Duk Soo Kim. "Intrusion detection method using bi-directional GPT for in-vehicle controller area networks." IEEE Access 9 (2021): 124931-124944.
[6]Li, Bai, et al. "How is BERT surprised? Layerwise detection of linguistic anomalies." arXiv preprint arXiv:2105.07452 (2021).
[7]Guo, Haixuan, Shuhan Yuan, and Xintao Wu. "Logbert: Log anomaly detection via bert." 2021 International Joint Conference on Neural Networks (IJCNN). IEEE, 2021.
[8]Lee, Yukyung, Jina Kim, and Pilsung Kang. "LAnoBERT: System Log Anomaly Detection based on BERT Masked Language Model." arXiv preprint arXiv:2111.09564 (2021).
[9]Le, Van-Hoang, and Hongyu Zhang. "Log-based anomaly detection without log parsing." 2021 36th IEEE/ACM International Conference on Automated Software Engineering (ASE). IEEE, 2021.
[10]Yu, Keping, et al. "Securing critical infrastructures: Deep-Learning-Based threat detection in IIoT." IEEE Communications Magazine 59.10 (2021): 76-82.
[11]Dang, Weixia, et al. "TS-Bert: Time Series Anomaly Detection via Pre-training Model Bert." International Conference on Computational Science. Springer, Cham, 2021.
[12]Vincent, Pascal, et al. "Extracting and composing robust features with denoising autoencoders." Proceedings of the 25th international conference on Machine learning. 2008.
[13]Hyunjae Kang, Byung Il Kwak, Young Hun Lee, Haneol Lee, Hwejae Lee, Huy Kang Kim, February 3, 2021, "Car Hacking: Attack & Defense Challenge 2020 Dataset", IEEE Dataport, doi: <https://dx.doi.org/10.21227/qvr7-n418>
[14]Abdi, Hervé, and Lynne J. Williams. "Principal component analysis." Wiley interdisciplinary reviews: computational statistics 2.4 (2010): 433-459.
[15]Liu, Fei Tony, Ting, Kai Ming and Zhou, Zhi-Hua. “Isolation-based anomaly detection.” ACM Transactions on Knowledge Discovery from Data (TKDD) 6.1 (2012): 3.
[16]Rumelhart, David E., Geoffrey E. Hinton, and Ronald J. Williams. "Learning representations by back-propagating errors." nature 323.6088 (1986): 533-536.
[17]Pytorch framework. <https://pytorch.org/.>
[18]Kingma, Diederik P., and Jimmy Ba. "Adam: A method for stochastic optimization." arXiv preprint arXiv:1412.6980 (2014).
[19]Hanselmann, Markus, et al. "CANet: An unsupervised intrusion detection system for high dimensional CAN bus data." Ieee Access 8 (2020): 58194-58205.
[20]Song, Hyun Min, Jiyoung Woo, and Huy Kang Kim. "In-vehicle network intrusion detection using deep convolutional neural network." Vehicular Communications 21 (2020): 100198.


###### Proof.

Consider a random function ${\bm{f}}$ where, for each $x\in\\{\pm 1\\}^{n}$,
${\bm{f}}(x)\sim\pi_{0.25}$. We’ll show that ${\bm{f}}$ meets the desired
criteria with a strictly positive probability, proving the existence of at
least one such $f$.

Let $\mu({\bm{f}})\coloneqq\mathop{{\mathds{E}}\/}_{\bm{x}\sim\\{\pm
1\\}^{n}}[{\bm{f}}(\bm{x})]$. Then $\mu({\bm{f}})$ is the average of $2^{n}$
independent samples of $\pi_{0.25}$. Applying Hoeffding’s inequality,

$\operatorname{{Pr}}[\mu({\bm{f}})>0.5]\leq\exp(-2\cdot(0.25)^{2}\cdot
2^{n})=\exp(-2^{n}/2).$

Similarly, for any $|S|\leq n/2$ and $x\in\\{\pm 1\\}^{n}$, let
$\mu({\bm{f}},S,x)\coloneqq\mathop{{\mathds{E}}\/}_{\bm{x}\sim\\{\pm
1\\}^{n}}[{\bm{f}}(\bm{x})\mid\bm{x}=x]$. $\mu({\bm{f}},S,x)$ the average of
at least $2^{n/2}$ independent samples of $\pi_{0.25}$. Once again, by
Hoeffding’s inequality,

$\operatorname{{Pr}}[\mu({\bm{f}},S,x)\leq 0]\leq\exp(-2\cdot(0.25)^{2}\cdot
2^{n/2})=\exp(-2^{n/2}/2).$

Union bounding over all $2^{n}$ choices of $S$ and $2^{n}$ choices for $x$, we
have that ${\bm{f}}$ meets the desired criteria with probability at least

$1-\exp(-2^{n}/2)-2^{2n}\exp(-2^{n/2}/2).$

When $n\geq 10$, the above probability is strictly positive, so such an $f$
must exist. ∎

###### Proof of Lemma A.7.

Let $f$ be a function with the properties of Proposition A.8, and
$g=\textsc{And}_{k}$ return $+1$ if and only if all $k$ of its inputs are
$+1$. By Proposition 6.10, for any $r\leq n/2$, $\tilde{f}_{r}$ is the
constant $+1$ function. Therefore, for any $r_{1}+\cdots+r_{k}=R$,
$g(\tilde{f}_{r_{1}},\ldots,\tilde{f}_{r_{k}})$ is the constant $+1$ function.
However,

$\operatorname{{Pr}}_{\bm{x}\sim\mathcal{D}^{k}}[(g\circ
f)(\bm{x})=+1]=(3/4)^{k}.\qed$

###  A.3 Counterexample to Conjecture 3

######  Lemma A.9.

There is $g:\\{\pm 1\\}^{k}\to\\{\pm 1\\}$, $f:\\{\pm 1\\}^{n}\to\\{\pm 1\\}$,
distribution $\mathcal{D}$ over $\\{\pm 1\\}^{n}$, and budget $R$ for which no
$R$-junta of composed form achieves optimal error among all $R$-Juntas for
$g\circ f$ with respect to $\mathcal{D}^{k}$.

###### Proof.

We’ll set $k=2$, $g=\textsc{And}_{2}$. Let $p:\\{\pm 1\\}^{2}\to[0,1]$ be
defined as

$p(x)\coloneqq\begin{cases}1&\text{if }x_{1}=x_{2}=1,\\\ 3/4&\text{if
}x_{1}\neq x_{2},\\\ 3/5&\text{if }x_{1}=x_{2}=-1.\end{cases}$

We begin by describing a probabilistic construction: Given the input $x$, the
value of ${\bm{f}}(x)$ will still be a random variable. In particular, we set
$n=2$, and ${\bm{f}}(x)$ is set to $+1$ with probability $p(x)$ and $-1$
otherwise. This probabilistic construction will later be derandomized. We
allow a junta budget of $R=4$.

Next, we construct an optimal approximator for $g\circ{\bm{f}}$. Given an
input $x^{(1)},x^{(2)}$, let $\bm{y}_{1}={\bm{f}}(x^{(1)})$ and
$\bm{y}_{2}={\bm{f}}(x^{(2)})$. For succinctness, we’ll use $p_{i}$ to refer
to the $\operatorname{{Pr}}[\bm{y}_{i}=1]$. Then, since $g=\textsc{And}_{2}$,
the optimal approximator will return $1$ iff $p_{1}p_{2}\geq 1/2$. For our
particular ${\bm{f}}$ the only choices for $p_{i}$ are $3/5,3/4,1$. As a
result,

$h^{(\mathrm{opt})}(p_{1},p_{2})=\begin{cases}1&\text{if }p_{1}=1\text{ or
}p_{2}=1,\\\ 1&\text{if }p_{1}=p_{2}=3/4,\\\ 0&\text{otherwise}.\end{cases}$

However, no composed form can achieve the above optimal approximator. Recall
that composed form approximators are of the form $h(q_{1},q_{2})$, where each
$q_{i}$ has range $\\{\pm 1\\}$. The fact that the size of this range is $2$,
but there are three possible choices $(3/5,3/4,1)$ for $p_{i}$, is the crux of
the issue.

In more detail, of the three choices $(3/5,3/4,1)$ for $p_{i}$, $q_{1}$ must
classify at least two of them the same way. This gives three cases.

  1. 1.

If $q_{1}$ classifies $3/4$ and $1$ the same way, $h(q_{1},q_{2})$ cannot
distinguish between $p_{1}=3/4,p_{2}=3/5$ and $p_{1}=1,p_{2}=3/5$, and so
cannot be optimal.

  2. 2.

If $q_{1}$ classifies $3/5$ and $3/4$ the same way, $h(q_{1},q_{2})$ cannot
distinguish between $p_{1}=3/4,p_{2}=3/4$ and $p_{1}=3/5,p_{2}=3/4$, and so
cannot be optimal.

  3. 3.

If $q_{1}$ classifies $3/5$ and $1$ the same way, $h(q_{1},q_{2})$ cannot
distinguish between $p_{1}=3/5,p_{2}=3/4$ and $p_{1}=1,p_{2}=3/4$, and so
cannot be optimal.

In all three cases composed form cannot achieve optimal error. It will always
be off by some constant.

To derandomize this construction, we set $n\gg 2$ sufficiently large. For each
$x\in\\{\pm 1\\}^{n}$, we sample the value $f(x)$ to be $+1$ with probability
$p(x_{1},x_{2})$ and $-1$ otherwise. Note that after randomly selecting the
value of $f$ on each input $x\in\\{\pm 1\\}^{n}$, $f$ is now a deterministic
function. Following the same arguments as in Proposition A.8, with high
probability over the random choices in defining $f$, the error of the optimal
$4$-junta and of the optimal composed form $4$-junta for $g\circ f$ are within
$\pm\varepsilon(n)$ of what they are for $g\circ{\bm{f}}$, where
$\varepsilon(n)$ goes to $0$ as $n\to\infty$. Therefore, for sufficiently
large $n$, there exists an $f$ meeting the desired criteria. ∎


# First principles determination of some static and dynamic properties of the
liquid 3$d$ transition metals near melting

B.G. del Rio, C. Pascual, O. Rodriguez, L.E. González, D.J. González

(Received January 31, 2020, in final form March 26, 2020)

###### Abstract

We report an ab initio molecular dynamics simulation study of several static
and dynamic properties of the liquid 3d transition metals. The calculated
static structure factors show qualitative agreement with the available
experimental data, and its second peak displays an asymmetric shape which
suggests a significant local icosahedral short-range order. The dynamical
structure reveals propagating density fluctuations whose dispersion relation
has been evaluated; moreover, its long wavelength limit is compatible with
their respective experimental sound velocity. Results are reported for the
longitudinal and transverse current spectral functions as well as for the
respective dispersion relations. We also analyze the possible appearance of
transverse-like low-energy excitations in the calculated dynamic structure
factors. Several transport coefficients have been evaluated and compared with
the available experimental data.

Key words: liquid metals, first principles calculations

###### Abstract

Ìè ïîâäîìëÿìî an initio äîñëäæåííÿ ñèìóëÿöÿìè ìîëåêóëÿðíî¿ äèíàìêè äåÿêèõ
ñòàòèíèõ äèíàìíèõ âëàñòèâîñòåé 3d ðäêèõ ïåðåõäíèõ ìåòàëâ. Îáèñëåí ñòàòèí
ñòðóêòóðí ôàêòîðè ïîêàçóþòü ÿêñíå óçãîäæåííÿ ç íàÿâíèìè åêñïåðèìåíòàëüíèìè
äàíèìè, à ¿õ äðóãèé ïê ìà àñèìåòðèíó ôîðìó, ÿêà ïåðåäáàà, ùî ñíó çíàíèé
ëîêàëüíèé êîñàåäðèíèé ïîðÿäîê. Äèíàìíà ñòðóêòóðà âèÿâëÿ ðîçïîâñþäæåííÿ
ôëóêòóàöé ãóñòèíè, äèñïåðñéíå ñïââäíîøåííÿ ÿêèõ áóëî îáèñëåíî; êðì òîãî ¿õ
äîâãîõâèëüîâà ãðàíèöÿ ñóìñíîþ ç âäïîâäíîþ ¿ì åêñïåðèìåíòàëüíîþ øâèäêñòþ çâóêó.
Ïðåäñòàâëåíî ðåçóëüòàòè äëÿ ñïåêòðàëüíèõ ôóíêöé äëÿ âäïîâäíèõ äèñïåðñéíèõ
ñïââäíîøåíü. Ìè òàêîæ àíàëçóìî ìîæëèâó ïîÿâó ïîïåðåíî-ïîäáíèõ
íèçüêîåíåðãåòèíèõ çáóäæåíü â îáèñëåíèõ äèíàìíèõ ñòðóêòóðíèõ ôàêòîðàõ. Îáèñëåíî
äåêëüêà êîåôöíòâ ïåðåíîñó ïîðâíÿíî ç íàÿâíèìè åêñïåðèìåíòàëüíèìè äàíèìè.

Ключов слова: ðäê ìåòàëè, ïåðøîïðèíöèïí ðîçðàõóíêè

##  1 Introduction

This paper presents recent ab initio molecular dynamics simulation studies for
a range of bulk static, dynamic and electronic properties of the liquid 3$d$
transition metals at thermodynamic conditions near melting. We review some of
our previous calculations [1, 2, 3, 4, 5] along with several new results.

Despite the practical relevance of the 3$d$ transition metals, it is
remarkable that the experimental information concerning their static and
dynamic properties in the liquid state is rather scarce.

The static structure factor, $S(q)$, encompasses the information about the
average short-range order of a liquid system, and was first obtained
experimentally by Waseda and coworkers [6] in the decade of the 1970s for the
3$d$ transition metals, by means of X-ray diffraction (XD) techniques. More
recently, the $S(q)$ of l-Ti, l-Fe, l-Ni and l-Zn have been determined again
by either XD or neutron diffraction (ND) experiments [7, 8, 9, 10, 11]. In
l-Zn, the new measurements showed a good agreement with the older ones, but
for l-Ni and l-Fe, the comparison with Waseda and coworkers’ XD (WXD) data for
the $S(q)$, revealed some discrepancies concerning the height of the main peak
as well as the shape of the second peak; nevertheless, the positions of the
maxima and minima of $S(q)$ were coincident in all the experiments. On the
other hand, the recent XD and ND data for the $S(q)$ of l-Ti [9, 10] yielded
results well compatible between each other, but when compared with WXD data,
they were displaced towards greater $q$-values by a significant amount of
$\approx$ 0.20 Å-1. Moreover, the recent XD data showed a second peak with a
marked shoulder on the high-$q$ side; actually, shoulders on the high-$q$ side
of the second peak of $S(q)$ have also been found experimentally in other
transition metals (Fe, Ni, Zr) [7, 8, 12]. We think that these disparities
call for new precise structural data for the other $3d$ transition metals,
both extracted from new additional experiments and/or from accurate
simulations as the ab initio ones presented here.

Several transport coefficients of most of the $3d$ liquid metals were reported
in the literature [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. The
dynamic structure factors of l-Ti, l-Fe, l-Ni and l-Zn near melting were
measured by both inelastic neutron (INS) and inelastic X-ray (IXS) scattering
[17, 26, 27, 28, 29, 30], whereas the self dynamic structure factors of l-Ni,
l-Ti and l-Cu were determined by quasi-elastic neutron scattering [32, 31,
33]. There is, however, no experimental information about the dynamic
structure of the rest of $3d$ liquid metals.

To our knowlege, only l-Cu, l-Ni and l-Fe had been previously studied by other
groups using ab initio simulation methods. For l-Cu, three ab initio studies
evaluated its static properties at melting and some undercooled states [34,
36, 37]. As for l-Ni, Jakse et al. [38, 39] studied a range of static and
dynamic properties near melting and calculated the self-diffusion and shear
viscosity coefficients. In the case of l-Fe, Ganesh and Widom [40] analyzed
several static properties, getting also into the undercooled region.

In our recent studies of the liquid 3$d$ transition metals [1, 2, 3, 4, 5], we
were interested in the calculation and analysis of their dynamic properties.
Besides its intrinsic interest, we were also motivated by the recent finding
of some remarkable features such as low-energy transverse-like excitations in
the dynamic structure factor and/or the appearance of a high-frequency branch
in the transverse current dispersion relation, in addition to the usual one
that appears at low frequencies. Indeed, the study of the microscopic
processes behind the dynamic properties of liquid systems has been a long
standing research field of Prof. Mryglod, and in this respect we spotlight his
contributions to the development and application of the generalized collective
modes (GCM) approach. As for the above mentioned features, there is still no
clear explanation of their physical origin, although we have recently advanced
some connection between the existence of a second branch in the transverse
dispersion relation and the coupling of the transverse current with density
fluctuations at all wavevectors [4, 41].

Our studies of the bulk liquid 3$d$ transition metals were performed with an
ab initio molecular dynamics (AIMD) simulation method based on the density
functional theory [42, 43]. The liquid metal is modelled as an interacting
system of ions and electrons where the ionic positions evolve according to
classical mechanics while the electronic subsystem follows adiabatically.
Table 1 provides information concerning the thermodynamic states as well as
some technical details. The AIMD simulations were performed with the Quantum-
ESPRESSO package [44, 45], excepting l-Fe and l-Zn for which we used the VASP
code [46, 47, 48, 49]. The generalized gradient approximation of Perdew-Burke-
Ernzerhof [50] was used to account for the electronic exchange-correlation
energy, with the exceptions of l-Ti and l-Cu for which we used Perdew and
Wang’s approximation [51] and Perdew and Zunger’s local density expression
[52], respectively.

Table 1:  Input data of the liquid 3$d$ transition metals considered in the present AIMD simulation study. $\rho$ is the total ionic number density, (taken from [57]), $T$ is the temperature, $N_{\rm part}$ is the number of particles, $\Delta t$ is the ionic timestep, $Z_{\rm val}$ is the number of valence electrons and $N_{\text{c}}$ is the total number of configurations.  |  $\rho$ (Å-3) |  $T$ (K) | $N_{\rm part}$ |  $\Delta t$ (ps) | $Z_{\rm val}$ | $N_{\text{c}}$
---|---|---|---|---|---|---
Sc | 0.0391 | 1875 | 120 | 0.0045 | 11 | 29000
Ti | 0.0522 | 2000 | 100 | 0.0045 | 12 | 30700
V | 0.0634 | 2173 | 120 | 0.0040 | 13 | 20000
Cr | 0.0713 | 2173 | 100 | 0.0045 | 14 | 27000
Mn | 0.0655 | 1550 | 120 | 0.0045 | 15 | 13200
Fe | 0.0746 | 1873 | 100 | 0.0020 | 8 | 22100
Co | 0.0787 | 1850 | 135 | 0.0045 | 9 | 13100
Ni | 0.0792 | 1773 | 120 | 0 0055 | 10 | 21000
Cu | 0.0755 | 1423 | 150 | 0 0050 | 11 | 25500
Zn | 0.0616 | 723 | 90 | 0 0040 | 12 | 36000

The ion-electron interaction was described by ultrasoft pseudopotentials [53],
excepting l-Fe and l-Zn for which the projector augmented wave all-electron
description [54, 55] was used. For all the systems, the cutoff for the plane-
wave representation of the orbitals was within the range $25.0$–$35.0$ Ryd and
the sampling of the Brillouin zone was performed by means of the single
$\Gamma$ point.

We recall that most previous studies of the liquid 3$d$ transition metals were
carried out by using semiempirical or more fundamental pair (or many-body)
potentials coupled with liquid state theories or classical molecular dynamics
simulations [59, 58, 56, 57]. Moreover, in these approaches the electronic
degrees of freedom are hidden into the effective potential and, therefore, no
electronic properties (i.e., density of states, conductivity, …) can be
consistently/reliably calculated. On the contrary, the AIMD simulation method
yields parameter-free interatomic interactions derived from first principles,
and the valence electrons are explicitly taken into account during the
calculation.

##  2 Results and discussion

The evaluation of static and dynamic properties of the bulk liquid metal was
performed by using the total number of equilibrium configurations listed in
table 1.

###  2.1 Static properties

The calculated static structure factors, $S(q)$, are plotted in figure 1 along
with the corresponding WXD data [6]. The figure also includes the more recent
XD and ND data [7, 8, 9, 10, 11] for l-Ti, l-Fe, l-Ni and l-Zn.

Aside from l-Sc and l-Ti, the calculated $S(q)$ shows a good qualitative
agreement with the WXD data. There are some slight disparities, namely: (i) a
very small phase shift in l-V, l-Mn and l-Co, (ii) the calculations predict an
asymmetric shape for the second peak with a more or less marked shoulder,
whereas the XD data display, for all systems, a symmetric second peak. Notice
that shoulders on the high-$q$ side of the second peak of $S(q)$ have also
been experimentally observed in l-Ti, l-Fe, l-Ni and l-Zr [7, 8, 9, 10, 12];
moreover, they become more marked upon undercooling [60, 61, 62].

As for l-Ti, comparison with the recent XD and ND data [9, 10] for the $S(q)$,
reveals an excellent agreement with the present resuts, including the shape of
the second peak as well as the shoulder on its high-$q$ side. Consequently, we
believe that the WXD data for l-Ti are somewhat unreliable and that something
similar might also happen with l-Sc, where an analogous (although opposite)
dephasage is visible when compared with the present AIMD result.

Figure 1: (Colour online) Static structure factor, $S(q)$, of the liquid 3$d$
transition metals. Continuous line: present AIMD calculations. Open circles:
WXD data from [6]. The blue circles are XD data [9, 10, 7, 8, 11] and red ones
are ND data [9, 10, 7, 8]. The inset shows the corresponding pair distribution
function.

From the calculated $S(q)$, within the range $q\leqslant 1.2$ Å-1, we used a
least squares fit, $S(q)=s_{0}+s_{2}q^{2}$, to obtain an estimate for $S(q\to
0)$ and the results are given in table 2. Then, we evaluated the isothermal
compressibility, $\kappa_{T}$, by resorting to the relation $S(q\to 0)$ =
$\rho k_{\text{B}}T\kappa_{T}$, where $k_{\text{B}}$ is Boltzmann’s constant.
Table 2 lists the obtained results along with the available experimental data.

In figure 1, we also depicted the associated pair distribution functions,
$g(r)$, which are compared with their respective XD data [6]. The associated
coordination number (CN) was evaluated by integrating the radial distribution
function, $4\piup\rho r^{2}g(r)$, up to the position of its first minimum,
$r_{\rm min}$. Table 2 shows the obtained values for the CN’s which are
typical of simple liquid metals near their respective triple points [35].

Figure 2: (Colour online) Variation of the most abundant bonded pairs.

A more detailed three-dimensional description of the short range order in
these liquid metals is achieved by applying the common neighbour analysis [63]
(CNA) method. This method allows one to characterize the local environment
surrounding each atomic pair that contributes to the peaks of the $g(r)$, in
terms of the number and properties of the common nearest neighbours of the
pair under consideration. Each bonded pair of atoms is characterized by four
indices, with the first index being 1 if the pair belongs to the first peak of
$g(r)$, the second index stands for the number of common first neighbours, the
third index is the number of bonds that connect those shared first neighbours
and the fourth index discriminates among those configurations having the same
first three indices but a different topology. The CNA method allows one to
discern different local structures such as FCC, HCP, BCC and ICOS. For
example, the close-packed (CP) structures, FCC and HCP, are composed of
142x-type pairs, the BCC is typified by 144x and 166x pairs whereas the
perfect ICOS is characterized by 155x pairs and the distorted ICOS is
described by the 154x and 143x pairs.

Table 2:  Calculated values for $r_{\rm min}$ (in Å), coordination numbers CN, $S(q\to 0)$ and isothermal compressibilities $\kappa_{T}$ (in 10-11$m^{2}$$N^{-1}$ units) for the liquid transition metals at the thermodynamic states given in table 1. The numbers in parenthesis are experimental and/or semiempirical data.  | $r_{\rm min}\;\;\;$ | CN | $S(q\to 0)\;\;\;\;$ | $\kappa_{T}\;\;$
---|---|---|---|---
Sc | 4.27 | 12.8 | 0.020 $\pm$ 0.002 | 2.00 $\pm$ 0.10 (3.25a )
Ti | 3.90 | 13.3 | 0.018 $\pm$ 0.001 | 1.28 $\pm$ 0.10 (1.67 $\pm$ 0.02b )
V | 3.62 | 12.6 | 0.021 $\pm$ 0.002 | 1.10 $\pm$ 0.10 (1.43a, 1.21b )
Cr | 3.47 | 12.5 | 0.020 $\pm$ 0.002 | 0.94 $\pm$ 0.10 (1.00a, 1.10b )
Mn | 3.48 | 12.3 | 0.069 $\pm$ 0.003 | 4.9$\pm$ 0.30 (3.00a, 3.74 $\pm$ 0.02)
Fe | 3.42 | 12.5 | 0.024 $\pm$ 0.002 | 1.24 $\pm$ 0.10 (1.21 $\pm$ 0.02)
Co | 3.25 | 11.6 | 0.022 $\pm$ 0.002 | 1.10 $\pm$ 0.10 (1.43a, 1.18 $\pm$ 0.02b )
Ni | 3.27 | 11.8 | 0.018 $\pm$ 0.002 | 0.95 $\pm$ 0.10 (1.45a, 1.04 $\pm$ 0.02b )
Cu | 3.28 | 11.8 | 0.016 $\pm$ 0.002 | 1.10 $\pm$ 0.10 (1.49b )
Zn | 3.69 | 12.4 | 0.009 $\pm$ 0.004 | 1.52 $\pm$ 0.60 (2.34b )
a [25], b [13]

The CNA study is performed on inherent structures in which the ions are
brought to local minima of the potential energy surface. For each metal, this
calculation was made for four/five ionic configurations, well separated in
time ($\approx 15.0{-}20.0$ ps), and the results were averaged. These are
plotted in figure 2 where it is noticed that the five-fold symmetry dominates,
as the sum of perfect and distorted ICOS structures varies between $\approx$
59% (l-Sc) and $\approx$ 72% (l-Cr, l-Co) of the pairs. The amount of local
BCC-type pairs is also remarkable as it ranges from $\approx$ 32% (l-Ti) to
$\approx$ 6% (l-Zn) of the pairs. Finally, the CP-type pairs are almost absent
in the early transition metals but become significant for Mn and for the late
elements (Co to Zn). This behaviour of the CP vs BCC local structures
correlates with the corresponding phase diagrams of the elements, since Sc to
Fe melt from a BCC phase, while Co to Zn do so from CP structures, either FCC
(Co, Ni, Cu) or HCP (Zn). The phase diagram of Mn is rather complex, and the
high temperature BCC phase only exists down to 100 K below the melting point,
where it transforms into FCC, and this latter phase appears to somehow
partially survive locally also after melting.

###  2.2 Dynamic properties

We evaluated several dynamic magnitudes and their associated time correlation
functions. Due to the macroscopically isotropic behaviour of the fluid, those
correlation functions with a $\vec{q}$-dependence are transformed into
$|\vec{q}|$-dependent magnitudes.

####  2.2.1 Single-particle dynamics

We report results for the self-diffusion coefficient, $D$, which, along with
the shear viscosity, plays an important role in the study of crystal
nucleation and growth in metallic melts. We evaluated the velocity
autocorrelation function, $Z(t)$, and the mean square displacement, $\delta
R^{2}(t)$, of a tagged ion in the fluid, for all the $3d$ liquid metals. The
diffusion coefficient can be determined through both of these functions, and
both approaches yielded practically the same result for $D$, which is given in
table 3 along with the available experimental/semiempirical data.

Table 3:  Calculated values of the self-diffusion coefficient, $D$ (Å2/ps), adiabatic sound velocity $c_{\text{s}}$ (m/s) and shear viscosity $\eta$ (GPa$\,\cdot\,$ps) for the liquid transition metals at the thermodynamic states given in table 1. The numbers in parenthesis are experimental and/or semiempirical data.  | $D\;\;\;$ | $c_{\text{s}}\;\;\;$ | $\eta\;\;\;\;$
---|---|---|---
Sc | 0.55 $\pm$ 0.01 | 4400 $\pm$ 200 |  2.10 $\pm$ 0.10
Ti | 0.49 $\pm$ 0.01 (0.53 $\pm$ 0.02a) | 4640 $\pm$ 100 (4407 $\pm$ 108b) |  2.85 $\pm$ 0.15 (2.20c)
V | 0.51 $\pm$ 0.01 | 4725 $\pm$ 200 (4742b) |  3.30 $\pm$ 0.20
Cr | 0.47 $\pm$ 0.01 | 4500 $\pm$ 200 (4298b ) |  3.70 $\pm$ 0.20
Mn | 0.46 $\pm$ 0.02 | 3100 $\pm$ 150 (3381d) |  2.50 $\pm$ 0.20
Fe | 0.37 $\pm$ 0.02 (0.355e) | 3950 $\pm$ 150 (3820 $\pm$ 150f) |  5.00 $\pm$ 0.30 (5.30g)
Co | 0.60 $\pm$ 0.01 | 3700 $\pm$ 150 (4048d) |  3.20 $\pm$ 0.20 (4.12h, 3.60i)
Ni | 0.38 $\pm$ 0.01 (0.380 $\pm$ 0.006j) | 4250 $\pm$ 150 (4040 $\pm$ 150g) |  4.10 $\pm$ 0.20 (4.37h)
Cu | 0.27 $\pm$ 0.01 (0.37 $\pm$ 0.006l) | 3550 $\pm$ 150 (3481 d) |  4.80 $\pm$ 0.20 (4.38m)
Zn | 0.23 $\pm$ 0.01 (0.24 $\pm$ 0.006n) | 3150 $\pm$ 250 (2849 d) |  2.40 $\pm$ 0.30 (3.50e,p)
a [17], b [14], c [18], d [13], e [19], f [20], g [21], h [15], i [16], j
[32], l [33], m [22], n [23], p [24]

####  2.2.2 Collective dynamics

In a liquid, the dynamics of density fluctuations is usually described by the
intermediate scattering function, $F(q,t)$, which is defined as
autocorrelation function of the microscopic wavevector-dependent density [35].
The Fourier transform (FT) of $F(q,t)$ into the frequency domain leads to the
dynamic structure factor, $S(q,\omega)$, which can be directly measured by
either INS and/or IXS experiments.

Another important collective magnitude is the current due to the overall
motion of the particles. It is a vectorial magnitude which is usually split
into its longitudinal and transverse components with respect to $\vec{q}$.
Then, the longitudinal, $J_{\text{L}}(q,t)$, and transverse
$J_{\text{T}}(q,t)$, current correlation functions are obtained as
autocorrelation functions of the corresponding quantities [35]. Their time FT
yields the associated spectra, $J_{\text{L}}(q,\omega)$ and
$J_{\text{T}}(q,\omega)$, respectively.

Figure 3: (Colour online) Intermediate scattering function, $F(q,t)/S(q)$, of
l-Mn at $T=1550$ K and l-Cr at $T=2173$ K for several $q$ values.

Figure 4: (Colour online) Same as the previous figure but for l-Co at $T=1850$
K and l-Cu at $T=1423$ K.

Figure 5: (Colour online) Dynamic structure factor, $S(q,\omega)/S(q)$, of
l-Mn at $T=1550$ K and l-Cr at $T=2173$ K for several $q$ values.

In figures 3–4, we plotted, for some metals and for some range of $q$-values,
the obtained AIMD results for the $F(q,t)$. For small $q$-values, the $F(q,t)$
show an oscillatory behaviour, indicative of wave propagation, which is
superposed on another decaying component, indicative of relaxation modes. The
oscillations become weaker with increasing $q$-values and disappear at around
$q\approx 4/5q_{p}$, where $q_{p}$ denotes the position of the main peak of
$S(q)$. For larger $q$ values, the decaying component becomes dominant leading
at $q\approx q_{p}$ to a very slow monotonous decay of the $F(q,t)$. The
relative strength of both components depends on the system and we observe that
l-Mn has the weaker oscillations superposed on a very strong diffusive
component.

Figures 5–6 show the associated $S(q,\omega)$. The behaviour of the calculated
$S(q,\omega)$ is qualitatively similar for all the 3$d$ metals studied in this
paper. Specifically, the $S(q,\omega)$ show side-peaks up to $q$-values
$\approx(3/5)q_{p}$, which thereafter evolve into shoulders lasting up to
$q\approx(4/5)q_{p}$. Therefrom, the $S(q,\omega)$ show a monotonous
decreasing behaviour. Notice that for small $q$-values, the $S(q,\omega)$ of
l-Mn takes relatively large values when $\omega\to 0$; this is due to the
important diffusive component of its $F(q,t)$ at those small $q$-values. From
the positions of the side-peaks, $\omega_{m}(q)$, in the $S(q,\omega)$, a
dispersion relation for the density excitations was obtained and is plotted in
figure 7. In the hydrodynamic region (small $q$), the slope of the dispersion
relation curve is the $q$-dependent adiabatic sound velocity
$c_{\text{s}}(q)=v_{\text{th}}\sqrt{\gamma/S(q)}$, with
$v_{\text{th}}=\sqrt{k_{\text{B}}T/m}$ being the thermal velocity and $\gamma$
being the ratio of specific heats. In the $q\rightarrow 0$ limit, the
$c_{\text{s}}(q)$ reduces to the bulk adiabatic sound velocity $c_{\text{s}}$.
However, the small size of the simulation box implies that the smallest
attainable $q$ value, namely $q_{\rm min}$, is not small enough so as to
permit an accurate determination of $c_{\text{s}}$; nevertheless, an
approximate estimation can be obtained from value of $\omega_{m}(q_{\rm
min})/q_{\rm min}$. From the values $q_{\rm min}$= 0.43 Å-1 (Sc), 0.506 Å-1
(Ti), 0.508 Å-1 (V), 0.561 Å-1 (Cr), 0.513 Å-1 (Mn), 0.570 Å-1 (Fe), 0.525 Å-1
(Co), 0.547 Å-1 (Ni), 0.500 Å-1 (Cu), and 0.553 Å-1 (Zn), we obtained the
values given in table 3 where they are compared with the available
experimental/semiempirical data.

Figure 6: (Colour online) Same as the previous figure but for l-Co at $T=1850$
K and l-Cu at $T=1423$ K.

We also evaluated the longitudinal and transverse currents,
$J_{\text{L}}(q,t)$ and $J_{\text{T}}(q,t)$, as well as their respective
spectra, $J_{\text{L}}(q,\omega)$ and $J_{\text{T}}(q,\omega)$. The obtained
$J_{\text{L}}(q,\omega)$ show one peak for each $q$-value and the frequencies
associated to those peaks, $\omega_{\text{L}}(q)$, constitute the longitudinal
dispersion relation for the associated collective modes. Figure 7 represents
the obtained dispersion relations, which show the typical behaviour found in
other liquid systems [35].

Figure 7: (Colour online) Dispersion relations for several 3$d$ liquid metals.
Red diamonds and stars: longitudinal dispersion obtained from the AIMD results
for the positions of the inelastic peaks in the $S(q,\omega)$ and from the
maxima in the spectra of the longitudinal current, $J_{\text{L}}(q,\omega)$,
respectively. Open circles and triangles: transverse dispersion from the
positions of the peaks in the spectra $J_{\text{T}}(q,\omega)$. The green
lozenges are the positions of the TA excitation modes found in the calculated
dynamic structure factors, $S(q,\omega)$.

As for the $J_{\text{T}}(q,\omega)$, it may exhibit peaks within some
$q$-range, which are related to propagating shear waves. For all the metals
studied in this paper, the associated $J_{\text{T}}(q=q_{\rm min},\omega)$
already showed a peak and with increasing $q$-values the associated frequency
of the peak, $\omega_{\text{T}}(q)$ also increased and reached a maximum value
for $q\approx(2/3)q_{p}$; therefrom it decreased and vanished at $\approx
3.0q_{p}$. An example of this behaviour is given in figure 8 which shows the
calculated $J_{\text{T}}(q,\omega)$ for l-V and l-Cr, as a function of
$\omega$, and for a range 0 $<q\leqslant 4.0$ Å-1.

Moreover, for some metals (l-V, l-Cr, l-Co, l-Cu, l-Ni and l-Zn) we found that
the $J_{\text{T}}(q,\omega)$ exhibited, within the $q$-region 0.8 $\leqslant
q/q_{p}\leqslant 1.2$, another peak with a higher frequency. Figure 8 shows,
for l-V and l-Cr, a closer view of that specific region where the two maxima
are visible. For the other metals, (i.e., l-Sc, l-Ti, l-Mn and l-Fe) the
results are not conclusive because, besides the lower frequency peak, we also
found higher frequency shoulders instead of peaks.

Figure 8: (Colour online) Transverse current spectra,
$J_{\text{T}}(q,\omega)$, of l-V at $T=2173$ K and l-Cr at $T=2173$ K at
several $q$-values. The inset gives a detailed view of the region where the
$J_{\text{T}}(q,\omega)$ displays two maxima.

The appearance of another, high frequency branch in the transvese dispersion
relation, had been reported first in some liquid metals at high pressure,
i.e., Li, Na, Fe, Al and Pb, [2, 65, 66, 64]. Subsequently, it was also found
in l-Tl and l-Pb [67, 66, 64] at ambient pressure.

From the calculated $J_{\text{T}}(q,t)$, we also estimated [35] the shear
viscosity coefficient $\eta$, as follows. The memory function representation
of $J_{\text{T}}(q,t)$, namely

$\tilde{J_{\text{T}}}(q,z)=\frac{k_{\text{B}}T}{m}\left[z+\frac{q^{2}}{\rho
m}\tilde{\eta}(q,z)\right]^{-1},$ (2.1)

where the tilde denotes the Laplace transform, introduces a generalized shear
viscosity coefficient $\tilde{\eta}(q,z)$. The area under the normalized
$J_{\text{T}}(q,t)$ gives $m\tilde{J_{\text{T}}}(q,z=0)/(k_{\text{B}}T)$, from
which $\tilde{\eta}(q,z=0)\equiv\tilde{\eta}(q)$ can be obtained.
Extrapolation to $q\rightarrow 0$ gives the standard shear viscosity
coefficient $\eta$. Different extrapolation functions lead to different values
of $\eta$, which are portrayed in the uncertainty assigned to the results.
Table 3 shows these results along with the available
experimental/semiempirical data.

Recently, it was suggested [68] that transverse-like low-energy excitations
might be observed as weak shoulders located in the region between the
quasielastic and the longitudinal inelastic peaks of the $S(q,\omega)$. These
excitations are usually visible within a small $q$-range around $q_{p}/2$,
because for smaller/greater $q$-values they are overcome by the
quasielastic/inelastic peaks. These transverse-like excitations were first
detected in the IXS spectra of l-Ga [68] at 313 K and subsequently in l-Cu,
l-Sn, l-Na, l-Fe, l-Ni and l-Zn [3, 69, 70, 29, 30, 71].

Figure 9: (Colour online) AIMD calculated dynamic structure factors
$S(q,\omega)/S(q)$ of the 3$d$ liquid transition metals. They are plotted for
different values of $q/q_{p}$ (top to bottom). The vertical scales are offset
for clarity. The arrows point to the locations of the transverse-like
excitations.

We analyzed the results for $S(q,\omega)$ and found that within the range 1.00
Å-1 $\leqslant q\leqslant$ 1.45 Å-1, some weak shoulders appear in the
$\omega$-region located between the quasielastic and the inelastic peaks. This
is shown in figure 9; moreover, the energies associated to these shoulders are
close to those corresponding to the peaks in the transverse current spectra.
The real physical origin of these excitations is still a matter of debate, and
some alternative interpretations, for instance heat waves, were proposed
within the framework of the GCM approach [72].

###  2.3 Electronic properties: density of states

The partial and total electronic density of states, $n(E)$, were obtained from
the self-consistently determined eigenvalues; they were averaged over four/six
ionic configurations well separated in time ($\approx 15.0$ ps), and the
sampling of the Brillouin zone was performed using an $8\times 8\times 8$
Monkhorst-Pack set.

Figure 10: Total electronic density of states. Notice the different scale in
the upper and lower panels.

Figure 10 shows the obtained results for the electronic total $n(E)$
associated to the outer valence electrons, which is dominated by the 3$d$
states, except in the early transition elements, where the $d$ band occupation
is low. Notice that for the metals Sc to Mn we are not showing the
contribution from the $3s$ and $3p$ atomic states, which were included in the
self-consistent calculation but lie way below the Fermi energy. In the figure
we observe the progresive filling of the $d$ band, which is complete at Cu,
together with a decrease of its width and an increase of its height as the
occupation rises.

##  3 Conclusions

An ab initio molecular dynamics simulation method was used to calculate
several static and dynamic properties of the bulk liquid $3d$ transition
metals near their respective melting points.

The results for the static structure, as characterized by the static structure
factor $S(q)$, show a good agreement with the available experimental data. The
only exception is l-Sc where a clear difference is observed in the phase of
the oscillations, but this may be due to inaccuracies in the old experimental
data, similar to the case of the old data for l-Ti. The $S(q)$ display, in
most cases, an asymmetric shape in the second peak which qualitatively agrees
with the recent ND and/or XD data. This feature is related to the appearance
of icosahedral short-range order in the liquid metal. This is confirmed by a
more detailed study of the liquid static structure by using the CNA method,
which revealed a marked abundance of five-fold structures. On the other hand,
secondary local structures also present in the liquid correlate rather well
with the solid phase from where the metals melt.

The calculated dynamic structure factors, $S(q,\omega)$, show side-peaks which
are indicative of collective density excitations. A detailed analysis of the
$S(q,\omega)$ revealed some type of excitations which have similar features as
the transverse-like excitation modes found in IXS and INS data for several
liquid metals near their respective melting points.

From the calculated longitudinal and transverse current correlation functions,
we evaluated the associated spectral functions and the corresponding
dispersion relations. For some metals, we obtained two branches of transverse
modes, with the high frequency branch appearing over a limited $q$-range which
is always located in the second pseudo-Brillouin zone.

Results are also reported for several transport coefficients such as the self-
diffusion, adiabatic sound velocity and shear viscosity. Although the
calculated values for the adiabatic sound velocity show a reasonable agreement
with the available experimental data, the lack of experimental data for the
self-diffusion and the shear viscosity provides a further relevance of these
results.

## Acknowledgements

This work has been supported by Junta de Castilla y Leon (project VA124G18).
LEG and DJG also acknowledge the funding from the Spanish Ministry of Economy
and Competitiveness (Project PGC2018-093745-B-I00), which is partially
supported by the EU FEDER program.

## References

  * [1] Del Rio B.G., Rodriguez O., González L.E., González D.J., Comput. Mater. Sci., 2017, 139, 243,
doi:10.1016/j.commatsci.2017.07.027.

  * [2] Marques M., González L.E., González D.J., Phys. Rev. B, 2015, 92, 134203, doi:10.1103/PhysRevB.92.134203.
  * [3] Del Rio B.G., González D.J., González L.E., J. Chem. Phys., 2017, 146, 034501, doi:10.1063/1.4973803.
  * [4] Del Rio B.G., González L.E., Phys. Rev. B, 2017, 95, 224201, doi:10.1103/PhysRevB.95.224201.
  * [5] Del Rio B.G., Pascual C., González D.J., González L.E., J. Phys.: Condens. Matter, 2020, 32, 214005,
doi:10.1088/1361-648X/ab6f16.

  * [6] Waseda Y., The Structure of Non-Crystalline Materials, McGraw-Hill, New York, 1980.
  * [7] Schenk T., Holland-Moritz D., Simonet V., Bellisent R., Herlach D.M., Phys. Rev. Lett., 2002, 89, 075507,
doi:10.1103/PhysRevLett.89.075507.

  * [8] Inui M., Maruyama K., Kajihara Y., Nakada M., Phys. Rev. B, 2009, 80, 180201,
doi:10.1103/PhysRevB.80.180201.

  * [9] Lee G.W., Gangopadhyay A.K., Kelton K.F., Hyers R.W., Rathz T.J., Rogers J.R., Robinson D.S., Phys. Rev. Lett., 2004, 93, 037802, doi:10.1103/PhysRevLett.93.037802.
  * [10] Holland-Moritz D., Heinen O., Bellissent R., Schenk T., Mater. Sci. Eng., A, 2007, 449–451, 42,
doi:10.1016/j.msea.2005.12.093.

  * [11] Lou H., Wang X., Cao Q., Zhang D., Zhang J., Hu T., Mao H.-K., Jiang J.-Z., Proc. Natl. Acad. Sci. U.S.A., 2013, 110, 10068, doi:10.1073/pnas.1307967110.
  * [12] Lee G.W., Gangopadhyay A.K., Hyers R.W., Rathz T.J., Rogers J.R., Robinson D.S., Goldman A.I., Kelton K.F., Phys. Rev. B, 2008, 77, 184102, doi:10.1103/PhysRevB.77.184102.
  * [13] Blairs S., Int. Mater. Rev., 2007, 52, 321, doi:10.1179/174328007X212490.
  * [14] Casas J., Keita N.M., Steinemann S.G., Phys. Chem. Liq., 1984, 14, 155, doi:10.1080/00319108408080806.
  * [15] Assael M.J., Kakosimos K., Banish R.M., Brillo J., Egry I., Brooks R., Quested P.N., Mills K.C., Nagashima A., Sato Y., Wakeham W.A., J. Phys. Chem. Ref. Data, 2006, 35, 285, doi:10.1063/1.2149380
  * [16] Egry I., Scr. Metall. Mater., 1993, 28, 1273, doi:10.1016/0956-716X(93)90467-7.
  * [17] Said A.H., Sinn H., Alatas A., Burns C.A., Price D.L., Saboungi M.L., Schirmacher W., Phys. Rev. B, 2006, 74, 172202, doi:10.1103/PhysRevB.74.172202.
  * [18] Agaev A.D., Kostikov V.I., Bobkovskii V.N., Izv. Akad. Nauk S.S.S.R., Ser. Met., 1980, 3, 43.
  * [19] Iida T., Guthrie R.I.L., The Thermophysical Properties of Metallic Liquids, Oxford University Press, Oxford, 2015.
  * [20] Nash P.M., Manghnani M.H., Secco R.A., J. Geophys. Res. Solid Earth, 1994, 99, 4285, doi:10.1029/93JB03111.
  * [21] Shimoji M., Itami T., Atomic Transport in Liquid Metals, Trans-Tech Publications, Brookfield, 1986.
  * [22] Hirai M., ISIJ Int., 1993, 33, 251, doi:10.2355/isijinternational.33.251.
  * [23] Nachtrieb N.H., Fraga E., Wahl C., J. Phys. Chem., 1963, 67, 2353, doi:10.1021/j100805a022.
  * [24] Assael M.J., Armyra I.J., Brillo J., Stankus S.V., Wu J., Wakeham W.A., J. Phys. Chem. Ref. Data, 2012, 41, 033101, doi:10.1063/1.4729873.
  * [25] Marcus Y., J. Chem. Thermodyn., 2017, 109, 11, doi:10.1016/j.jct.2016.07.027.
  * [26] Hosokawa S., Inui M., Matsuda K., Ishikawa D., Baron A.Q.R., Phys. Rev. B, 2008, 77, 174203,
doi:10.1103/PhysRevB.77.174203.

  * [27] Johnson M.W., McCoy B., March N.H., Page D.I., Phys. Chem. Liq., 1977, 6, 243,
doi:10.1080/00319107708084143.

  * [28] Bermejo F.J., Saboungi M.L., Price D.L., Alvarez M., Roessli B., Cabrillo C., Ivanov A., Phys. Rev. Lett., 2000, 85, 106, doi:10.1103/PhysRevLett.85.106.
  * [29] Hosokawa S., Inui M., Kajihara Y., Tsutsui S., Baron A.Q.R., J. Phys.: Condens. Matter, 2015, 27, 194104,
doi:10.1088/0953-8984/27/19/194104.

  * [30] Zanatta M., Sacchetti F., Guarini E., Orecchini A., Paciaroni A., Sani L., Petrillo C., Phys. Rev. Lett., 2015, 114, 187801, doi:10.1103/PhysRevLett.114.187801.
  * [31] Horbach J., Rozas R.E., Unruh T., Meyer A., Phys. Rev. B, 2009, 80, 212203, doi:10.1103/PhysRevB.80.212203.
  * [32] Chathoth S.M., Meyer A., Koza M.M., Juranyi F., Appl. Phys. Lett., 2004, 85, 4881, doi:10.1063/1.1825617.
  * [33] Meyer A., Phys. Rev. B, 2010, 81, 012102, doi:10.1103/PhysRevB.81.012102.
  * [34] Pasquarello A., Laasonen K., Car R., Lee C., Vanderbilt D., Phys. Rev. Lett., 1992, 69, 1982,
doi:10.1103/PhysRevLett.69.1982.

  * [35] Balucani U., Zoppi M., Dynamics of the Liquid State, Clarendon, Oxford, 1994.
  * [36] Kresse G., Hafner J., Phys. Rev. B, 1993, 48, 13115, doi:10.1103/PhysRevB.48.13115.
  * [37] Jakse N., Pasturel A., Condens. Matter Phys., 2015, 18, 43603, doi:10.5488/CMP.18.43603.
  * [38] Jakse N., Pasturel A., J. Chem. Phys., 2004, 120, 6124, doi:10.1063/1.1651054.
  * [39] Jakse N., Wax J.F., Pasturel A., J. Chem. Phys., 2007, 126, 234508, doi:10.1063/1.2741521.
  * [40] Ganesh P., Widom M., Phys. Rev. B, 2008, 77, 014205, doi:10.1103/physrevb.77.014205
  * [41] Del Rio B.G., Chen M., González L.E., Carter E.A., J. Chem. Phys., 2018, 149, 094504, doi:10.1063/1.5040697.
  * [42] Hohenberg P., Kohn W., Phys. Rev., 1964, 136, B864, doi:10.1103/PhysRev.136.B864
  * [43] Kohn W., Sham L.J., Phys. Rev., 1965, 140, A1133, doi:10.1103/PhysRev.140.A1133.
  * [44] Giannozzi P., Baroni S., Bonini N., Calandra M., Car R., Cavazzoni V., Ceresoli D., Chiarotti G.L., Cococcioni M., Dabo I., J. Phys.: Condens. Matter, 2009, 21, 395502, doi:10.1088/0953-8984/21/39/395502.
  * [45] Giannozzi P., Andreussi O., Brumme T., Bunau O., Buongiorno Nardelli M., Calandra M., Car R., Cavazzoni C., Ceresoli D., Cococcioni M., J. Phys.: Condens. Matter, 2017, 29, 465901, doi:10.1088/1361-648X/aa8f79.
  * [46] Kresse G., Hafner J., Phys. Rev. B, 1993, 47, 558(R), doi:10.1103/PhysRevB.47.558.
  * [47] Kresse G., Hafner J., Phys. Rev. B, 1994, 49, 14251, doi:10.1103/physrevb.49.14251
  * [48] Kresse G., Furthmuller J., Comput. Mater. Sci., 1996, 6, 15, doi:10.1016/0927-0256(96)00008-0.
  * [49] Kresse G., Furthmuller J., Phys. Rev. B, 1996, 54, 11169, doi:10.1103/PhysRevB.54.11169.
  * [50] Perdew J.P., Burke K., Ernzerhof M., Phys. Rev. Lett., 1996, 77, 3865, doi:10.1103/PhysRevLett.77.3865.
  * [51] Perdew J.P., Wang Y., Phys. Rev. B, 1992, 45, 13244, doi:10.1103/PhysRevB.45.13244.
  * [52] Perdew J.P., Zunger A., Phys. Rev. B, 1981, 23, 5048, doi:10.1103/PhysRevB.23.5048.
  * [53] Vanderbilt D., Phys. Rev. B, 1990, 41, 7892, doi:10.1103/PhysRevB.41.7892.
  * [54] Blöchl P.E., Phys. Rev. B, 1994, 50, 17953, doi:10.1103/PhysRevB.50.17953.
  * [55] Kresse G., Joubert D., Phys. Rev. B, 1999, 59, 1758, doi:10.1103/PhysRevB.59.1758.
  * [56] Bhuiyan G.M., Bretonnet J.L., Gonzalez L.E., Silbert M., J. Phys.: Condens. Matter, 1992, 4, 7651,
doi:10.1088/0953-8984/4/38/002.

  * [57] Bhuiyan G.M., Silbert M., Stott M.J., Phys. Rev. B, 1996, 53, 636, doi:10.1103/PhysRevB.53.636.
  * [58] Hausleitner C., Kahl G., Hafner J., J. Phys.: Condens. Matter, 1991, 3, 1589, doi:10.1088/0953-8984/3/11/017.
  * [59] Todd J.R., Brown J.S., Phys. Lett. A, 1976, 59, 302, doi:10.1016/0375-9601(76)90801-X.
  * [60] Kim T.H., Kelton K.F., J. Chem. Phys., 2007, 126, 054513, doi:10.1063/1.2431173.
  * [61] Jakse N., Pasturel A., Phys. Rev. Lett., 2003, 91, 195501, doi:10.1103/PhysRevLett.91.195501
  * [62] Jakse N., Le Bacq O., Pasturel A., J. Non-Cryst. Solids, 2007, 353, 3684, doi:10.1016/j.jnoncrysol.2007.05.131.
  * [63] Honeycutt J.D., Andersen H.C., J. Phys. Chem., 1987, 91, 4950, doi:10.1021/j100303a014.
  * [64] Bryk T., Jakse N., J. Chem. Phys., 2019, 151, 034506, doi:10.1063/1.5099099.
  * [65] Bryk T., Ruocco G., Scopigno T., Seitsonen A.P., J. Chem. Phys., 2015, 143, 104502, doi:10.1063/1.4928976.
  * [66] Bryk T., Demchuk T., Jakse N., Phys. Rev. B, 2019, 99, 014201, doi:10.1103/PhysRevB.99.014201.
  * [67] Bryk T., Demchuk T., Jakse N., Wax J.F., Front. Phys., 2018, 6, 00006, doi:10.3389/fphy.2018.00006.
  * [68] Hosokawa S., Inui M., Kajihara Y., Matsuda K., Ichitsubo T., Pilgrim W.C., Sinn H., González L.E., González D.J., Tsutsui S., Baron A.Q.R., Phys. Rev. Lett., 2009, 102, 105502,
doi:10.1103/PhysRevLett.102.105502.

  * [69] Hosokawa S., Munejiri S., Inui M., Kajihara Y., Pilgrim W.C., Baron A.Q.R., Shimojo F., Hoshino K., AIP Conf. Proc., 2013, 1518, 695, doi:10.1063/1.4794661.
  * [70] Hosokawa S., Munejiri S., Inui M., Kajihara Y., Pilgrim W.C., Ohmasa Y., Sutsui T., Baron A.Q.R., Shimojo F., Hoshino K., J. Phys.: Condens. Matter, 2013, 25, 112101, doi:10.1088/0953-8984/25/11/112101.
  * [71] Giordano V.M., Monaco G., Proc. Natl. Acad. Sci. U.S.A., 2010, 107, 21985, doi:10.1073/pnas.1006319107.
  * [72] Bryk T., Wax J.F., J. Chem. Phys., 2016, 144, 194501, doi:10.1063/1.4948709.

Ukrainian 0

Âèçíàåííÿ ç ïåðøèõ ïðèíöèïâ äåÿêèõ ñòàòèíèõ äèíàìíèõ âëàñòèâîñòåé 3$d$ ðäêèõ
ïåðåõäíèõ ìåòàëâ ïîáëèçó ïëàâëåííÿ Á.Ã. äåëü Ðo, Ê. Ïàñêóàëü, Î. Ðîäðãåñ, Ë..
Ãîíçàëåñ, Ä.Äæ. Ãîíçàëåñ

Âääë òåîðåòèíî¿ ôçèêè, óíâåðñèòåò Âàëüÿäîëäà, 47011 Âàëüÿäîëä, ñïàíÿ