# Application of Explainable Machine Learning in Detecting and Classifying
Ransomware Families Based on API Call Analysis
††thanks: This work is funded by NetApp.

Rawshan Ara Mowri  Department of Computer Science
North Carolina A&T State University
Greensboro, USA
<EMAIL_ADDRESS>Madhuri Siddula  Department of Computer Science
North Carolina A&T State University
Greensboro, USA
<EMAIL_ADDRESS>Kaushik Roy  Department of Computer Science
North Carolina A&T State University
Greensboro, USA
<EMAIL_ADDRESS>
###### Abstract

Ransomware has appeared as one of the major global threats in recent days. The
alarming increasing rate of ransomware attacks and new ransomware variants
intrigue the researchers to constantly examine the distinguishing traits of
ransomware and refine their detection strategies. Application Programming
Interface (API) is a way for one program to collaborate with another; API
calls are the medium by which they communicate. Ransomware uses this strategy
to interact with the OS and makes a significantly higher number of calls in
different sequences to ask for taking action. This research work utilizes the
frequencies of different API calls to detect and classify ransomware families.
First, a Web-Crawler is developed to automate collecting the Windows Portable
Executable (PE) files of 15 different ransomware families. By extracting
different frequencies of 68 API calls, we develop our dataset in the first
phase of the two-phase feature engineering process. After selecting the most
significant features in the second phase of the feature engineering process,
we deploy six Supervised Machine Learning models: Naïve Bayes, Logistic
Regression, Random Forest, Stochastic Gradient Descent, K-Nearest Neighbor,
and Support Vector Machine. Then, the performances of all the classifiers are
compared to select the best model. The results reveal that Logistic Regression
can efficiently classify ransomware into their corresponding families securing
99.15% overall accuracy. Finally, instead of relying on the ‘Black box’
characteristic of the Machine Learning models, we present the post-hoc
analysis of our best-performing model using ’SHapley Additive exPlanations’ or
SHAP values to ascertain the transparency and trustworthiness of the model’s
prediction.

###### Index Terms:

Ransomware Classification, Machine Learning, Explainable AI, Cyber Security

##  I Introduction

Recently, ransomware has become one of the biggest global challenges that are
agitating peoples’ normal lives. Being harmful software, it applies symmetric
and asymmetric cryptography to inscribe user information and poses a Denial-
of-Service (DoS) attack on the intended user [1]. The unique functional
process of ransomware attacks makes it more harmful than any malware attacks
and causes irreversible losses. According to [2], Fig. 1 illustrates the
number of publicized ransomware attacks in 2021, with inflation of 25% than
the same time in the previous year. Although the report does not include the
number of supply chain attacks, it is creating a big interference in providing
healthcare, purchasing groceries, and even loading fuel in vehicles. Examples
of these attacks are the Kaseya attack, the colonial pipeline attack, etc. In
addition, in the first six months of 2021, the FBI’s Internet Crime Complaint
Center documented 2,084 ransomware attacks [3], and the U.S. Treasury’s
Financial Crimes Enforcement Network (FinCEN) recorded the cost of around $590
million related to ransomware activities during that period [4]. Moreover,
distinct ransomware variants are being detected regularly and more than 130
different ransomware variants have been identified from 2020 till this year
causing an inevitable disturbance in day-to-day lives. [5].

Figure 1: Number of worldwide ransomware attacks in different sectors in 2021
(*Till November 2021)

Due to the increasing number of ransomware variants and ransomware attacks,
researchers have been earnestly involving themselves to look for efficient
ways to improve the scenarios. While some researchers are analyzing the
distinctive behaviors of ransomware by executing it in a secure environment
called Dynamic Analysis [6]-[11], some researchers are analyzing the
ransomware without any execution, referred to as Static Analysis [12], [13].
However, a good number of researchers are combining these two approaches and
adopting a Hybrid Analysis Approach [14], [15]. In this research, we have
opted for the dynamic analysis approach for its ability to detect and classify
ransomware based on behavioral patterns regardless of the code obfuscation
techniques deployed by the ransomware programmers [16], [17]. The main
contributions of this paper are:

  * •

Develop a Web-Crawler, ‘GetRansomware’ to automate collecting the Windows
Portable Executable (PE) files of 15 different ransomware families from the
ransomware repository. The Web-Crawler is essential to automate searching and
downloading the samples and to cut down the manual workload, but no prior
works targeted this scenario.

  * •

Develop our dataset and conduct feature selection through a two-phase feature
engineering process that includes- ‘Feature Extraction’ from the sample
binaries, and ‘Feature Selection’ to select the most important features for
each ML classifier.

  * •

Develop, evaluate and compare the performance of six State-of-the-art
Supervised Machine Learning models. Our approach includes utilizing Recursive
Feature Elimination with Cross-Validation (RFECV) for selecting the
significant features and RandomSearchCV for selecting the optimum
hyperparameter values for each ML classifier. Thereby we attempt to optimize
each model’s performance before the comparison is made.

  * •

Present the post-hoc analysis of the best-performing model using ‘SHapley
Additive exPlanations’ or SHAP values to ascertain the transparency and
trustworthiness of the model’s prediction. This insight presents a better idea
about which features are more dominant in detecting and classifying the
ransomware families. While explainability has been widely presented in malware
detection scenarios, to the best of the authors’ knowledge, till today, no
prior works presented their models’ explainability that considered only the
ransomware families.

The rest of this paper is structured as follows: Section II discusses the
related works. Section III presents our proposed method. The experimental
results and discussion are made in Section IV. Section V presents our model’s
explainability. Section VI concludes the paper with the direction for future
works.

##  II Related Works

Most researchers prefer the dynamic analysis approach because it can delineate
the behaviors of the ransomware in a more explicit manner. Maniath et al. [6]
analyzed the API call sequence of 157 ransomware and presented an LSTM-based
ransomware detection method. Despite securing 96.67% accuracy, this work lacks
complete information about the ransomware families/variants, and the number of
benign software used for the experiment. VinayaKumar et al. [7] proposed an
MLP-based ransomware detection method focusing on API call frequency and
secured 100%, and 98% accuracy for binary and multi-class classification
respectively. However, they deployed a simple MLP network that failed to
distinguish CryptoWall and Cryptolocker ransomware. Z. Chen et al. [8] used
the API Call Flow Graph (CFG) generated from the extracted API sequence of 83
ransomware and 83 benign software. Regardless of securing 98.2% exactness
using the Logistic Regression model, the work is based on a smaller dataset
that includes only four ransomware families. Also, graph-similarity analysis
requires higher computational power that some systems may fail to provide.
Takeuchi et al. [9] used API call sequences extracted from 276 ransomware, and
312 benign files to identify zero-day ransomware attacks. Although the work
secured 97.48% accuracy by deploying the Support Vector Machine, the accuracy
of this work decreases while using standardized vector representation because
of the less diverse dataset. Using the Intel Pin Tool, Bae et al. [10]
extracted the API call sequences from 1000 ransomware, 900 malware, and 300
benign files. Their sequential process includes generating an n-gram sequence,
input vector, and Class Frequency Non-Class Frequency (CF-NCF) for every
sample before fitting their model. Regardless of obtaining 98.65% accuracy
using the Random Forest classifier, the model’s performance can be improved
with the help of deception-based techniques. Hwang et al. [11] analyzed the
API call sequence of 2507 ransomware and 3886 benign files. They used two
Markov chains, one for ransomware and another for benign software to capture
the API call sequence patterns. By using Random Forest, they compensate Markov
Chains and control FPR and FNR to achieve better performance. Despite securing
97.3% accuracy, their model produces high FPR that can be improved with the
help of signature-based techniques.

A good number of researchers chose the static analysis approach to detect
ransomware. Baldwin and Dehghantanha [12] analyzed the opcode characteristics
of 5 crypto-ransomware families and 350 benign samples. Their experiment
involved the WEKA AI toolset, and the experimental results showed an accuracy
of 96.5% while recognizing five crypto-ransomware families and benign software
by using the Support Vector Machine classifier. However, their work could be
improved by extending the dataset and extracting those groups of opcodes
identified during the evaluation of attribute selection. Zhang et al. [13]
analyzed the opcode-based characteristics of 1787 ransomware of 8 different
ransomware families and 100 benign software. Their technique included moving
opcode groupings to the N-gram sequence and afterward Term Frequency Inverse
Document Frequency (TF-IDF). Five ML classifiers were used with 10-fold cross-
validation among which the Random Forest classifier achieved the highest
91.43% exactness. However, their model could not distinguish Reveton,
CryptoWall, and Locky.

Some researchers adopted a hybrid analysis approach that combines the features
extracted from the dynamic and static analyses. Subedi et al. [14] used both
dynamic and static analysis on the library, assembly, and function calls.
Moreover, they came up with a new analysis tool, namely, CRSTATIC which was
deployed to build signatures that could classify ransomware families with the
help of reverse engineering. However, they analyzed only 450 samples of
ransomware for their experiment and their experimental results showed some
important correlations among various code components instead of detecting
ransomware. Ferrante et al. [15] proposed a crossover way to deal with android
ransomware recognition. While the dynamic detection strategy considered memory
utilization, system call insights, CPU, and network utilization, the static
identification technique utilized the recurrence of opcodes. They analyzed
2386 benign android applications and 672 applications containing ransomware
for their experiment and used 3 ML classifiers- J48, Naïve Bayes and Logistic
Regression. Although their experimental result showed 100% accuracy, their
test set for the hybrid model included only 9 ransomware applications that
were not detected in the static analysis phase but were detected in the
dynamic analysis.

Regardless of the relentless effort of the prior research works, we have found
some common gaps and challenges, such as a) manual data collection process, b)
not having a large or diverse dataset that focuses on both Crypto and Locker
types of ransomware, c) using the traditional sandbox environment that often
falls short providing in-depth analysis report (i.e., Cuckoo Sandbox), d)
modeling binary classification models to detect ransomware from malware or
benign software, and e) not presenting the models’ explainability. Therefore,
the goal of our research work lies in a) automating the data collection
process, b) constructing a balanced dataset consisting of the most potential
features extracted through a methodical dynamic analysis of both Crypto and
Locker types of ransomware, c) utilizing an efficient sandbox environment that
conducts large scale integration with other antivirus engines and provides in-
depth analysis reports (i.e., Falcon Sandbox), d) developing efficient and
faster Supervised Machine Learning model that can attain a higher degree of
accuracy in case of multi-class classification, and e) presenting the
explainability of our model to affirm the trustworthiness of the model’s
prediction.

##  III Proposed Method

Fig. 2 illustrates the process overview of our proposed method which consists
of three steps: Data Collection, Feature Engineering, and Classification.

Figure 2: Process overview of our proposed method

###  III-A Data Collection

We have developed a Web-Crawler, namely, ‘GetRansomware’ to automate searching
and downloading Windows Portable Executable (PE) files of ransomware from
VirusShare [18]. We have also shared our Web-Crawler on our GitHub repository
for public access [19]. About 95% of our sample ransomware was collected from
VirusShare using GetRansomware. The rest of the ransomware samples were
collected from theZoo [20] and Hybrid-Analysis.com [21]. Our samples included
both Crypto and Locker types of ransomware and every sample was downloaded as
a password-protected compressed file from the repositories. In addition, we
have collected benign software from the PORTABLE FREEWARE collection [28]
which includes portable executables of benign media players, file explorers,
file compression tools, file managers, image editors, etc. Table I presents
the description of our dataset.

Table I: Description of Dataset Ransomware | Number of samples
---|---
Cerber(c1) | 95
CryptoLocker (c2) | 93
CryptoWall (c3) | 97
Eris (c4) | 98
Hive (c5) | 100
Jigsaw (c6) | 93
Locky (c7) | 95
Maze (c8) | 100
Mole (c9) | 100
Sage (c10) | 100
Satan (c11) | 100
Shade (c12) | 97
TeslaCrypt (c13) | 96
Virlock (c14) | 93
WannaCry (c15) | 93
Benign (c0) | 100
Total | 1550

###  III-B Feature Engineering

The scarcity of the ransomware dataset is one of the major challenges that
hinder the research work in this area [22]. Therefore, we aim to construct our
dataset through a thorough feature engineering process to overcome this
challenge. The feature engineering step is composed of two phases. The phases
are:

  1. 1.

Phase 1: Feature Extraction

  2. 2.

Phase 2: Feature Selection

####  III-B1 Phase 1: Feature Extraction

From the wide range of distinct behavioral features, we have considered
utilizing API call frequencies for our classifiers. Application Programming
Interface (API) calls are made by the application or program running at a user
level to request system services or gain access to the system resources, such
as network connection, files, registry, processes, etc. The OS performs the
requested services by issuing these calls, and the outcomes are returned to
the caller user applications. Thus, API calls made by the ransomware program
allow the attackers to obtain the control of the system and perform malicious
activities. Since API calls provide a complete picture of a particular
program, analyzing API call behavior leads the researchers to better
understand the program’s behavior [23], [24]. Therefore, we have considered
extracting the API call frequencies by executing the Windows PE files of the
ransomware and benign software in a secure virtual environment and utilizing
them as distinguishing features for our classification models.

We have analyzed the PE files with the help of Hybrid-Analysis.com [21],
powered by the CrowdStrike Falcon Sandbox [25]. To automate submitting
malicious binaries, pull the analysis report after the analysis, and perform
advanced or required search queries on the database, Falcon Sandbox provides a
free, convenient, and efficient API key that one can obtain from an authorized
user account. For analysis, we have used our API key and Falcon Sandbox Python
API Connector- VxAPI wrapper [26] to automatically submit the binaries from
our system. After submission, Falcon Sandbox runs the binaries in a Virtual
Machine (VM) and captures the run-time behaviors as illustrated in Fig. 3.
Later, it shows the analysis results on the web interface.

Contrary to the prior research works where the analysis tasks were done using
the Cuckoo Sandbox [6], [7], [9], [11], we have analyzed the PE files using
the Falcon Sandbox that used a VM (Windows 7 64-bit) to execute the PE files.
Falcon Sandbox can analyze a wide variety of files with or without
dependencies, such as PE (i.e., .com, .exe, .dll, etc.), Office (i.e., .doc,
.docx, .ppt, .pptx, etc.), and PDF, etc. Also, Falcon Sandbox incorporates
many other services, such as VirusTotal, Thug honeyclient, OPSWAT
Metadefender, TOR, NSRL (Whitelist), Phantom, and a large number of antivirus
engines to provide an integrated and accurate analysis report compared to
other sandboxes. While executing the binaries, we have set run-time to the
maximum available duration in the Falcon Sandbox to deal with the delayed
execution techniques deployed by the attackers. The total time for the
analysis was (1550 PE files * 7 minutes) = 181 hours = 8 days approximately.
Next, we obtained the analysis report by using our API key. As we aim to focus
on the API call behaviors, we have only sorted and computed the frequency of
each API call. At the end of the Feature Extraction phase, we obtained our
dataset ’RansomClass’ consisting of the frequencies of 68 distinct API calls
commonly called by the 15 different ransomware families and benign software as
presented in Table II.

Figure 3: Block diagram of the PE file execution process Table II: List of API
Calls Name of the API Calls
---
1\. FindWindowExW | 2\. LdrGetDllHandle
3\. NtAdjustPrivilegesToken | 4\. NtAlertThread
5\. NtAllocateVirtualMemory | 6\. NtAlpcSendWaitReceivePort
7\. NtConnectPort | 8\. NtCreateEvent
9\. NtCreateFile | 10\. NtCreateKey
11\. NtCreateKeyEx | 12\. NtCreateMutant
13\. NtCreateSection | 14\. NtCreateThreadEx
15\. NtCreateUserProcess | 16\. NtDelayExecution
17\. NtDeleteValueKey | 18\. NtDeviceIoControlFile
19\. NtEnumerateKey | 20\. NtEnumerateValueKey
21\. NtFsControlFile | 22\. NtGetContextThread
23\. NtMapViewOfSection | 24\. NtNotifyChangeKey
25\. NtOpenDirectoryObject | 26\. NtOpenEvent
27\. NtOpenFile | 28\. NtOpenKey
29\. NtOpenKeyEx | 30\. NtOpenMutant
31\. NtOpenProcess | 32\. NtOpenProcessToken
33\. NtOpenSection | 34\. NtOpenThreadToken
35\. NtProtectVirtualMemory | 36\. NtQueryAttributesFile
37\. NtQueryDefaultLocale | 38\. NtQueryDirectoryFile
39\. NtQueryInformationFile | 40\. NtQueryInformationProcess
41\. NtQueryInformationToken | 42\. NtQueryKey
43\. NtQueryObject | 44\. NtQuerySystemInformation
45\. NtQueryValueKey | 46\. NtQueryVirtualMemory
47\. NtQueryVolumeInformationFile | 48\. NtReadFile
49\. NtReadVirtualMemory | 50\. NtRequestWaitReplyPort
51\. NtResumeThread | 52\. NtSetContextThread
53\. NtSetInformationFile | 54\. NtSetInformationKey
55\. NtSetInformationProcess | 56\. NtSetInformationThread
57\. NtSetSecurityObject | 58\. NtSetValueKey
59\. NtTerminateProcess | 60\. NtTerminateThread
61\. NtUnmapViewOfSection | 62\. NtWaitForMultipleObjects
63\. NtWriteFile | 64\. NtWriteVirtualMemory
65\. NtYieldExecution | 66\. OpenSCManager
67\. OpenServiceW | 68\. SetWindowsHookEx

####  III-B2 Phase 2: Feature Selection

Feature selection is one of the potential steps toward building an efficient
ML model because it reduces the probability of overfitting by removing
redundant data, improves the model’s accuracy by eliminating irrelevant
features, and thereby reduces training time. At the beginning of the Feature
Selection phase, we have evenly divided (stratified train-test split) our
dataset into train data (75%) and test data (25%) to avoid data leakage. Next,
we have applied Recursive Feature Elimination with Cross-Validation (RFECV) to
our train data for selecting the most significant features for each ML
classifier. RFECV is a wrapper-style feature selection method that wraps a
given ML model and selects the optimal number of features for that model by
recursively eliminating $0-n$ features in each loop. Next, it selects the best
performing subset of features based on the accuracy or the score of cross-
validation. Besides, RFECV also removes the dependencies and collinearities
existing in the model. By using RFECV, we have selected 6 distinct subsets of
features for 6 ML classifiers. Table III presents the list of features that
have been selected by RFECV for each ML classifier. To represent the names of
the selected features, the serial numbers from Table II are used in Table III.
These features have been selected by setting ‘min_features_to_select’ as 34
(half of the features), cv=5, and ‘scoring’= ‘accuracy’ so that RFECV would
select at least half of the features based on the optimum accuracy over the
5-fold cross-validation.

Table III: List of features selected by RFECV for each ML classifier Classifier | Selected Features | Total
---|---|---
LR |  1,2,4-8,10,11,14-16,18-22,24,28,29,31-35, 37-39,42,43,45,46,48,53,54,56,59,61,63,64,67 | 41
SGD |  1,2,4-7,10,11,14-16,19,20,23-39,43, 45-48,51,53,54,56,59,61,63,64,67 | 44
KNN |  4-8,12-16,19,23-25,27-36,39-41, 43-48,50,51,53,54,56,59,61,63,64 | 42
NB |  1-4,6,10-12,14,15,17,19-22,24,28,30-34,36-38, 41-44,48,49,51,53,54,56,58,59,61-65,67,68 | 44
RF |  2,3,5,6,8-11,13,14,16,18-20,24,27-36,38-48, 50,51,54-56,59,62,63 | 44
SVM |  1,2,4-7,9,10,13-16,19-20,23-25,27-40,43, 45-48,50,51,53,54,56,59,61,63,64 | 44

###  III-C Classification

We have employed Supervised Machine Learning algorithms to detect and classify
15 ransomware families into corresponding categories. Supervised Machine
Learning algorithms are trained on the labeled dataset to make a decision in
response to the unseen test dataset. As our work focuses on detecting and
classifying 15 ransomware families, we have employed the following 6 state-of-
the-art Supervised Machine Learning algorithms that are widely used for both
binary and multi-class classification problems as per requirement:

Logistic Regression (LR): a type of statistical analysis that predicts the
probability of a dependent variable from a set of independent variables using
their linear combination.

Stochastic Gradient Descent (SGD): is an optimization algorithm to find the
model parameters by updating them for each training data so that the best fit
is reached between predicted and actual outputs.

K-Nearest Neighbor (KNN): estimates the likelihood of a new data point being a
member of a specific group by measuring the distance between neighboring data
points and the new data point.

Naïve Bayes (NB): is based on Bayes’ theorem and predicts the probability of
an instance belonging to a particular class.

Random Forest (RF): constructs multiple decision trees and finally determines
the class selected by the maximum number of trees.

Support Vector Machine (SVM): takes one or more data points from different
classes as inputs and generates hyperplanes as outputs that best distinguish
the classes.

Since we aim to do multi-class classification and some classifiers are only
designed for binary classification problems (i.e., Logistic Regression,
Support Vector Machine, etc.), we cannot directly apply these classifiers for
multi-class classification problems. Therefore, heuristic methods can be
applied to divide a multi-class classification problem into several binary
classification problems. There are two types of heuristic methods, such as
One-vs-Rest (OvR) which splits the dataset into one class against all other
classes each time, and One-vs-One (OvO) which splits the dataset into one
class against every other class each time. We have applied the OvR method for
our experiment to reduce the time and computational complexities. All these
classifiers are built along with ‘RandomSearchCV’- a hyperparameter
optimization technique, to find the best combination of hyperparameters for
maximizing the performance of the models’ output in a reasonable time. Instead
of exhaustively searching for the optimal values of the hyperparameters
through a manually determined set of values (i.e., Grid Search),
RandomSearchCV randomly searches the grid space and selects the best
combination of hyperparameter values based on the accuracy or the score of
cross-validation. Since we have used RFECV for feature selection and
RandomSearchCV for hyperparameter optimization, the nested cross-validation
technique has been applied to build our model.

##  IV Experimental Results

We have evaluated our models in terms of Precision, Recall, F1-score,
Accuracy, and Area Under the Receiver Operator Characteristic (AUROC). These
performance metrics are measured as follows:

$Precision=\frac{TP}{TP+FP}$ (1)

$Recall=\frac{TP}{TP+FN}$ (2)

$F1{\text{-}}score=\frac{2\times Precision\times Recall}{Precision+Recall}$
(3)

$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}\times 100$ (4)

where, TP = True Positives, FP = False Positives, TN = True Negative, FN =
False Negatives.

Table IV describes the experimental results showing that LR outperforms other
classifiers securing 99.15% overall accuracy. Although KNN achieves 85.65%
accuracy, the processing time of this classifier is less than other
classifiers as presented in Table V. To examine the degree of separability of
our best-performed machine learning model, we have constructed Area Under the
Receiver Operator Characteristic (AUROC) curve of the LR classifier as shown
in Fig. 4. We construct the AUROC curve by plotting the True Positive Rate or
Sensitivity against the False Positive Rate or (1-Specificity). For threshold
= 0.5, the average AUROC of LR is 99.55. This is an ideal situation for LR
which indicates that the model has a 99.55% probability of predicting the
correct class against the rest. Therefore, we can ascertain that LR is our
best-performing classifier.

Table IV: Performance ($P_{avg}$) comparison of LR, SGD, KNN, NB, RF, and SVM Classifiers $P_{avg}$ | LR | SGD | KNN | NB | RF | SVM
---|---|---|---|---|---|---
Accuracy | 99.15 | 85.26 | 85.65 | 94.96 | 82.38 | 94.23
Precision | 99.22 | 95.94 | 93.47 | 98.04 | 99.86 | 99.13
Recall | 99.15 | 87.62 | 85.65 | 94.96 | 82.38 | 94.23
F1-score | 99.15 | 91.14 | 88.04 | 96.28 | 89.31 | 96.37
AUROC | 99.55 | 93.67 | 92.55 | 97.41 | 91.19 | 97.09
Table V: Classifier’s processing time comparison Classifier | Processing Time (in seconds)
---|---
LR | 52.44
SGD | 53.65
KNN | 51.34
NB | 53.77
RF | 54.19

Figure 4: AUROC curve of Logistic Regression

##  V Model Explainability

Although some predictive models may not require explanation because of their
usage in a low-risk real-world environment, well-studied problem statements,
and the availability of sufficient practical solutions, some models that deal
with the real-world high-risk environment (i.e., ransomware
detection/classification) need explanation. When it comes to the performance
evaluation of any model, knowing both ‘What’ and ‘Why’ the models have taken
these decisions is equally important. These insights lead us to dive deeper
into the problems, explore the most important features that affect the model’s
prediction both locally and globally and thereby build the trustworthiness of
the model’s predictions. We present the explainability of our model’s
prediction with the help of ‘SHapley Additive exPlanations’ or SHAP values
[27]. SHAP is based on the coalition game theory that measures each feature’s
individual contribution to the final output while conserving the sum of
contributions being the same as the final result. Contrary to the other
explanation techniques, such as permutation feature importance, LIME, etc.,
SHAP provides both global and local explainability. Besides, we can use SHAP
values to explain a wide variety of models, such as DeepExplainer to explain
Deep Neural Networks (i.e., Multi-Layer Perceptron, Convolutional Neural
Network, etc.), TreeExplainer to explain tree-based models (i.e., Random
Forest, XGBoost, etc.), and KernelExplainer to explain any model, etc.

To illustrate both local and global explainability, we have used
KernerExplainer to explain our Logistic Regression model.

###  V-A Local Interpretation

For our classification model, the SHAP value is regarded as a 2-D array where
the columns represent the features used in the model and the rows represent
individual predictions predicted by the model. So, each SHAP value in this
array indicates a specific feature’s contribution to the output of that row’s
prediction as shown in Fig. 5. Here, a positive SHAP value specifies that a
feature is positively pushing the base value or expected value to the model
output. On the other hand, the negative SHAP value specifies that the feature
is negatively pushing the base value to the model output. The base value or
the mean model output is computed based on the train data. By using the Force
plot, we can visualize this explanation for an individual prediction as shown
in Fig. 6. For one instance, the features with higher SHAP values (red) are
positively pushing the base value to the model output, and the features with
lower SHAP values (blue) are negatively pushing the base value to the model
output.

Figure 5: Array of SHAP values

Figure 6: Force plot for single instance

###  V-B Global Interpretation

Passing the array of SHAP values to a ‘summary_plot’ function creates a
feature importance plot as shown in Fig. 7. It illustrates 20 highly
contributing features to the model output. Here, the x-axis denotes the mean
of the absolute SHAP value for each feature which indicates the total
contribution of the feature to the model and the y-axis denotes the features
used for the classification. The features are organized in descending order
from top to bottom by how strongly they influence the model’s decision. For
our LR model, the frequency of NtAlpcSendWaitReceivePort is the most important
feature with a higher mean absolute SHAP value followed by NtFsControlFile,
NtWriteVirtualMemory, etc.

Figure 7: Summary plot showing feature importance

In addition to the summary plot, we have created a ‘beeswarm’ plot as
presented in Fig. 8 which shows the distribution of the features’ impact on
the model output. The x-axis denotes the SHAP value and the y-axis denotes
each feature used for the classification. Every point on this plot stands for
one SHAP value for one feature and prediction. For our LR model, the lower
values (blue) of NtAlpcSendWaitReceivePort have a negative impact on the model
output, while the lower values of NtFsControlFile have a positive impact on
the model output.

Figure 8: Beeswarm plot showing directionality impact

After obtaining these insights from the explanation of our LR model, we
conducted another distinct experiment to reevaluate the performance of our LR
model including only 20 highly contributing features as shown in Fig. 7 and
Fig. 8. As illustrated in Table VI, our model could distinguish the ransomware
families with 95.48% overall accuracy based on those features. This infers
that the rest 21 features selected by RFECV for the LR classifier are also
contributing to the model achieving its best 99.15% overall accuracy.

Table VI: Performance evaluation of LR model including 20 highly contributing features Performance | LR (including only the top 20 features)
---|---
Accuracy | 95.48
Precision | 96.06
Recall | 95.48
F1-score | 95.45
AUROC | 97.59

To evaluate our model’s performance concerning the existing classification
methods, we have compared the accuracy of the previous API call-based
ransomware detection and classification approaches with the proposed method.
Both API call frequency and sequence-based detection methods were surveyed
from the literature. Table VII illustrates the comparison of the proposed
model with those previous works. It is evident that the proposed method
provides higher accuracy while having the advantage of multi-class
classification with model explainability. In addition, we have developed a
Web-Crawler to automate collecting 15 different ransomware families unlike [7]
where only 7 ransomware families were utilized for the experiment. Besides, we
have conducted feature engineering and utilized hyperparameter optimization
techniques to optimize each model’s performance contrary to [6]-[11].
Therefore, our approach is unique and efficient in multi-class classification.

Table VII: Performance comparison with existing ransomware classification models Reference |  Feature |  Classifier |  Multiclass |  Accuracy (%)
---|---|---|---|---
Maniath et al. [6] |  API call |  Long Short-Term Memory |  No |  96.67
VinayaKumar et al. [7] |  API call |  Multilayer Perceptron |  Yes |  98.00
Z.-G. Chen et al. [8] |  API call |  Logistic Regression |  No |  98.20
Takeuchi et al. [9] |  API call |  Support Vector Machine |  No |  97.48
Bae et al. [10] |  API call |  Random Forest |  Yes |  98.65
Hwang et al. [11] |  API call |  Markov Chain, Random Forest |  No |  97.30
Proposed Method |  API call |  Logistic Regression |  Yes |  99.15

##  VI Conclusion

Ransomware has become the most popular weapon to cybercriminals recently for
its easy way to earn money through untraceable cryptocurrency. In this
research work, we have developed a Web-Crawler to automate searching and
downloading both Crypto and Locker types of ransomware. By executing the
binaries of ransomware and benign software using an advanced sandbox tool, we
extracted API call frequencies, constructed our dataset, and selected features
through a two-phase feature engineering process. Next, we have applied six
state-of-the-art Supervised Machine Learning models to detect and classify
ransomware families. Our approach includes utilizing Recursive Feature
Elimination with Cross-Validation (RFECV) for selecting the significant
features and RandomSearchCV for selecting the optimum hyperparameter values
for each ML classifier. Thereby we attempt to optimize each model’s
performance. The results reveal that Logistic Regression can efficiently
classify ransomware into their corresponding families securing 99.15% overall
accuracy. Finally, instead of relying on the ‘Black box’ characteristic of the
Machine Learning models, we also present the explainability of our best-
performing model using SHAP values to ascertain the transparency and
trustworthiness of the model’s prediction. In the future, we aim to enlarge
our dataset by including more samples and potential distinguishing features.
In addition, we plan to apply state-of-the-art deep neural networks with
explainability for ransomware classification.

## References

  * [1] M. M. Hasan and M. M. Rahman, ”RansHunt: A support vector machines based ransomware analysis framework with integrated feature set,” 2017 20th International Conference of Computer and Information Technology (ICCIT), 2017, pp. 1-7, doi: 10.1109/ICCITECHN.2017.8281835.
  * [2] Florian Zandt, “The Industries Most Affected by Ransomware” statista.com, https://www.statista.com/chart/26148/number-ofpublicized-ransomware-attacks-worldwide bysector/ (accessed on April 6, 2022).
  * [3] “Ransomware Awareness for Holidays and Weekends — CISA,” www.cisa.gov. https://www.cisa.gov/uscert/ncas/alerts/aa21-243a
  * [4] “FinCEN: 2021 ransomware activity outpaces 2020 in 6 months,” SearchSecurity. https://www.techtarget.com/searchsecurity/news/252508296/FinCEN-2021-ransomwareactivity-outpaces-2020-in-6-months (accessed May 10, 2022).
  * [5] “VIRUSTOTAL RANSOMWARE ACTIVITY REPORT,” storage.googleapis.com. https://storage.googleapis.com/vtpublic/vt-ransomware-report-2021.pdf (accessed Sep 18, 2022)
  * [6] S. Maniath, A. Ashok, P. Poornachandran, V. G. Sujadevi, P. Sankar A.U. and S. Jan, ”Deep learning LSTM based ransomware detection,” 2017 Recent Developments in Control, Automation & Power Engineering (RDCAPE), 2017, pp. 442-446, doi: 10.1109/RDCAPE.2017.8358312.
  * [7] R. Vinayakumar, K. P. Soman, K. K. Senthil Velan and S. Ganorkar, ”Evaluating shallow and deep networks for ransomware detection and classification,” 2017 International Conference on Advances in Computing, Communications and Informatics (ICACCI), 2017, pp. 259-265, doi: 10.1109/ICACCI.2017.8125850.
  * [8] Z.-G. Chen, H.-S. Kang, S.-N. Yin, S.-R. Kim, ‘Automatic Ransomware Detection and Analysis Based on Dynamic API Calls Flow Graph’, In Proceedings of the International Conference on Research in Adaptive and Convergent Systems, New York, NY, USA, 196–201.
  * [9] Y. Takeuchi, K. Sakai, S. Fukumoto, ‘Detecting Ransomware Using Support Vector Machines’, In Proceedings of the 47th International Conference on Parallel Processing Companion, Eugene, OR, USA, 2018.
  * [10] S. I. Bae, G. B. Lee, E. G. Im, ‘Ransomware detection using machine learning algorithms’, Concurrency and Computation: Practice and Experience, Hoboken, NJ, USA, 2019; Volume 32.
  * [11] J. Hwang, J. Kim, S. Lee, and K. Kim, “Two-stage ransomware detection using dynamic analysis and Machine Learning Techniques,” Wireless Personal Communications, vol. 112, no. 4, pp. 2597–2609, 2020.
  * [12] J. Baldwin and A. Dehghantanha, “Leveraging Support Vector Machine for Opcode Density Based Detection of Crypto-Ransomware,” Advances in Information Security, pp. 107–136, 2018, doi: 10.1007/978-3-319- 73951-9 6.
  * [13] J. Hwang, J. Kim, S. Lee, and K. Kim, “Two-stage ransomware detection using dynamic analysis and Machine Learning Techniques,” Wireless Personal Communications, vol. 112, no. 4, pp. 2597–2609, 2020.
  * [14] K. P. Subedi, D. R. Budhathoki and D. Dasgupta, ”Forensic Analysis of Ransomware Families Using Static and Dynamic Analysis,” 2018 IEEE Security and Privacy Workshops (SPW), 2018, pp. 180-185, doi: 10.1109/SPW.2018.00033.
  * [15] A. Ferrante, M. Malek, F. Martinelli, F. Mercaldo, and J. Milosevic, “Extinguishing Ransomware - A Hybrid Approach to Android Ransomware Detection,” Foundations and Practice of Security, pp. 242–258, 2018, doi: 10.1007/978-3-319-75650-9 16.
  * [16] K. A. Roundy and B. P. Miller, “Binary-code obfuscations in prevalent Packer Tools,” ACM Computing Surveys, vol. 46, no. 1, pp. 1–32, 2013.
  * [17] K. Coogan, S. K. Debray, T. Gregg and M. Townsend, “Automatic static unpacking of malware binaries,” in Proc. 16th Working Conf. Reverse Eng., Lille, France, Oct. 2009, pp. 13–16.
  * [18] VirusShare.com. Available: http://virusshare.com/. (accessed Sep 18 2022).
  * [19] R Mowri, “RMOWRI/GetRansomware,” GitHub. https://github.com/rmowri/GetRansomware. (accessed Sep 18 2022).
  * [20] tisf, “theZoo - A Live Malware Repository,” GitHub. http://github.com/ytisf/theZoo (accessed May 10, 2022).
  * [21] “Free Automated Malware Analysis Service - powered by Falcon Sandbox,” Hybridanalysis.com, 2019. https://www.hybrid-analysis.com/ (accessed Sep 18, 2022)
  * [22] B. A. Al-rimy, M. A. Maarof, and S. Z. Shaid, “Ransomware threat success factors, taxonomy, and countermeasures: A survey and research directions,” Computers Security, vol. 74, pp. 144–166, 2018.
  * [23] A. Al-Bakri, H. Hussein, ‘Static Analysis Based Behavioral API for Malware Detection using Markov Chain’, Computer Engineering and Intelligent Systems, Jan 2014.
  * [24] S. Al amro, A. Cau, “Behavioural API based Virus Analysis and Detection”, International Journal of Computer Science and Information Security, Vol. 10, No. 5, 2012.
  * [25] “Falcon Sandbox: Automated Malware Analysis Tool — CrowdStrike,” crowdstrike.com. https://www.crowdstrike.com/products/threatintelligence/falcon-sandbox-malware-analysis/ (accessed May 10, 2022).
  * [26] PayloadSecurity, “Payloadsecurity/vxapi: A generic interface and CLI for all endpoints of the Falcon Sandbox API,” GitHub. https://github.com/PayloadSecurity/VxAPI. (accessed on Sep 18, 2022).
  * [27] S. Lundberg and S.-I. Lee, “A unified approach to interpreting model predictions,” arXiv.org, May 2017. (accessed Sep 18,2022).
  * [28] “Latest entries - The Portable Freeware Collection,” www.portablefreeware.com. https://www.portablefreeware.com/ (accessed May 10, 2022).


# Hafnia HfO2 as a Proper Ferroelectric

Aldo Raeliarijaona<EMAIL_ADDRESS>R. E. Cohen
<EMAIL_ADDRESS>Extreme Materials Initiative, Earth and Planets
Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW,
Washington, DC 20015, USA

###### Abstract

We clarify the nature of hafnia as a proper ferroelectric and show that there
is a shallow double well involving a single soft polar mode as in well-known
classic ferroelectrics. Using symmetry analysis, density-functional theory
(DFT) structural optimizations with and without epitaxial strain, and density
functional perturbation theory (DFPT), we examine several important possible
hafnia structures derived ultimately from the cubic fluorite structure,
including baddeleyite (P$2_{1}/c$), tetragonal antiferroelectric P42nmc, Pbca
(OI and brookite), and ferroelectric rhombohedral ($R3m$), Pmn21 and Pca21
structures. The latter is considered to be the most likely ferroelectric phase
seen experimentally, and has an antiferroelectric parent with space group
Pbcn, with a single unstable polar mode and a shallow double well with a well
depth of 24 meV/atom. Strain is not required for switching or other
ferroelectric properties, nor is coupling of the soft-mode with any other
modes within the ferroelectric Pca21 phase.

## Introduction

Hafnia has great potential as a technologically important material due to
experimentally observed ferroelectricity in thin films[1, 2] and its
compatibility with silicon for use in integrated electronics[3, 4]. The nature
of the observed ferroelectricity is complicated by a complex set of potential
phases and distortions from the parent cubic fluorite structure, and by the
dependence of its properties and phase relations due to doping, strain, and
applied electric fields that give rise to the so-called wake-up effect. Here
we clarify that ferroelectric hafnia can be understood as a simple, proper,
ferroelectric, from which perturbations from doping, strain, and applied
fields can be explored straightforwardly.

Complications in understanding hafnia stem from the fact that there are no
known stable ground state polar phases, even varying pressure and temperature
[5]. At ambient pressure, the phase diagram shows the fluorite cubic phase
($Fm\bar{3}m$) above 2600K, tetragonal phase (P42nmc) between 1800K and 2600K,
and monoclinic baddeleyite (P21/c) below 1800 K[5]. Increasing pressure at
room temperature gives non-polar orthorhombic Pbca for 5 GPa$\leq p\leq$ 15
GPa, and Pnma for $p>$ 15 GPa.

The situation is clear if we compare with well known classic perovskite
ferroelectrics which are now well understood [6, 7]. In the perovskites, the
parent paraelectric phase has simple cubic perovskite structure with five
atoms per unit cell, but this does not mean that the atoms all sit on their
ideal sites. Rather the atoms hop among symmetric off-centered sites, and only
the average structure has the simple cubic structure. As the temperature is
lowered at ambient pressure, three phase transitions are observed in BaTiO3,
for example, cubic to tetragonal ($P4mm$) around 400K, tetragonal to
orthorhombic ($Amm2$) around 280K, and orthorhombic to rhombohedral ($R3m$)
around 200K[8, 9]. Although the transitions are understood to be distortions
of the parent average structure due to a soft phonon that is polar, and only
the ground state rhombohedral phase is an ordered structure with atoms
vibrating around their ideal sites in the $R3m$ structure. In contrast, the
classic ferroelectric LiNbO3, has ten atoms per unit cell, from a zone
boundary instability in the parent cubic perovskite. The resulting
paraelectric R$\bar{3}$c phase has a single polar soft mode giving a proper
double well to a polar, ferroelectric $R3c$ ground state. Just like in the
simple perovskite ferroelectrics, the atoms do not sit on their ideal sites in
the ferroelectric phase, but hop [10, 11]. Unlike BaTiO3 etc., in LiNbO3,
cubic perovskite never forms below the melting point. As we show below, hafnia
is akin to LiNbO3 in that the ferroelectric phase is a proper ferroelectric
with a single polar mode, but unlike LiNbO3, at high temperatures the parent
tetragonal and cubic phases do form below the melting point, partly due to its
very high melting temperature of 3031K. Indeed, similar to hafnia, originally
Megaw suggested that LiNbO3 switched through the cubic perovskite as its
paraelectric parent[12], but Barker and Loudon showed that the parent is
R$\bar{3}c$[13], which was confirmed in later work [10, 14].

At high temperatures, cubic fluorite hafnia has a soft zone-boundary mode,
$X^{-}_{2}$, leading to tetragonal [15] at 2600K. However, tetragonal hafnia
is dynamically stable with a first-order phase transition to monoclinic
baddeleyite at 1800K. The absence of a polar, stable ferroelectric phase made
the discovery of ferroelectricity in thin film hafnia a great surprise[1].

Previous experiments and theory show that ferroelectric hafnia is most likely
orthorhombic with space group Pca21[16, 17]. The energy barriers between
ferroelectric Pca21, the non-polar ground state baddeleyite P21/c,
orthorhombic antiferroelectric $Pbcn$, and tetragonal antiferroelectric P42nmc
are all large, larger than kT at room temperature for pure hafnia[18].
However, once formed, ferroelectric Pca21 is metastable and switchable. There
have a been a number of proposals for formation of the ferroelectric phase,
including doping and vacancies[19], strain[18], kinetics[20], and electric
field[21, 22]. Other polar phases such as Pmn21 [23, 24] and rhombohedral
$R3m$[25, 26] are also metastable albeit with energies larger than that of
Pca21. The former has been only proposed theoretically[23, 24], whereas the
latter has been experimentally observed in (111)-oriented Hf0.5Zr0.5O2[25,
26]. Ref. 27 presented the tree of subgroups from fluorite in hafnia (and
zirconia), mostly consistent with our work, though we concentrate on the
important phases leading to ferroelectricity, and neglect any phases with
partially occupied sites, which seem improbable. We concentrate on the
symmetry modes (soft and hard modes) over group/subgroup relations. We discuss
differences from Nentwich’s analysis below in the discussion section.

Cheema et al. found that very thin films of only a couple unit cells thick
have enhanced ferroelectric response and lattice parameters that differ
significantly from thicker films or bulk[28], but we do not address such thin
films here using bulk computations. Thicker films of 3-10 nm have lattice
parameters consistent with our computations even at zero stress, consistent
with previous bulk computations[29, 30] and what we report below. Ref. 31
considers the large strains in bulk computations and find that they help
promote stability of the ferroelectric Pca21 phase. Here we show that such
large strains are not necessary for ferroelectric behavior, consistent with
experimental observations.

Ferroelectricity in hafnia has been called improper [15, 32, 31] and for
ferroelectric Pca21 several parent structures have been proposed, such as
fluorite cubic[15], $Ccce$[32], $Pbcm$[33], and tetragonal[23, 31]. For
ferroelectric Pmn21, some correctly considered tetragonal P42nmc as the parent
structure[23, 24], which would be a proper soft-mode transition. Since polar
$R3m$ has a higher symmetry than most phases except the cubic, it must have
the fluorite cubic as a parent structure, and it would also be a proper
transition from cubic.

Ref. 15 concluded that ferroelectricity in hafnia is improper, and is due to
antipolar-polar coupling and nonlinear interactions with the primary phonon
instability in cubic $X^{-}_{2}$. In a similar but different way, Zhou et al.
showed that tensile uniaxial strain and antipolar-polar coupling leads to the
ferroelectric Pca21 and claimed that hafnia is an improper ferroelectric[31].
Using a different parent structure $Ccce$ (equivalently $Bbab$) Ref. 32
claimed that three modes, $\Gamma_{3-}$, $Y_{3+}$, and $Y_{2+}$, are coupled
to generate the Pca21 structure. In contrast, Ref. 33, using $Pbcm$ as the
parent structure, stated correctly that ferroelectricity in hafnia is proper,
but $Pbcm$ is a higher energy parent, and gives deeper double wells than we
show below.

In this paper we clarify symmetry relationships among the different phases of
hafnia and the phonon modes that connect them. We find that antiferroelectric
$Pbcn$ has a zone center unstable mode $\Gamma^{-}_{2}$, whose distortion
leads to Pca21. We further find that the double-well potential resulting from
this $\Gamma^{-}_{2}$ distortion is shallow, making $Pbcn$ key in the
polarization switching of Pca21 hafnia. Having a correct model for hafnia’s
ferroelectric modes and double wells is important for learning to control its
properties by doping, strain, and field, as well for proper understanding of
its material physics. In order to optimize its switching behavior and design
domain engineering and best epitaxial growth, a correct understanding of the
ferroelectric instabilities are crucial.

## Methods

We performed symmetry analyses starting with the fluorite structure and the
many potential structures from subgroups of the fluorite structure using
isodistort[34, 35] and the Bilbao Crystallographic Server[36, 37, 38]. We
performed first-principles calculations using Quantum Espresso (QE) [39, 40,
41] using Garrity-Bennett-Rabe-Vanderbilt (GBRV)[42] pseudopotentials from the
standard solid-state pseudopotential (SSSP)[43] with PBEsol [44] exchange-
correlation. The plane-wave expansion cutoff energy E${}_{\textbf{cutoff}}$
was 544 eV, and the Brillouin zone was sampled using an 6$\times$6$\times$6
Monkhorst-Pack grid [45]. Convergence tests were performed for all parameters.
Phonon frequencies were obtained using Density Functional Perturbation Theory
(DFPT) [46] implemented in the QE/PH package[39, 40, 41] and the phonon
dispersion curves along the $\Gamma$–X–M–$\Gamma$ high symmetry line of the
Brillouin zone (Fig.3) are an interpolation of the interatomic force constants
(IFC) computed using 4$\times$4$\times$4 q-vector grid. The computed energy
differences, the relaxed lattice lattice parameters of the various polymorphs,
and their Raman spectra compare well with previous studies[47, 48].

We parameterized the studied distortions using the displacive order parameter
$Q$ as the mode amplitude:

$Q=\sqrt{\sum_{i}{\textbf{u}_{i}}^{2}}$ (1)

where ${\textbf{u}}_{i}$ is the displacement of atom $i$ in Angstroms (Å).

## Results

### Group-subgroup relations of different hafnia polymorphs

First we consider instabilities of the high temperature tetragonal P42nmc
phase because it is also a precursor to antiferroelectric and ferroelectric
orthorhombic phases and the ground state baddeleyite. Using group theory we
investigate how baddeleyite, Pca21, Pmn21, and the $Pbcn$ can be obtained from
a transformation of the tetragonal phase.

Figure 1: Summary of the transitions to different hafnia polymorphs starting
from cubic fluorite Fm$\bar{3}$m(topmost) to tetragonal P42nmc, Pmn21 (top
right), Pbcn (middle), Pca21 (bottom right), and baddeleyite (bottom left).
The symmetry designations of the order parameter responsible for the
distortion are labelled near the arrow, and are all written in terms of the
parent phase (tail of the arrow), as it is done throughout this paper. The
blue symmetry labels indicate soft modes and the outlines depict the unit cell
relationships. Solid lines with arrows show soft-mode relationships, whereas
dashed lines must be first-order transitions. The dark blue arrows on the
atoms of the FE and antiferroelectric (AFE) structures show the direction of
the polarization. The angle $\beta$ of the monoclinic structure is 99.06∘.

We find that tetragonal P42nmc is related to antiferroelectric $Pbcn$ by a
displacive order parameter with $M1$ symmetry (Fig.1). $M1$ is a soft mode of
tetragonal P42nmc. Pca21 is related to tetragonal P42nmc through the coupling
of the $M1$ and $\Gamma^{-}_{5}$ modes, and baddeleyite results from the
coupled $M1$ and $\Gamma^{+}_{5}$ distortions. We use simplified symmetry
designations of the parent structure modes (tetragonal P42nmc in this case)
according to isodistort[34, 35]. We found no group-subgroup relation between
Pca21 and baddeleyite, so any transition from ground state baddeleyite to
ferroelectric Pca21 must be through nucleation and growth.

$Pbcn$ is a parent of Pca21 and baddeleyite P21/c. A distortion of $Pbcn$ with
$\Gamma^{-}_{2}$ or $\Gamma^{+}_{2}$ symmetry generates Pca21 and baddeleyite
P21/c respectively (Fig.1). The $\Gamma^{-}_{2}$ mode is soft for relaxed
$Pbcn$ and for epitaxially strained $Pbcn$, a${}_{\text{epi}}$$\leq$ 5.2Å,
where a${}_{\text{epi}}$ is the lattice constant of a square substrate. To lay
an orthorhombic cell on a square substrate the diagonals of the film and the
substrate were matched and $a_{\text{epi}}=\sqrt{\frac{a^{2}+b^{2}}{2}}$ with
$a$ and $b$ the in-plane lattice parameters of the film. The layout of $Pbcn$
is such that the [110] and [1$\bar{1}$0] of $Pbcn$ are along the xy and
x$\bar{y}$ of the substrate respectively. Similarly, tetragonal is laid out
such that the [110] and [1$\bar{1}$0] are along the x and y of the substrate.

Although not observed experimentally, Pmn21 is the subject of many
computational studies[24, 23, 49]. Tetragonal P42nmc is distorted into the
Pmn21 structure by its $\Gamma^{-}_{5,X}$ mode. The $\Gamma^{-}_{5,XY}$ mode
leads from tetragonal to symmetry $Aba2$, but this structure is a saddle point
and relaxes to centrosymmetric Bbab. Bbab is stable and is 29 meV/atom above
baddeleyite in energy, similar to Pbcn but above Pca21 and $Pbca$. There are
two different Pbca structures, and they are often confused. They were properly
distinguished by Kersch and Falkowski, who called them nonpolar uP61 and
antipolar aP61, which is the Brookite structure [50]. The designation “OI”, as
denoted in Ref. [30], is sometimes applied to brookite and sometimes to the
nonpolar structure. They have the same Wyckoff positions, but different atomic
coordinates (Table 1 and Fig. 2). The o-phase in Ref.22, 50 referred to the
nonpolar structure, and in Ref.51, 50 to the antipolar brookite. Although
brookite might be considered antiferroelectric (Fig.2), there is no phonon
mode that connects it to Pca21 nor is there a group-subgroup relationship
between them.

There are lots of local minima, which complicates the relationships among the
metastable phases. Any of these metastable and saddle point structures might
be stabilized entropically with temperature by anharmonic thermal vibrations
or hopping.

Figure 2: Summary of the transitions to different hafnia polymorphs starting
from cubic fluorite Fm$\bar{3}$m(center) to rhombohedral $R3m$ (left), Pca21
(bottom right), orthorhombic Pbca (bottom left), and baddeleyite (right). The
symmetry designations of the order parameter responsible for the distortion
are labelled near the arrow, and are all written in terms of the parent phase
(tail of the arrow). The outlines depict the unit cell relationships. The
polarization for the rhombohedral phase is out of the page as indicated, along
the Cartesian z-direction. Possible soft-mode transitions are indicated by
solid arrows, and possible first-order transition involving hard modes are all
dashed. The dashed-dotted lines to the two Pbca structures indicate either a
first order transition from cubic to $Pbca$ or a complicated combination of 6
optical modes (cubic modes $X^{+}_{5}$, $X^{-}_{5}$, $X^{1}_{3}$, $W_{3}$,
$W_{4}$, and $\Delta_{5}(0,\frac{1}{2},0)$) and strain. Both Pbca structures
(brookite antipolar and the nonpolar structure) have the same cubic modes
involved in the coupling but their difference is in their linear combinations.

Tetragonal P42nmc results from a distortion of cubic fluorite Fm$\bar{3}$m
according to the $X^{-}_{2}$ phonon mode (Fig. 1). Complicated coupling among
six modes plus strain leads to Pbca from cubic, which would be first-order or
have intermediate structures (as in Ref. 27). Coupled distortions (mainly
cubic modes $X^{+}_{5}$ and $X^{+}_{3}$) lead to baddeleyite. Rhombohedral
$R3m$, proposed by [52], is a distortion of the high-symmetry fluorite
following the $\Gamma^{-}_{4,XYZ}$ mode (Fig.2). $\Gamma^{-}_{4}$ is three-
dimensional (Tu), and the linear combination of the degenerate modes:

$\Gamma^{-}_{4,XYZ}=a\Gamma^{-}_{4,X}+b\Gamma^{-}_{4,Y}+c\Gamma^{-}_{4,Z}$ (2)

is denoted by (a,b,c), where $a$,$b$, and $c$ are coefficients of the linear
combination. For (a,a,a), the distortion leads to the rhombohedral $R3m$
phase; this is a hard mode of fluorite which we find to have $\omega$= 235
cm-1. There is no soft-mode that leads to the proposed [52] rhombohedral
ferroelectric. Coupling $\Gamma^{-}_{4,X}$ to $X^{+}_{5}$, on the other hand,
leads to Pca21.

We stress that in cubic fluorite Fm$\bar{3}$m, the only soft mode is the zone
boundary $X^{-}_{2}$ that leads to tetragonal P42nmc. The other distortions,
which lead to rhombohedral, Pca21, $Pbca$, and baddeleyite, are all hard
modes.

a b c Pbca (oI) 5.13 Å 5.25 Å 10.07 Å Hf 0.042 0.157 0.862 O 0.169 0.666 0.466
O 0.253 0.589 0.224 Pbca[22] 5.15 Å 5.29 Å 10.14 Å Hf 0.041 0.343 0.138 O
0.325 0.160 0.032 O 0.248 0.589 0.275 Brookite 5.03 Å 9.95 Å 5.20 Å Hf 0.246
0.116 0.536 O 0.002 0.023 0.239 O 0.128 0.710 0.875 Pbca[51] 5.06 Å 10.02 Å
5.23 Å Hf 0.244 0.385 0.034 O 0.372 0.290 0.372 O 0.000 0.023 0.243

Table 1: Lattice parameters and atomic positions of different Pbca structures.
a, b, and c are the crystallographic direction of the orthorhombic Bravais
Lattice. The atomic positions are given relative to the Bravais lattice
vectors. All the atoms sit on Wyckoff position 8c.

### Phonon dispersion

The tetragonal phase is dynamically stable at zero strain and under
compressive strain ($0\leq\eta\leq$ -3.75%) (Fig.3(a) and (b)). However, under
tensile strain ($\eta\geq 1.5\%$ or a${}_{\text{epi}}$$\geq$ 5.11Å) tetragonal
hafnia becomes unstable. Indeed for $\eta=3.75\%$, corresponding to
a${}_{\text{epi}}$= 5.22 Å (Fig.3(c)), we find soft modes at the Brillouin
zone center ($\Gamma^{-}_{5}$), the zone boundary points M (M1), X, and A. The
lowest branch is soft along the whole M–A–$\Gamma$ line, indicating the
possibility of incommensurate or disordered ferroelectric/antiferroelectric
ordering possibly explaining previous observations[15]. Along the A–$\Gamma$
portion of the phonon dispersion curve around $\omega=$ 400 cm-1, we observe
avoided crossing of the optical phonon branches. This involves optical
branches of the same symmetry, which tracks back to the $E_{u}$ at the
$\Gamma$ point. The frequency difference between these two branches decreases
as the tensile strain increases, i.e. the strain becomes more positive.

Figure 3: Phonon dispersion of tetragonal P42nmc at strain (a) $\eta=$ -3.75
%, (b) $\eta=$0%, and (c) $\eta=$3.75%. The dotted lines on panel (c) shows
the zero frequency level. The epitaxial strain $\eta$ is a biaxial in-plane
strain with $c$ allowed to relax, so that $a_{\eta}=(1+\eta)a_{0}$ where
$a_{0}$ is the lattice constant the stress-free tetragonal ($a_{0}$ =5.00Å).
Shown are the harmonic frequencies, i.e. the square root of the curvature of
the potential surface for atoms in their ideal positions.

### Soft phonon eigendisplacements

The symmetry mode for the double well in the ferroelectric Pca21 phase has
$\Gamma_{2}^{-}$ symmetry in the antiferroelectric parent $Pbcn$. It is an
unusual case, though, compared with what one is used to seeing for polar modes
in perovskite ferroelectrics. At high temperatures the atoms can hop between
the up and down positions in $Pbcn$, which is the tetragonal P42nmc M1 mode,
leading to the tetragonal structure. This situation holds for the relaxed
lattice and for epitaxial strains up to a${}_{\text{epi}}$= 5.19 Å.

Figure 4: The arrows show the polar, soft $\Gamma_{2}^{-}$ mode of $Pbcn$ that
leads to Pca21. The $Pbcn$ antipolar direction is top to bottom, and the blue
and green arrows indicate displacements of front and back oxygen atoms
respectively. The top oxygen atoms are moving out of the page and the bottom
ones are moving into the page as indicated by the symbols adjacent to the
structure. Pca21 can be thought of as polar in one direction ($z$), and
antipolar in the perpendicular $x$ and $y$ directions. The displacements of
the hafnium atoms are small in comparison to the oxygen atoms, so are not
shown.

The $\Gamma^{-}_{5}$ mode of tetragonal P42nmc becomes unstable at 1.5%
epitaxial strain (a${}_{\text{epi}}$= 5.11 Å). It is a two-dimensional mode,
with two symmetrically equivalent displacive order parameters, denoted
$Q_{\Gamma^{-}_{5,X}}$ and $Q_{\Gamma^{-}_{5,Y}}$, represented as (a,0) and
(0,a), respectively using isodistort[34, 35] notation. We find that the
eigendisplacements of the $\Gamma^{-}_{5}$ modes are related by screw symmetry
($4_{2}$), i.e. the combination of 90∘ rotation about the c-axis followed by a
translation by half the lattice constant along the c (Fig.5). The distortion
of either $Q_{\Gamma^{-}_{5,X}}$, or $Q_{\Gamma^{-}_{5,Y}}$ by itself leads to
the Pmn21 ferroelectric phase, and its linear combination:

$Q_{\Gamma^{-}_{5,XY}}=a~{}Q_{\Gamma^{-}_{5,X}}+b~{}Q_{\Gamma^{-}_{5,Y}},$ (3)

leads to the structure with $Aba2$ symmetry when the coefficients $a=b$.

Figure 5: The doubly degenerate $\Gamma_{5}^{-}$ soft modes of tetragonal
P42nmc, denoted by $\Gamma^{-}_{5,X}$ (panels a, b , and c); and
$\Gamma^{-}_{5,Y}$ (panels d, e, and f). The structures were projected on the
ab plane (panels a and d), ac plane (panels b and e), and bc plane (panels c
and f). The blue arrows on the atoms indicate the directions of their
displacements.

$Pbcn$, the centrosymmetric parent of ferroelectric Pca21, is generated from
tetragonal P42nmc from the $M1$ distortion (Fig.1). The order parameter
associated with the doubly degenerate soft $M1$ mode is denoted by $Q_{M1}$
and represented as (a,0), and their eigendisplacements (Fig.6) are related by
$4_{2}$ screw symmetry.

Figure 6: The doubly degenerate soft $M1$ modes of tetragonal P42nmc, denoted
by $M1_{1}$ (panels a, b , and c); and $M1_{2}$ (panels d, e, and f). The
structures were projected on the ab plane (panels a and d), ac plane (panels b
and e), and bc plane (panels c and f). The blue arrows on the atoms indicate
the directions of their displacements.

A fourth order expansion of the energy per atom in terms of the $M1$
distortion gives:

$\mathcal{E}(Q_{M1})=a~{}Q^{2}_{M1}+b~{}Q^{4}_{M1}+...,,$ (4)

with a = -7.28 meV/(Å2 atom), and b = 1.12 meV/(Å4 atom) obtained from a
polynomial fitting of the M1 curve (orange star) (Fig.9). The lattice strain
is fixed for the distortion, with lattice parameters relaxed at $Q=0$.

With $Pbcn$ as the parent structure, there is a single polar, soft mode
$\Gamma^{-}_{2}$ that leads to Pca21 (Fig. 1) ($\omega(\Gamma^{-}_{2})$ =
128$i$). The expression of the energy to the fourth order of $\Gamma^{-}_{2}$
is:

$\mathcal{E}(Q_{\Gamma^{-}_{2}})=a^{{}^{\prime}}~{}{Q^{2}_{\Gamma^{-}_{2}}}+b^{{}^{\prime}}~{}{Q^{4}_{\Gamma^{-}_{2}}}+...,$
(5)

where $Q_{\Gamma^{-}_{2}}$ is the order parameter associated to mode
$\Gamma^{-}_{2}$. A polynomial fit to the first-principles data in Fig. 7
yields a${}^{{}^{\prime}}$ = -36.55 meV/(Å2 atom),and b${}^{{}^{\prime}}$ =
14.77 meV/(Å4 atom).

In $Pbcn$ the oxygen displacements can be considered as up and down along the
antipolar axis (Fig.7) and the polar mode (Fig.4) brings one of these
displacements (or local polarizations) for half the oxygen atoms in the unit
cell to zero, relative to the surrounding hafnia atoms, while the displacement
for the other set of oxygens grows to double that in the parent $Pbcn$ phase.
When the polarization is reversed, the displaced oxygens move to the center of
their polyhedra, while the other set of oxygens gain a large polarization. So
the parent phase has “up, down” polarization, where one side of the double
well has “zero, double down” and the other side has “double up, zero”. The
parent structure can be considered as an ordered antiferroelectric.

Figure 7: Potential well for the $\Gamma^{-}_{2}$ symmetry mode leading from
$Pbcn$ to Pca21for relaxed lattice and different biaxial epitaxial strains.
The structures corresponding to the three extrema of the double well are drawn
above them, and the polarization switching is shown by the polarization
configurations: $\textbf{0}\uparrow$ (right well), $\downarrow\uparrow$ (top
of double well), and $\downarrow\textbf{0}$ (left well).

For completeness, we also discuss the stabilization of ferroelectric Pmn21.
For an extensional epitaxial strain $\eta=1.5\%$ (a${}_{\text{epi}}$= 5.11Å),
Pmn21 is more stable than P42nmc along $Q_{\Gamma^{-}_{5,X}}$. The fourth
order expansion of the energy of this P42nmc to Pmn21 phase transition,
without any coupling is:

$\mathcal{E}=a~{}{Q^{2}_{\Gamma^{-}_{5,X}}}+b~{}{Q^{4}_{\Gamma^{-}_{5,X}}}+...,$
(6)

a polynomial fit to the $\Gamma^{-}_{5,X}$ curve in Fig. 8 gives a=-2.46
meV/(Å2 atom), b = 0.43 meV/(Å4 atom). The Pmn21 structure is further
stabilized by coupling the $\Gamma^{-}_{5,X}$ distortion to the
$\Gamma^{+}_{2}$ and shear strain (Fig.9). The linear combination of
$Q_{\Gamma^{-}_{5,X}}$ with its degenerate mode $Q_{\Gamma^{-}_{5,Y}}$ gives a
structure with slightly lower energy with symmetry $Aba2$ as determined by
findsym[53, 54]. $Aba2$ is a saddle point, and relaxes to centrosymmetric
$Bbab$.

Figure 8: Energy as a function of $\Gamma^{-}_{5,X}$ and $\Gamma^{-}_{5,XY}$
distortions for relaxed tetragonal P42nmc(blue curve), and for
a${}_{\text{epi}}$= 5.21 Å, and at a${}_{\text{epi}}$= 5.285 Å. Figure 9:
Energy as a function of distortion for Pmn21. The curves are labelled
according to the irreducible representation, and $\eta$ indicates the
tetragonal lattice parameters are relaxed. The symmetry of the parent and
distorted structure are labelled near the top and bottom of the double-well,
respectively. The red, green and blue curves were obtained for
a${}_{\text{epi}}$= 5.21 Å, whereas the orange curve was obtained for
a${}_{\text{epi}}$= 5.25 Å.

###  .1 Discussion

We have shown that there are single polar soft modes that distort the
tetragonal and the $Pbcn$ parent phases to the polar Pmn21 and Pca21 phases,
respectively, so that hafnia is a proper ferroelectric. The centrosymmetric
parent for the experimentally observed Pca21, is $Pbcn$, which is dynamically
unstable.

We now consider possible ways to access $Pbcn$ and thus Pca21, since they are
not known to be stable at any conditions of pressure and temperature. Hafnia-
based samples are usually annealed at high temperature ($>$773K)[55].
Although, the temperature–pressure phase diagram of pure hafnia shows that
tetragonal is the stable phase at ambient pressure within the temperature
range of 1800K to 2600K[5], the transition temperature usually decreases with
doping concentration[19, 56] so that the high temperatures ($>$773K) used in
experiments[55] in doped hafnia are possibly large enough to stabilize the
tetragonal phase. The strain imposed by the substrate (e.g. yttria-stabilized
zirconia) destabilizes the tetragonal, which distorts into $Pbcn$.

Although Ref. 31 indicated $Pbcn$ as a key intermediate phase leading to
Pca21, the parent considered was the tetragonal phase which led to their
conclusion about the improper nature of ferroelectricity in hafnia.

In considering $Pbcm$ as the reference structure instead of cubic, Ref. 33
showed the existence of a single soft polar mode that gives Pca21, also
consistent with proper ferroelectricity in hafnia. $Pbcm$, however, is not a
subgroup of the tetragonal phase. The double well depth from $Pbcn$ to Pca21
is shallower than $Pbcm$, ( 72 meV/f.u vs 117 meV/f.u.[33]). The correct
switching pathway is the lowest energy pathway.

The lowest energy parent structure we find for the ferroelectric Pca21 phase
is $Pbcn$, which is thus the low energy pathway for switching the polarization
direction in hafnia. There may not be a thermodynamically stable phase field
for $Pbcn$–that is not required. However, even if formed metastably, it would
probably consist of disordered off-centered cells with atoms hopping back and
forth between double wells, and have at least partly order-disorder character.

The Bärnighausen tree of the HfxZr1-xO3 solid solution[27] gives the group-
subgroup relations between different hafnia polymorphs. The phases and paths
shown there are not complete, since the parent of different structures and the
sequence of transitions depends on which mode is first frozen. The fact that
Bärnighausen trees are not unique may not have been generally recognized. For
example in the cubic to $Pbca$ transition, the freezing of the $W3$ or
$X^{-}_{3}$ distorts the cubic into $I4_{1}/acd$ or $P4/nmm$, which are
possible distortions not shown in Ref.[27]. On the other hand, the number of
subgroup possible is infinite as pointed out by Bärnighausen[57] so there is
always a choice of which to include. The goal here is to concentrate on low
energy structures and the most probable routes to them from the high
temperature cubic and tetragonal phases. We focus on distortions involving
optical phonons, and provide physical bases of phase transition through DFT
computations. Additionally, the soft phonon modes and the ensuing double well
potential of the $Pbcn$ indicates that there is no need for reconstructive
transitions as proposed by Ref. 27. We find the maximum ion displacements from
$Pbcn$ to baddeleyite or Pca21 are only 0.45 Å and 0.25 Å respectively, much
smaller than shown in Ref. 27. We also do not consider static structures which
partially occupied sites included in Ref. 27 to be likely. For example, Pbcm
included there as a parent of $Pbca$ and Pca21was refined with half-occupied
sites from in situ diamond anvil cell X-ray data for ZrO2 [58], but neutron
diffraction on quenched samples showed $Pbca$ [59, 60]. However, the existence
of transitions at finite temperatures that hop among partially occupied sites
cannot be excluded.

##  I Conclusions

The number of possible phases, local minima, and saddle points for the ABO2
fluorite aristotype are truly astounding. We have shown the symmetry relations
between the high temperature phases of hafnia to the relevant low symmetry
structures, baddeleyite, Pca21, Pmn21, AO $Pbca$, $R3m$ and the $Pbcn$.
Tetragonal P42nmc hafnia is dynamically unstable under tensile strain, and the
soft mode $M1$ leads to the antiferroelectric $Pbcn$ structure. $Pbcn$ is the
best centrosymmetric parent for ferroelectric Pca21. Ferroelectricity in
hafnia is proper and follows the conventional soft mode theory. The shallow
depth of the $Pbcn$ soft mode potential well is responsible for polarization
switching under applied electric field.

###### Acknowledgements.

A.R. and R.E.C thank P. Zubko for helpful discussions. This work is supported
by U. S. Office of Naval Research Grant N00014-20-1-2699, and the Carnegie
Institution for Science. Computations were supported by high-performance
computer time and resources from the DoD High Performance Computing
Modernization Program, Carnegie computational resources, and REC gratefully
acknowledges the Gauss Centre for Supercomputing e.V. (https://www.gauss-
centre.eu/) for funding this project by providing computing time on the GCS
Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (LRZ, www.lrz.de).

## References

  * Böscke  _et al._ [2011] T. S. Böscke, J. Müller, D. Bräuhaus, U. Schröder, and U. Böttger, Ferroelectricity in hafnium oxide thin films, Applied Physics Letters 99, 102903 (2011).
  * Boscke  _et al._ [2011] T. Boscke, J. Muller, D. Brauhaus, U. Schroder, and U. Bottger, Ferroelectricity in hafnium oxide: Cmos compatible ferroelectric field effect transistors, 2011 International Electron Devices Meeting , 24.5.1 (2011).
  * Robertson [2005] J. Robertson, High dielectric constant gate oxides for metal oxide si transistors, Reports on Progress in Physics 69, 327 (2005).
  * Hubbard and Schlom [1996] K. J. Hubbard and D. G. Schlom, Thermodynamic stability of binary oxides in contact with silicon, Journal of Materials Research 11, 2757–2776 (1996).
  * Ohtaka  _et al._ [2001] O. Ohtaka, H. Fukui, T. Kunisada, T. Fujisawa, K. Funakoshi, W. Utsumi, T. Irifune, K. Kuroda, and T. Kikegawa, Phase relations and volume changes of hafnia under high pressure and high temperature, Journal of the American Ceramic Society 84, 1369 (2001).
  * Cohen [1992] R. E. Cohen, Origin of ferroelectricity in perovskite oxides, Nature 358, 136 (1992).
  * Cohen [2000] R. Cohen, Theory of ferroelectrics: a vision for the next decade and beyond, Journal of Physics and Chemistry of Solids 61, 139 (2000).
  * Kwei  _et al._ [1993] G. H. Kwei, A. C. Lawson, S. J. L. Billinge, and S. W. Cheong, Structures of the ferroelectric phases of barium titanate, The Journal of Physical Chemistry 97, 2368 (1993).
  * Ishidate  _et al._ [1997] T. Ishidate, S. Abe, H. Takahashi, and N. Môri, Phase diagram of bati${\mathrm{o}}_{3}$, Phys. Rev. Lett. 78, 2397 (1997).
  * Inbar and Cohen [1996] I. Inbar and R. E. Cohen, Comparison of the electronic structures and energetics of ferroelectric ${\mathrm{linbo}}_{3}$ and ${\mathrm{litao}}_{3}$, Phys. Rev. B 53, 1193 (1996).
  * Bakker  _et al._ [1993] H. J. Bakker, S. Hunsche, and H. Kurz, Quantum-mechanical description of the ferroelectric phase transition in ${\mathrm{litao}}_{3}$, Phys. Rev. B 48, 9331 (1993).
  * Megaw [1968] H. D. Megaw, The thermal expansion of interatomic bonds, illustrated by experimental evidence from certain niobates, Acta Crystallographica Section A 24, 589 (1968).
  * Barker and Loudon [1967] A. S. Barker and R. Loudon, Dielectric properties and optical phonons in linb${\mathrm{o}}_{3}$, Phys. Rev. 158, 433 (1967).
  * Scrymgeour  _et al._ [2005] D. A. Scrymgeour, V. Gopalan, A. Itagi, A. Saxena, and P. J. Swart, Phenomenological theory of a single domain wall in uniaxial trigonal ferroelectrics: Lithium niobate and lithium tantalate, Phys. Rev. B 71, 184110 (2005).
  * Lee  _et al._ [2020] H.-J. Lee, M. Lee, K. Lee, J. Jo, H. Yang, Y. Kim, S. C. Chae, U. Waghmare, and J. H. Lee, Scale-free ferroelectricity induced by flat phonon bands in HfO2, Science 369, 1343 (2020).
  * Sang  _et al._ [2015] X. Sang, E. D. Grimley, T. Schenk, U. Schröder, and J. M. LeBeau, On the structural origins of ferroelectricity in HfO2 thin films, Applied Physics Letters 106, 162905 (2015).
  * Shimizu  _et al._ [2015] T. Shimizu, K. Katayama, T. Kiguchi, A. Akama, T. J. Konno, and H. Funakubo, Growth of epitaxial orthorhombic YO1.5-substituted HfO2 thin film, Applied Physics Letters 107, 10.1063/1.4927450 (2015).
  * Liu and Hanrahan [2019] S. Liu and B. M. Hanrahan, Effects of growth orientations and epitaxial strains on phase stability of HfO2 thin films, Physical Review Materials 3, 054404 (2019).
  * Hoffmann  _et al._ [2015] M. Hoffmann, U. Schroeder, C. Künneth, A. Kersch, S. Starschich, U. Böttger, and T. Mikolajick, Ferroelectric phase transitions in nanoscale hfo2 films enable giant pyroelectric energy conversion and highly efficient supercapacitors, Nano Energy 18, 154 (2015).
  * Xu  _et al._ [2021] X. Xu, F.-T. Huang, Y. Qi, S. Singh, K. M. Rabe, D. Obeysekera, J. Yang, M.-W. Chu, and S.-W. Cheong, Kinetically stabilized ferroelectricity in bulk single-crystalline HfO2:Y, Nature Materials 20, 826 (2021).
  * Batra  _et al._ [2017] R. Batra, T. D. Huan, J. L. Jones, G. Rossetti, and R. Ramprasad, Factors Favoring Ferroelectricity in Hafnia: A First-Principles Computational Study, The Journal of Physical Chemistry C 121, 4139 (2017).
  * Materlik  _et al._ [2015] R. Materlik, C. Künneth, and A. Kersch, The origin of ferroelectricity in Hf1-xZrxO2: A computational investigation and a surface energy model, Journal of Applied Physics 117, 134109 (2015).
  * Huan  _et al._ [2014] T. D. Huan, V. Sharma, G. A. Rossetti, and R. Ramprasad, Pathways towards ferroelectricity in hafnia, Physical Review B 90, 064111 (2014).
  * Qi  _et al._ [2020] Y. Qi, S. Singh, C. Lau, F.-T. Huang, X. Xu, F. J. Walker, C. H. Ahn, S.-W. Cheong, and K. M. Rabe, Stabilization of competing ferroelectric phases of ${\mathrm{hfo}}_{2}$ under epitaxial strain, Phys. Rev. Lett. 125, 257603 (2020).
  * Wei  _et al._ [2018] Y. Wei, P. Nukala, M. Salverda, S. Matzen, H. J. Zhao, J. Momand, A. S. Everhardt, G. Agnus, G. R. Blake, P. Lecoeur, B. J. Kooi, J. Íñiguez, B. Dkhil, and B. Noheda, A rhombohedral ferroelectric phase in epitaxially strained hf0.5zr0.5o2 thin films, Nature Materials 17, 1095 (2018).
  * Kaiser  _et al._ [2023] N. Kaiser, Y.-J. Song, T. Vogel, E. Piros, T. Kim, P. Schreyer, S. Petzold, R. Valentí, and L. Alff, Crystal and electronic structure of oxygen vacancy stabilized rhombohedral hafnium oxide, ACS Applied Electronic Materials 5, 754 (2023).
  * Nentwich [2022] M. Nentwich, Structure relations in the family of the solid solution hfxzr(1-x)o2, Zeitschrift für Kristallographie - Crystalline Materials 237, 141 (2022).
  * Cheema  _et al._ [2020] S. S. Cheema, D. Kwon, N. Shanker, R. Dos Reis, S. L. Hsu, J. Xiao, H. Zhang, R. Wagner, A. Datar, M. R. McCarter, C. R. Serrao, A. K. Yadav, G. Karbasian, C. H. Hsu, A. J. Tan, L. C. Wang, V. Thakare, X. Zhang, A. Mehta, E. Karapetrova, R. V. Chopdekar, P. Shafer, E. Arenholz, C. Hu, R. Proksch, R. Ramesh, J. Ciston, and S. Salahuddin, Enhanced ferroelectricity in ultrathin films grown directly on silicon, Nature 580, 478 (2020).
  * Hyuk Park  _et al._ [2013] M. Hyuk Park, H. Joon Kim, Y. Jin Kim, W. Lee, T. Moon, and C. Seong Hwang, Evolution of phases and ferroelectric properties of thin Hf0.5Zr0.5O2 films according to the thickness and annealing temperature, Applied Physics Letters 102, 10.1063/1.4811483 (2013), 242905, https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.4811483/14021942/242905_1_online.pdf .
  * Schroeder  _et al._ [2019] U. Schroeder, C. Hwang, and H. Funakubo,  _Ferroelectricity in Doped Hafnium Oxide: Materials, Properties and Devices_, Woodhead Publishing Series in Electronic and Optical Materials (Elsevier Science, 2019).
  * Zhou  _et al._ [2022] S. Zhou, J. Zhang, and A. M. Rappe, Strain-induced antipolar phase in hafnia stabilizes robust thin-film ferroelectricity, Sci Adv 8, eadd5953 (2022), zhou, Songsong Zhang, Jiahao Rappe, Andrew M eng 2022/11/26 Sci Adv. 2022 Nov 25;8(47):eadd5953. doi: 10.1126/sciadv.add5953. Epub 2022 Nov 25.
  * Delodovici  _et al._ [2021] F. Delodovici, P. Barone, and S. Picozzi, Trilinear-coupling-driven ferroelectricity in ${\mathrm{hfo}}_{2}$, Phys. Rev. Materials 5, 064405 (2021).
  * Aramberri and Íñiguez [2023] H. Aramberri and J. Íñiguez, heoretical approach to ferroelectricity in hafnia and related materials (2023), arXiv:2302.00688 .
  * Campbell  _et al._ [2006] B. J. Campbell, H. T. Stokes, D. E. Tanner, and D. M. Hatch, isodisplace: a web-based tool for exploring structural distortions, Journal of Applied Crystallography 39, 607 (2006).
  * Stokes and Hatch [a] H. T. Stokes and D. M. Hatch, isodistort, isotropy Software Suite, iso.byu.edu (a).
  * Aroyo  _et al._ [2006a] M. I. Aroyo, A. Kirov, C. Capillas, J. M. Perez-Mato, and H. Wondratschek, Bilbao crystallographic server. ii. representations of crystallographic point groups and space groups, Acta Crystallographica Section A 62, 115 (2006a).
  * Aroyo  _et al._ [2006b] M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, A. Kirov, and H. Wondratschek, Bilbao crystallographic server: I. databases and crystallographic computing programs, Zeitschrift für Kristallographie - Crystalline Materials 221, 15 (2006b).
  * Aroyo  _et al._ [2006c] M. I. Aroyo, A. Kirov, C. Capillas, J. M. Perez-Mato, and H. Wondratschek, Bilbao Crystallographic Server. II. Representations of crystallographic point groups and space groups, Acta Crystallographica Section A 62, 115 (2006c).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum espresso, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum espresso toward the exascale, The Journal of Chemical Physics 152, 154105 (2020).
  * Garrity  _et al._ [2014] K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Pseudopotentials for high-throughput dft calculations, Computational Materials Science 81, 446 (2014).
  * Prandini  _et al._ [2018] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations, npj Computational Materials 4, 72 (2018).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Raeliarijaona and Cohen [2022] A. Raeliarijaona and R. E. Cohen, First-principles calculations of raman and infrared spectroscopy for phase identification and strain calibration of hafnia, Applied Physics Letters 120, 242903 (2022), https://doi.org/10.1063/5.0092247 .
  * Raeliarijaona and Cohen [2023] A. Raeliarijaona and R. E. Cohen, Origin of wake-up effect in hafnia (2023), arXiv:2302.01981 .
  * Barabash [2017] S. V. Barabash, Prediction of new metastable HfO2: toward understanding ferro- and antiferroelectric films, Journal of Computational Electronics 16, 1227 (2017).
  * Kersch and Falkowski [2021] A. Kersch and M. Falkowski, New low-energy crystal structures in zro2 and hfo2, physica status solidi (RRL) – Rapid Research Letters 15, 2100074 (2021).
  * Ohtaka  _et al._ [1995] O. Ohtaka, T. Yamanaka, S. Kume, N. Hara, H. Asano, and F. Izumi, Structural analysis of orthorhombic hafnia by neutron powder diffraction, Journal of the American Ceramic Society 78, 233 (1995).
  * Noheda and Íñiguez [2020] B. Noheda and J. Íñiguez, A key piece of the ferroelectric hafnia puzzle, Science 369, 1300 (2020).
  * Stokes and Hatch [2005] H. T. Stokes and D. M. Hatch, findsym: program for identifying the space-group symmetry of a crystal, Journal of Applied Crystallography 38, 237 (2005).
  * Stokes and Hatch [b] H. T. Stokes and D. M. Hatch, findsym, isotropy Software Suite, iso.byu.edu (b).
  * Joh  _et al._ [2022] H. Joh, M. Jung, J. Hwang, Y. Goh, T. Jung, and S. Jeon, Flexible Ferroelectric Hafnia-Based Synaptic Transistor by Focused-Microwave Annealing, ACS Applied Materials & Interfaces 14, 1326 (2022).
  * Park  _et al._ [2017] M. H. Park, T. Schenk, M. Hoffmann, S. Knebel, J. Gärtner, T. Mikolajick, and U. Schroeder, Effect of acceptor doping on phase transitions of hfo2 thin films for energy-related applications, Nano Energy 36, 381 (2017).
  * Bärnighausen [1980] H. Bärnighausen, Group-subgroup relations between space groups : a useful tool in crystal chemistry, MATCH, Communications in Mathematical Chemistry 9, 139 (1980).
  * Kudoh  _et al._ [1986] Y. Kudoh, H. Takeda, and H. Arashi, In situ determination of crystal structure for high pressure phase of zro2 using a diamond anvil and single crystal x-ray diffraction method, Physics and Chemistry of Minerals 13, 233 (1986).
  * Ohtaka  _et al._ [1990] O. Ohtaka, T. Yamanaka, S. Kume, N. Hara, H. Asano, and F. Izumi, Structural analysis of orthorhombic zro2 by high resolution neutron powder diffraction, Proceedings of the Japan Academy, Series B 66, 193 (1990).
  * Haines  _et al._ [1997] J. Haines, J. M. Léger, S. Hull, J. P. Petitet, A. S. Pereira, C. A. Perottoni, and J. A. H. da Jornada, Characterization of the cotunnite-type phases of zirconia and hafnia by neutron diffraction and raman spectroscopy, Journal of the American Ceramic Society 80, 1910 (1997).