# A Machine Learning Approach to Understanding the Physical Properties of
Magnetic Flux Ropes in the Solar Wind at 1 AU

Hameedullah Farooki Institute for Space Weather Sciences, New Jersey Institute
of Technology, University Heights, Newark, NJ, USA  Department of Computer
Science, New Jersey Institute of Technology, University Heights, Newark, NJ,
USA  Center for Solar-Terrestrial Research, New Jersey Institute of
Technology, University Heights, Newark, NJ, USA  Yasser Abduallah Institute
for Space Weather Sciences, New Jersey Institute of Technology, University
Heights, Newark, NJ, USA  Department of Computer Science, New Jersey Institute
of Technology, University Heights, Newark, NJ, USA  Sung-Jun Noh Center for
Solar-Terrestrial Research, New Jersey Institute of Technology, University
Heights, Newark, NJ, USA  Now at Los Alamos National Laboratory, Los Alamos,
NM  Hyomin Kim Institute for Space Weather Sciences, New Jersey Institute of
Technology, University Heights, Newark, NJ, USA  Center for Solar-Terrestrial
Research, New Jersey Institute of Technology, University Heights, Newark, NJ,
USA  George Bizos  Department of Computer Science, New Jersey Institute of
Technology, University Heights, Newark, NJ, USA  Youra Shin Center for Solar-
Terrestrial Research, New Jersey Institute of Technology, University Heights,
Newark, NJ, USA  Jason T. L. Wang Institute for Space Weather Sciences, New
Jersey Institute of Technology, University Heights, Newark, NJ, USA
Department of Computer Science, New Jersey Institute of Technology, University
Heights, Newark, NJ, USA  Haimin Wang Institute for Space Weather Sciences,
New Jersey Institute of Technology, University Heights, Newark, NJ, USA
Center for Solar-Terrestrial Research, New Jersey Institute of Technology,
University Heights, Newark, NJ, USA  Big Bear Solar Observatory, New Jersey
Institute of Technology, 40386 North Shore Lane, Big Bear City, CA 92314, USA

###### Abstract

Interplanetary magnetic flux ropes (MFRs) are commonly observed structures in
the solar wind, categorized as magnetic clouds (MCs) and small-scale MFRs
(SMFRs) depending on whether they are associated with coronal mass ejections.
We apply machine learning to systematically compare SMFRs, MCs, and ambient
solar wind plasma properties. We construct a dataset of 3-minute averaged
sequential data points of the solar wind’s instantaneous bulk fluid plasma
properties using about twenty years of measurements from _Wind_. We label
samples by the presence and type of MFRs containing them using a catalog based
on Grad-Shafranov (GS) automated detection for SMFRs and NASA’s catalog for
MCs (with samples in neither labeled non-MFRs). We apply the random forest
machine learning algorithm to find which categories can be more easily
distinguished and by what features. MCs were distinguished from non-MFRs with
an AUC of 94% and SMFRs with an AUC of 89% and had distinctive plasma
properties. In contrast, while SMFRs were distinguished from non-MFRs with an
AUC of 86%, this appears to rely solely on the $\langle
B\rangle>$5\text{\,}\mathrm{nT}$$ threshold applied by the GS catalog. The
results indicate that SMFRs have virtually the same plasma properties as the
ambient solar wind, unlike the distinct plasma regimes of MCs. We interpret
our findings as additional evidence that most SMFRs at 1 au are generated
within the solar wind, and furthermore, suggesting that they should be
considered a salient feature of the solar wind’s magnetic structure rather
than transient events.

Solar Wind (1534) — Random Forests (1935) — Heliosphere (711)

††journal: ApJ

tablenum SIXtablenum

##  1 Introduction

Magnetic flux ropes (MFRs) are structures of plasma in space characterized by
their spiraling magnetic field lines observed throughout the heliosphere with
a wide range of shapes and sizes. They are believed to be the core structure
of many phenomena in the atmosphere of the Sun such as prominences/filaments
(Gibson, 2018; Liu, 2020). As observed in the solar wind near Earth, MFRs are
usually classified as either magnetic clouds (MCs) or small-scale magnetic
flux ropes (SMFRs) (Hu et al., 2018). Both MCs and SMFRs have been detected in
all locations of the heliosphere at which the necessary measurements are
available, from near the Sun (Zhao et al., 2020; Chen et al., 2020), to near
the Earth (Moldwin et al., 2000; Hu et al., 2018; Nieves-Chinchilla et al.,
2018), to above the ecliptic plane (Chen et al., 2019) all the way out to 8 au
(Chen & Hu, 2020). Although they have similar magnetic structures (hence they
are both referred to as MFRs) and are both observed in the solar wind near
Earth, MCs and SMFRs are very different and may have different origins. MCs
are larger and are the internal magnetic structure of interplanetary coronal
mass ejections (ICMEs; coronal mass ejections traveling through the solar
wind), whereas SMFRs are smaller, found scattered throughout the bulk of the
solar wind, and have an uncertain origin (or multiple origins). MCs and SMFRs
differ significantly in their plasma properties, e.g. MCs have very low
temperature and plasma beta whereas SMFRs do not always have these properties.

MCs were originally reported with a strict definition by Burlaga et al.
(1981), although the usage of the term has become more general over time
(Nieves-Chinchilla et al., 2018). Typically observed a few times a month with
durations on the order of tens of hours at 1 AU, they are widely accepted to
approximately be the internal magnetic structure of many, if not all,
interplanetary coronal mass ejections (ICMEs) (Hu et al., 2014; Chen, 2017;
Nieves-Chinchilla et al., 2018), which are coronal mass ejections departed
from the sun propagating through interplanetary space (Howard & Tappin, 2009;
Kilpua et al., 2017)). ICMEs have a significant impact on space weather (Baker
et al., 2004) and thus there is much interest in their study.

In contrast, SMFRs occur frequently near Earth, with a few hundred observed
monthly on average and a strong solar cycle dependency (Hu et al., 2018).
Unlike MCs, the origin of SMFRs is not firmly established (Borovsky, 2008; Hu
et al., 2018; Chen & Hu, 2022; Sanchez-Diaz et al., 2017b, a, 2019; Rouillard
et al., 2010, 2011). They are also much smaller with a typical duration of
under an hour at 1 AU (Hu et al., 2019). Moreover, their scale sizes range
multiple orders of magnitude and their physical properties vary with size (Hu
et al., 2018). Furthermore, detecting SMFRs is challenging as the
observational signatures are not as clear as those of ICMEs. Understanding
SMFRs is important due to their ubiquitous presence in the solar wind,
potential significance for larger-scale phenomena such as solar wind
acceleration, and their close relationship with MCs. There has been increased
interest in detecting and analyzing SMFRs in recent years, particularly near
the Sun using new data from the Parker Solar Probe (Chen & Hu, 2022).

SMFRs were first presented by (Moldwin et al., 1995, 2000). They identified
one SMFR from measurements of the _Ulysses_ spacecraft and six SMFRs from two
other spacecraft, _Wind_ and IMP 8. They found that SMFRs were very similar to
MCs but much smaller and with some different plasma properties. Multiple
statistical studies were conducted in later years, but they were typically
limited to no more than a few hundred samples (Hu et al., 2018). More
recently, (Zheng et al., 2017) developed a novel automated detection algorithm
based on the Grad-Shafranov (GS) reconstruction technique that could detect
tens of thousands of SMFRs from about twenty years of _Wind_ observations.
This event list was then formally published by (Hu et al., 2018). Since then,
the GS-based automated detection methodology has been extended to various
spacecraft throughout the heliosphere such as _Ulysses_ and the Parker Solar
Probe (Chen et al., 2019, 2020; Chen & Hu, 2022).

Machine learning is a part of artificial intelligence with algorithms for
fitting models to data that allow computers to learn patterns in the data and
produce predictions for new data or provide insights about the current data
(Alpaydin, 2020). Machine learning has seen increased use in a wide range of
disciplines including space physics (e.g. Abduallah et al. (2022); Liu et al.
(2020); Vech & Malaspina (2021); Zewdie et al. (2021); Richardson & Cohen
(2021); Roberts et al. (2020); Raheem et al. (2021)). This is because machine
learning is a versatile tool useful when dealing with large amounts of data.
However, there have been limited attempts at applying machine learning to
studying MFRs in the solar wind. Recently, there have been a few machine
learning applications to MCs: dos Santos et al. (2020) used a deep
convolutional neural network to detect signatures of simple MCs in ICMEs, and
Narock et al. (2022) also used a deep convolutional network in order to
identify the axis of a known MC. Nguyen et al. (2019) developed a deep
learning method to detect ICMEs in general. Reiss et al. (2021) used machine
learning algorithms to predict the $B_{z}$ component of an ICME. However,
there is still much room for further application of machine learning to MCs.
Furthermore, to our knowledge, there have been no attempts to apply machine
learning to SMFRs.

In this paper, we apply a machine learning technique to better understand the
importance of various physical properties during the presence (or lack
thereof) of MCs and SMFRs. We do this by using a machine learning algorithm
(known as random forests) to learn the probability distribution of a sample
corresponding to various MFR categorizations (such as SMFR, MC, and non-MFR)
based only on point-in-time physical properties. We use the detected events in
the existing catalogs as training data. The purpose of our application of
machine learning is not to detect MFRs, as point-in-time data is insufficient
to properly detect MFRs (for example, detection algorithms use multiple points
in time to detect a rotation in the magnetic field direction, not just looking
at the direction at a single point in time). Furthermore, both MC and SMFR
catalogs and detection algorithms are already available (see Section 2).
Instead, our machine learning approach makes it possible to systematically
compare the physical properties of the solar wind under various MFR
conditions.

The contribution of our paper is different from previous studies. Hu et al.
(2018) compared the physical properties of SMFRs in the fast and slow solar
wind and analyzed the overall statistics of SMFR properties. But a detailed
and systematic comparison of SMFR properties to the background solar wind
(when no SMFR is detected), as well as between SMFRs and MCs, is still needed
and addressed in this paper. Although we use a machine learning algorithm to
distinguish between MFRs and non-MFRs, our ultimate purpose is not to detect
MFRs but to analyze the differences between MFR categorizations. Previous
machine learning-based studies (e.g. Camporeale et al. (2017) and Li et al.
(2020)) have introduced various methods to accurately classify solar wind
regimes based on in-situ measurements. However, our study differs in its focus
on MFR categorizations (especially SMFR) and its focus on the strength of the
physical properties of the solar wind as probability estimators of the
categorizations rather than the overall accuracy of the classifier.

The rest of this paper is structured as follows. First, in Section 2, we
describe in detail the input data and MFR event lists that we utilize. We also
provide a brief statistical analysis of the dataset. Then, in Section 3, we
describe the machine learning techniques that we use, including the
classification algorithm and feature selection procedure. Next, in Section 4,
we present the results of our experiment. Finally, in Section 5, we discuss
the results and prospectives for future work.

##  2 Data

###  2.1 Input Features

The fundamental physical aspects of the solar wind that we analyze are the
magnetic field and plasma properties. Firstly, the magnetic field is the most
relevant due to its fundamental role in defining what an MFR is (that is, a
plasma structure with twisted magnetic field lines), although the single
timestamp measurements that we use don’t individually contain information
about the spatial structure or temporal change. The atomic composition of the
solar wind plasma is mostly of protons and then of alpha particles, so we
narrow down our focus to those two. Since the data availability and count
statistics for protons (which dominate the solar wind) is higher than those of
alpha particles, we only use proton parameters rather than alpha particle
parameters, other than the alpha/proton number ratio to indicate the abundance
of alpha particles. Because of data gaps in the electron measurements, we
opted to leave out electron data for this study.

All input features come from observations from the _Wind_ spacecraft provided
by NASA’s Space Physics Data Facility. The time period of the data is from
1996 to 2016. The magnetic field parameters are retrieved from the Wind
Magnetic Field Instrument (MFI) dataset (Lepping et al., 1995) at 1 minute
cadence, and the plasma parameters are retrieved from the Soar Wind Experiment
(SWE) dataset (Ogilvie et al., 1995). We use vector parameters in the
Geocentric Solar Ecliptic (GSE) coordinate system. The magnetic field
parameters are the three components of the magnetic field vector $\vec{B}$
($B_{x}$, $B_{y}$, and $B_{z}$) converted to altitude and azimuth
($B_{\theta}$ and $B_{\phi}$), as well as the magnetic field strength $B$,
which is provided by the MFI data averaged from higher cadence measurements.
The plasma parameters are the three components of the proton bulk flow
velocity vector $\vec{u_{p}}$ ($u_{p,x}$, $u_{p,y}$, and $u_{p,z}$) converted
to altitude and azimuth ($u_{p,\theta}$ and $u_{p,\phi}$), proton bulk speed
$u_{p}$, proton number density $n_{p}$, and the alpha/proton ratio
$n_{\alpha}/n_{p}$. Rather than including the absolute value of the proton
temperature, we include the ratio between the observed and expected
temperatures based on the proton speed (Richardson & Cane, 1995). Although we
ignore temporal variations, we included the directional components of the
vector quantities as there may be an overall tendency to different directions.
We converted to spherical coordinates to separate the magnitude from the
direction. Since the Wind/SWE data has an inconsistent cadence, all of the
data is resampled to 3-minute averages. We excluded all data points in which
any of the features had no data in the resampling bins.

In addition to the fundamental parameters, relevant derived parameters are
added to the dataset in case they provide a stronger relationship. Proton gas
pressure is added, as it is important property of the solar wind plasma:

$P_{\mathrm{gas},p}=n_{p}k_{B}T_{p}$ (1)

where $k_{B}$ is the Boltzmann constant. Another form of plasma pressure is
the dynamic pressure from the motion of the protons:

$P_{\mathrm{dyn},p}=\frac{1}{2}m_{p}n_{p}u_{p}^{2}$ (2)

where $m_{p}$ is proton mass. The proton beta is a significant parameter
because one of the most important properties of MFRs, especially larger ones,
is that they tend to have a low beta (Klein & Burlaga, 1982; Hu et al., 2018).
Using the gas pressure and magnetic field strength, the proton beta is
calculated:

$\beta_{p}=\frac{P_{\mathrm{gas,p}}}{B^{2}/(2\mu_{0})}$ (3)

Finally, two other plasma properties important in the solar wind are
Alfv$\acute{\mathrm{e}}$n velocity $v_{A}$ and Alfv$\acute{\mathrm{e}}$n Mach
number $M_{A}=u_{p}/v_{A}$. $v_{A}$ plays an important role in
magnetohydrodynamic fluctuations (Mullan & Smith, 2006). It is conceivable
that differences in the density and magnetic field of MFRs may result in
significant differences in $v_{A}$. This may also affect $M_{A}$, which is an
important indicator of the level of dominance of the magnetic field, plays an
important role in shocks (Sundberg et al., 2017), and has an impact on the
magnetosphere (Lavraud et al., 2013). We include them both, calculating the
Alfv$\acute{\mathrm{e}}$n velocity with the equation:

$v_{A}=\frac{B}{\sqrt{2m_{p}n_{p}\mu_{0}}}$ (4)

###  2.2 Output Labels

The labels are based on separate sources for MCs and SMFRs. MC labels come
from the ICME catalog developed by (Nieves-Chinchilla et al., 2018), which is
available at (https://wind.nasa.gov/ICME_catalog/ICME_catalog_viewer.php).
This catalog is based on previous catalogs as well as additional visual
inspection. Most of the ICMEs have either well-defined MCs or similar
structures (Nieves-Chinchilla et al., 2018), so we simply use the magnetic
object start and end times provided by the catalog, regardless of the magnetic
object type.

SMFRs come from the event list on https://fluxrope.info developed by (Hu et
al., 2018) using the GS automated detection algorithm. We use it because it is
the most comprehensive list of SMFRs from the _Wind_ spacecraft. However, we
found that most MCs were observed as a series of shorter SMFRs due to the
limited window sizes employed by Hu et al. (2018). Therefore, we exclude
events overlapping with MCs (approximately 5% of all SMFR events).

The SMFRs were identified by (Hu et al., 2018) by trying to identify flux rope
structures using segments of the time series observations of the solar wind
from the _Wind_ spacecraft. In contrast, this paper focuses on the properties
of individual data points, not structures measured over time, in order to
focus on the physical properties associated with MFRs. However, there is a
caveat: In order to avoid detecting mere fluctuations as SMFRs, Hu et al.
(2018) excluded SMFR candidates with an average magnetic field strength below
$5\text{\,}\mathrm{nT}$. As will become apparent later, this poses an issue
for analyzing the physical properties of SMFRs.

We distinguish between SMFRs of durations under one hour (labeled short SMFR
or SSMFR) and over one hour (labeled long SMFR or LSMFR). The reason for the
distinction is that most SMFRs are under an hour in duration, so long SMFRs
are unusual. It is possible that LSMFRs and SSMFRs have different origins and
impacts on the geospace system. Furthermore, Hu et al. (2018) have noted
physical differences between short SMFRs and long SMFRs. Finally, long SMFRs
may be expected to have more in common with MCs than short SMFRs due to their
more comparable size. Therefore, it is interesting to compare our results for
short SMFRs and long SMFRs.

In order to generate the binary labels from the combined event lists, we mark
each timestamp contained within any MFR from either list as positive and the
rest as negative. We also add categorical labels describing whether the MFR of
a positive timestamp was an SMFR or an MC and a duration label providing the
MFR’s duration. These labels are not directly used by the machine learning
models, but rather they are used to generate the various binary classification
tasks described in Section 3 before being removed from the dataset.

All samples not included in either catalog were marked as non-MFRs. Therefore,
a limitation of this study is that we assume that the event lists that we use
are comprehensive and exclusive. That is, we assume that all of the events
really were MFRs with correct and precise time boundaries and that all of the
time periods with no events in either list contained no MFR whatsoever.
Nevertheless, the overall differences in statistical properties should be
reasonably reliable assuming that the event lists are mostly correct and
complete. As we will see, the main impact is that the difference between non-
MFRs and SMFRs is difficult to determine.

Figure 1 (left) shows SMFRs on the date 2016-01-18. The shaded regions contain
SMFRs according to the catalog. Each feature is plotted in its own subplot.
Our dataset contains not entire regions, but individual data points, evenly
spaced. For a given sample, there are 15 values corresponding to the 15
features and the aforementioned labels based on what region the sample’s
timestamp falls under. Figure 1 (right) shows an MC that was observed between
dates 2016-01-19 and 2016-01-20, in the blue-shaded region. Before and after
the MC are SMFRs. In between the MC and SMFRs are unshaded regions, which were
not included in either catalog. The samples falling under these regions are
assumed to contain no MFR and are thus marked non-MFRs.

Figure 1: (left) Data for the date 2016–01-18. Shaded regions are SMFR
intervals. Data gaps in individual values are represented by breaks in the
line plots. (right) Data for dates 2016-01-19 and 2016-01-20. Same format as
(left), except the region shaded blue contains an MC.

###  2.3 Statistical Analysis

The final dataset contains 15 features (described previously) and $2.8\times
10^{6}$ samples. Of these, $1.9\times 10^{6}$ (69%) contain no MFR, $7.7\times
10^{5}$ (31%) contain any MFR, $7.3\times 10^{5}$ (26% of the full dataset)
contain an SMFR (including $3.0\times 10^{5}$ SSMFR and $4.3\times 10^{5}$
LSMFR), and $1.4\times 10^{5}$ (5% of the full dataset) contain an MC. The
total number of MFRs is 63807. 63461 of them are SMFRs (from the GS-based
list, not overlapping with an MC) and the remaining 346 are MCs (from the ICME
list). Of the SMFRs, 51,574 are SSMFRs and 11,887 are LSMFRs.

Figure 2: Correlation heatmap. Shows Spearman correlation between features
including raw measurements and derived parameters. Higher positive values
indicate a co-increasing relationship, whereas lower negative values indicate
an inverse relationship. Near-zero indicates a lack of any monotonously
decreasing or increasing relationship.

Considering correlations within the dataset is important both for
understanding the data and because correlations can impact feature importance
scores (since a machine learning model could use a threshold on either of two
variables that have a strong correlation to get the same split of data
samples). Figure 2 shows a heatmap of the pairwise correlation matrix of all
the features. The correlations are calculated using the Spearman correlation
metric, which is similar to the Pearson correlation metric but shows any
monotonically increasing (positive correlation)/decreasing (negative
correlation) relationship rather than only linear relationships. The
correlation matrix is symmetric, so it is the same across the diagonal. The
diagonal contains all ones since each feature is perfectly correlated with
itself. While many feature pairs have almost no correlation, there are some
significant correlations. The strongest ones should be kept in mind when
interpreting the results.

##  3 Methodology

###  3.1 Classification Tasks

We define six binary classification tasks to be solved by our model based on
the data: SSMFR-NMFR, LSMFR-NMFR, SSMFR-LSMFR, MC-NMFR, MC-SSMFR, and MC-
LSMFR. Each task consists of estimating the likelihood of a sample’s true
label being positive based on a subset of the data marked as positive merged
with another subset of the data marked as negative. For example, SSMFR-NMFR
has samples in SMFRs of duration under 1 hour as the positive class and non-
MFR samples as the negative class. Unmarked labels are excluded from the
dataset for the corresponding task. The different tasks extensively explore
the MFR categorization of a point in time and are defined so that every pair
is compared.

###  3.2 Algorithm

We adopt the random forest classification algorithm (Breiman et al., 1984) to
generate our machine learning model. This algorithm is well-suited because it
tries random subsets of the features on random subsets of the data to build an
ensemble model. The random sampling of features is to thoroughly explore
different possible feature combinations to give them a chance to demonstrate
their predictive power and importance with respect to the labels. On the other
hand, the random subsetting of the data, or bootstrapping, is so that the
trees are more diverse. It is widely used in many domains and has been used
successfully for space weather science in particular. For example, (Liu et
al., 2017) used a random forest in order to forecast solar flares in solar
active regions.

In simple terms, the random forest algorithm learns by using a training set of
input data vectors and corresponding output categorizations, or labels, to
generate an ensemble of randomized decision trees that collectively vote on
the correct label of new input. The model is trained as follows: Given a
training set consisting of $N$ vectors with $M$ features
$X=x_{1},x_{2},\cdots,x_{N-1},x_{N}$ with labels
$Y=y_{1},y_{2},\cdots,y_{N-1},y_{N}$, generate $B$ binary trees
$T_{1},T_{2},\cdots,T_{B-1},T_{B}$ according to the following algorithm:

  1. 1.

Sample $N$ elements (with replacement) of the training set, denoted
$(X_{i},Y_{i})$.

  2. 2.

Train binary decision tree $T_{i}$ on $(X_{i},Y_{i})$. At each split, a subset
of $\sqrt{M}$ features are selected to compare. The optimal feature for
splitting is selected by maximizing the resulting information gain or
minimizing the Gini impurity measure.

In our case, we limit the number of trees to 64, the tree depth to 8, and the
maximum number of samples used for training each tree to 10%. Tree splitting
is based on the Gini impurity. When training the model, we undersample the
negative class, i.e. we select a random subset of the negative samples to
match the number of positive samples so that the model is not biased toward
negative samples.

###  3.3 Hyperparameter Tuning

There are various parameters that must be set before training a random forest
model, otherwise known as hyperparameters. Although random forests are
generally not overly sensitive to the chosen hyperparameters, the most
important are the number of trees and number of features used per tree.
Furthermore, although not quite a hyperparameter, a tradeoff between the model
complexity and preventing overfitting must be made. This is done by limiting
the size of the trees, such as by limiting the tree depth. Lower tree depth
also entails faster training.

To optimize the hyperparameters, we temporarily set aside 10% of the training
sets as validation sets. We evaluated how the varying parameters affected all
of the tasks. We found that the default recommendation for the number of
features per tree was optimal. Additionally, we experimented with numbers of
trees including 32, 64, and 128. We found that 64 trees provided optimal
performance without unnecessary computational overhead. Furthermore, we found
that a maximum tree depth of 8 provided optimal performance without
overfitting. Finally, we found that only 10% of the samples (randomly selected
with replacement) per tree were needed due to the large size of our dataset.
Therefore, we applied this limit to decrease training time. An added benefit
is that there is more disparity between the samples used to train each tree of
the random forests, thus increasing their diversity.

###  3.4 Feature Importances

After building a model using the random forest algorithm, we use Gini
importance to rank the importance of each feature for that model Breiman et
al. (1984). In order to calculate the importance of each feature, we start by
calculating the non-normalized importance $\hat{I}_{f,i}$ of each feature $f$
for each individual tree $i$. We calculate the sum of the gain $G_{j}$
(calculated by the random forest algorithm) weighted by the number of samples
$n_{j}$ reaching each node $j$ out of the tree’s nodes that split on feature
$f$:

$\hat{I}_{f,i}=\sum_{j}n_{j}G_{j}$

Then, we normalize the feature importances by dividing them by the sum of all
of the importances:

$I_{f,i}=\frac{\hat{I}_{f,i}}{\sum_{i}\hat{I}_{f,i}}$

Finally, we take the mean of the normalized importances from each tree to get
the overall forest’s importance levels for each feature:

$I_{f}=\frac{1}{B}\sum_{i}\hat{I}_{f,i}$

###  3.5 Feature Selection

In order to evaluate the relevance of the features to each task, we evaluate
the model’s performance for each task on a subset of the dataset set aside for
testing, which is not used for training or hyperparameter tuning. An important
consideration is how to split the dataset into testing and training sets. If
we use random sampling, then that results in the test set being drawn from the
same time period as the training set, which can potentially allow the model to
overfit to the time periods without reducing the performance score on the test
set. On the other hand, the solar wind properties and flux rope occurrence
rate change throughout the solar cycle, so the test set should include samples
from throughout the dataset. Therefore, we split the data into 20 sequential
segments and take the first 20% of each of those 20 segments for the test set,
and use the remaining to build the training set, which contains 80% of the
full dataset.

Using the Gini metric of feature importance, we perform feature selection by
iteratively eliminating features in order of least to most important. That is,
for each task, we remove the least important feature from the dataset, then
re-train and re-evaluate the model with the remaining features. The
performance is evaluated using the Area Under the Receiver-Operator Curve
(AUC) score, which is equivalent to the probability that a positive sample
will receive a higher model output than a negative sample (Hanley & McNeil,
1982). We continue to do this (without recalculating the feature importances)
until there is only the most important feature remaining. Then, we plot this
result to show the AUC score for each number of features, demonstrating how
the model performance decreases when increasingly important features are
removed.

A well-known limitation of using Gini importance for feature ranking is that
if some of the features are correlated, then all of the importance can be
given to one of the features at the expense of the other. A commonly used
alternative is permutation importance, which measures the decrease in the
model’s performance if a given feature’s values are shuffled randomly. We have
opted to use Gini importance for the following reasons: (1) Permutation
importance is only meaningful if the model has high performance, but for some
of the cases in this paper, the random forest classifier was not able to
distinguish between the two categories with good performance. In contrast,
Gini importance is still meaningful for a poorly performing model because it
measures the relative information gained from a particular feature. (2)
Permutation performance can sometimes make the issue of correlations even
worse by reducing the importance of both correlated features. (3) By
experimenting with the use of permutation importance instead of Gini
importance, we have found that for the data used in this study, the results do
not differ significantly and the conclusions are unaffected. While more
sophisticated methods exist to determine feature importances for correlated
features, the number of features in our study is not that large and the
correlations between them can be understood physically and kept in mind when
interpreting the results. Therefore, we have opted not to include the results
with permutation importance in the paper. Instead, we discuss the effect of
the correlations by physical reasoning and use them to reach conclusions about
the most significant features in Section 4.

One issue that may arise when interpreting the feature selection results is
determining how many features one needs to avoid losing a statistically
significant amount of performance. It is difficult to evaluate this visually
when the slope is not steep. Our solution to this is as follows: We split the
training set into 10 disjoint training folds of equal size using the same time
segmenting procedure as before. Furthermore, we split the test set into 10
disjoint testing folds of equal size. Then, we iterate 10 times. For each
iteration $i$, we train the model on the $i$th training fold. Then, once the
model is trained, we perform a nested iteration 10 times. For each nested
iteration $j$, we test the model on the $j$th testing fold. We calculate the
AUC score for iteration $ij$ and store it. At the end, there are $10^{2}=100$
AUC scores. This is repeated with the top A features and top B features,
yielding two sets of 100 AUC scores for statistical comparison. We begin with
$A=1$ and $B=2$, then $A=2$ and $B=3$, and so on. Each time, we perform a
Wilcoxon signed-rank test (Wilcoxon, 1947) to determine whether there is a
statistically significant improvement in performance between using the top $A$
features and the top $B$ features. We continue doing this until we encounter a
p-value of greater than 5%. For example, if this occurs when comparing the top
5 vs the top 6 features for a particular task, then the top 5 features are
considered significant for that task.

##  4 Results

Figure 3: Examples of a decision tree from SSMFR-NMFR. Only the first few
levels of the tree are shown. Each task’s random forest has many trees; this
was selected because it started with the top feature for the respective test.
The gini metric shows the impurity after the previous split. Branches to the
left result when the condition is satisfied. The “value” shows the percentage
of NMFRs (left) and percentage of SMFRs (right) assuming that the samples
start off balanced.

We trained the classification model for each of the six tasks using the random
forest algorithm using their respective training sets. The models themselves
consist of many fairly deep trees, so it would be very difficult to visualize
them directly. In Figure 3, we show an example of an individual decision tree
from task SSMFR-NMFR. This illustrates how the decision tree process works for
individual trees. The overall random forest, however, is based on an ensemble
of many trees. Hence, we use statistical methods below to understand the
results.

Table 1: Cross-Task Performance Class | SSMFR-NMFR | LSMFR-NMFR | SSMFR-LSMFR | MC-NMFR | MC-SSMFR | MC-LSMFR
---|---|---|---|---|---|---
NMFR | 66.8% | 66.5% | - | 88.2% | - | -
SSMFR | 97.0% | - | 64.0% | - | 87.6% | -
LSMFR | - | 93.1% | 63.9% | - | - | 79.0%
MC | - | - | - | 83.5% | 82.8% | 78.5%
AUC Score | 86.1% | 85.2% | 66.2% | 94.1% | 92.7% | 86.4%

We present a detailed breakdown of the performance of each task across the
different MFR classes in the rest of Table 1. The table is generated as
follows. First, we iterate through each class (NMFR, SSMFR, LSMFR, SMFR, MC,
MFR). Then, using the initial models trained on the training sets with all of
the features, predict the class of each sample in the test set belonging to
the current class. We calculate the average predicted class $\langle
y^{\prime}\rangle=\frac{1}{N}\sum y^{\prime}_{i}$ where $y^{\prime}_{i}$ is
$1$ if the model classifies a sample as positive and $0$ if it classifies a
sample as negative. For tasks where the current class is a subset of the
positive class, we use $\langle y^{\prime}\rangle$ as the accuracy. If the
current class is a subset of the negative class, we use $1-\langle
y^{\prime}\rangle$ as the accuracy. If the current class is not a subset of
either one, then there is no meaningful accuracy, so we put dashed lines. The
results are recorded in Table 1. Additionally, the resulting AUC scores on the
test sets are displayed in the last row of Table 1 With the exception of
SSMFR-LSMFR, all of them had around 90% AUC, with MC-NMFR exceeding 94%. The
MC-related tasks performed better than the SMFR-related tasks. Higher AUC
scores indicate that the model is able to assign a higher prediction value to
a true sample over a negative sample more often, and thus the features given
to it have more predictive power. The baseline AUC score for random guessing
is 0.5, so in all cases, the model had significant (if not good) performance.
The fact that the model can give meaningful probability estimates indicates
that there is a significant correspondence between the features that it found
and the labels of each task.

The poor performance of SSMFR-LSMFR is important. If big SMFRs tended to be
more similar to MCs, there should have been a big physical difference between
LSMFR and SSMFR and SSMFR-LSMFR should have scored well. In fact, MC-LSMFR
does perform worse than MC-SSMFR, but SSMFR-LSMFR performs worse than either
one of them, suggesting big SMFRs have more in common with small SMFRs than
MCs. The physical properties exhibit only slight differences, as we will see
below. We will also explore the reasons that MC-NMFR, SSMFR-NMFR, and LSMFR-
NMFR perform well.

Table 2: Feature Importances (Significant Features in Red) Rank | SSMFR-NMFR | LSMFR-NMFR | SSMFR-LSMFR | MC-NMFR | MC-SSMFR | MC-LSMFR
---|---|---|---|---|---|---
1 | $B$ | $B$ | $\beta_{p}$ | $\beta_{p}$ | $\beta_{p}$ | $\beta_{p}$
2 | $M_{A}$ | $M_{A}$ | $T_{p}/T_{ex}$ | $B$ | $T_{p}/T_{ex}$ | $T_{p}/T_{ex}$
3 | $v_{A}$ | $v_{A}$ | $P_{\mathrm{gas},p}$ | $M_{A}$ | $P_{\mathrm{gas},p}$ | $M_{A}$
4 | $P_{\mathrm{gas},p}$ | $P_{\mathrm{dyn},p}$ | $u_{p}$ | $v_{A}$ | $M_{A}$ | $B$
5 | $P_{\mathrm{dyn},p}$ | $n_{p}$ | $B$ | $T_{p}/T_{ex}$ | $B$ | $P_{\mathrm{gas},p}$
6 | $n_{p}$ | $\beta_{p}$ | $n_{\alpha}/n_{p}$ | $P_{\mathrm{gas},p}$ | $v_{A}$ | $v_{A}$
7 | $T_{p}/T_{ex}$ | $P_{\mathrm{gas},p}$ | $v_{A}$ | $n_{p}$ | $P_{\mathrm{dyn},p}$ | $n_{p}$
8 | $\beta_{p}$ | $T_{p}/T_{ex}$ | $n_{p}$ | $P_{\mathrm{dyn},p}$ | $n_{p}$ | $P_{\mathrm{dyn},p}$
9 | $u_{p}$ | $u_{p}$ | $M_{A}$ | $u_{p}$ | $n_{\alpha}/n_{p}$ | $n_{\alpha}/n_{p}$
10 | $B_{\theta}$ | $u_{p,\phi}$ | $P_{\mathrm{dyn},p}$ | $n_{\alpha}/n_{p}$ | $u_{p}$ | $u_{p}$
11 | $u_{p,\theta}$ | $B_{\theta}$ | $B_{\phi}$ | $B_{\theta}$ | $u_{p,\phi}$ | $B_{\phi}$
12 | $u_{p,\phi}$ | $n_{\alpha}/n_{p}$ | $u_{p,\theta}$ | $B_{\phi}$ | $u_{p,\theta}$ | $u_{p,\theta}$
13 | $n_{\alpha}/n_{p}$ | $B_{\phi}$ | $u_{p,\phi}$ | $u_{p,\phi}$ | $B_{\phi}$ | $u_{p,\phi}$
14 | $B_{\phi}$ | $u_{p,\theta}$ | $B_{\theta}$ | $u_{p,\theta}$ | $B_{\theta}$ | $B_{\theta}$

Figure 4: (a) Feature elimination results generated using each individual
task’s feature ranking. (b) Same as (a), except rather than dropping features
based on each task’s ranking, drop in order of average importance across
tasks.

After training each task’s model, we calculate the feature importances for
each task, then use it to rank the features, recorded in Table 2. Finally, we
generated the feature selection curves for each task. The least important
feature is dropped iteratively until only the most important feature remains,
and with each subset of features, the model is re-trained and the AUC score is
calculated. Figure 4 (a) shows the plot of the results for each task. In order
to determine which features were significant for each task, we performed the
Wilcoxon tests in the manner described in Section 3.5. Using this information,
we colored the significant features red in Table 2. We also sorted the
features by their average importance across tasks and then conducted the same
process except using this new cross-task ordering. The resulting feature
elimination plot is Figure 4 (b). The top 5 features across tasks were $B$,
$\beta_{p}$, $T_{p}/T_{\mathrm{exp}}$, $M_{A}$, and $P_{\mathrm{gas}}$.

Figure 5: Histogram and partial dependency plots for the top 5 most important
features across tasks.

The feature ranking tells us which features are significant for distinguishing
a sample’s class, but not what values correspond to a higher probability.
These can be understood with histograms and partial dependency plots (PDPs),
plotted in Figure 5. To generate the PDPs, each model is trained with all
features using the task’s training set, then 10,000 random samples from the
task’s test set are selected. Take magnetic field strength $B$ as an example.
For each point on the PDP, the magnetic field strength $B$ of each selected
sample is changed to a new value $B^{\prime}$, and then all the samples are
passed through the task’s model. The mean resulting probability prediction is
represented by a point on the PDP. The change of the value of the PDP shows
how the model’s output depends on a given parameter.

For both SSMFR-NMFR and LSMFR-NMFR, the top feature is the magnetic field
strength $B$. This implies that the magnetic field strength stands out the
most during the SMFRs in the event list compared to other features. The
comparison of their distributions in Figure 5 reveals a sharp change centered
at $5\text{\,}\mathrm{nT}$: Only 2.7% of samples with
$B<$5\text{\,}\mathrm{nT}$$ contain an SMFR, whereas 47.6% of samples with
$B\geq$5\text{\,}\mathrm{nT}$$ contain an SMFR. Apparently, this is due to the
fact that Hu et al. (2018) excluded any MFR event candidates with an average
magnetic field strength below a $5\text{\,}\mathrm{nT}$ threshold in order to
exclude random fluctuations. We have also included the distribution for the
non-MFR samples with the same threshold applied, which results in a very
similar distribution to the SMFR distribution. The PDP shows a very sharp
increase after $5\text{\,}\mathrm{nT}$ from approximately 0.1 or less to a
nearly constant value for SMFR tasks. In contrast to the PDP for the
previously discussed tasks, the MC tasks have a smooth increase of probability
with magnetic field strength. Interestingly, MC-SSMFR and MC-LSMFR have an
_increase_ in probability below 5 nT, suggesting that at values below the
threshold, the model gives a higher chance of the sample being an MC than an
SMFR even though MCs tend to have higher magnetic field strength than SMFRs.
Overall, it seems that there may be many SMFRs with lower magnetic field
strength and that the threshold has imposed a selection bias.

For the MC-related tasks, the most important feature is proton beta. In fact,
judging by the PDP, $\beta_{p}$ is the main factor in distinguishing MCs from
other categories. It is well-known that the proton beta tends to be very low
for ICMEs (e.g. Liu et al. (2005)). MCs tend to have near-0 values which only
occur rarely in the background solar wind. However, the tail of the
distribution cannot be ignored, since observational studies show that not all
MCs have low proton beta for their entire duration (e.g. Pal et al. (2022)).

SSMFRs and LSMFRs also tend to have lower proton beta than NMFRs, although the
difference is not nearly as pronounced as in the case of MCs. The distribution
of SMFR proton beta has been analyzed by Hu et al. (2018). The smaller
difference between MC and SMFR proton beta compared to MC and non-MFR may
explain the lower AUC score for MC-SSMFR and MC-LSMFR. Likewise, it appears
that LSMFRs have slightly lower $\beta_{p}$ than SSMFRs, so MC-LSMFR performs
worse than MC-SSMFR. However, the shape of the distribution for LSMFRs,
SSMFRs, and NMFRs is virtually the same, whereas MC-LSMFR has a completely
different distribution shape, so this does not necessarily indicate a
similarity between the plasma properties of LSMFRs and MCs. In fact, both
SSMFRs and LSMFRs have almost the same $\beta_{p}$ distribution as NMFRs that
have $B>$5\text{\,}\mathrm{nT}$$, suggesting they both share the distribution
of the background solar wind.

Most features besides $B$ do not appear to be significantly different for
SMFRs based on the histograms compared to NMFRs (with the 5 nT threshold
applied). $M_{A}$ was given high importance for SMFR-related tasks, but does
not appear to affect the PDP much compared to $B$. Since it is highly
correlated with $B$ Figure 2, that is probably what resulted in its high
ranking. In contrast, MCs have significantly different shape of the $M_{A}$
distribution, but since the overall distribution has too much overlap in terms
of the range of values, it was not ranked highly in terms of Gini importance.

Beta, temperature, and gas pressure appear to play an important role for
distinguishing MCs from NMFRs and LSMFRs from SSMFRs. As Hu et al. (2018)
pointed out, smaller SMFRs have a wide range of temperatures whereas larger
SMFRs have lower temperatures. That finding is consistent with the PDP, which
shows that SSMFRs are more likely to have higher temperatures than LSMFRs. But
even though LSMFRs have low temperature, their temperature distribution is
shaped completely different from MCs. Unlike MCs, the distribution of LSMFR
temperature is similar to both SSMFR and NSMFR temperature. Accordingly,
temperature only plays a very strong role for MC-related tasks in the PDP. Gas
pressure, related to temperature by $P=nk_{B}T$, tends to be very low for MCs
compared to SMFRs, LSMFRs, or NMFRs, which all have a similar distribution of
gas pressure. However, larger SMFRs appear to have slightly lower gas
pressure.

##  5 Discussion and Conclusions

In this paper, we ranked the informativeness of various physical properties of
the solar wind properties with regards to various MFR-related categorizations
(SSMFR, LSMFR, MC, and NMFR) using a machine learning feature selection
technique. First, we generated a dataset using measurements from the _Wind_
spacecraft. Then, we used 6 binary classification to compare the point-in-time
physical properties of each unique pair of categories. For each task, we used
the random forest algorithm to build an ensemble of decision trees to
distinguish the positive from negative samples and then used the normalized
average decrease in the impurity of each feature to rank the features. After
that, we applied feature selection by plotting the change in model performance
with increasingly reduced sets of top features used for training and
classification.

Our results demonstrated that MCs tend to have very distinctive solar wind
properties compared to the absence of any MFR, such that an excellent
probability estimate of a given timestamp corresponding to an MC instead of a
non-MFR can be given based on plasma and magnetic field measurements from a
single point in time. The proton beta alone is able to give an AUC score of
over 0.9 for distinguishing MCs from non-MFRs whereas the top 3 features
together ($\beta_{p}$, $B$, and $v_{A}$) result in an AUC score of close to
0.94. This makes sense because it is known that intervals of both low
$\beta_{p}$ and $v_{A}$ are most common in the solar wind during ICMEs
(Gosling & Phan, 2013), and that MCs tend to have elevated magnetic field
strength $B$ (Nieves-Chinchilla et al., 2018).

We found that SMFRs in the catalog can be distinguished from non-MFRs with
significant performance. However, the primary distinguishing feature is just
the $B>$5\text{\,}\mathrm{nT}$$ threshold used by Hu et al. (2018) to make the
catalog. As a result, the plasma properties that differ significantly in
distribution from the non-MFRs—$v_{A}$, $M_{A}$, and $P_{\mathrm{gas}}$ in
particular—are strongly correlated with $B$. Therefore, it seems that SMFRs
share the same distribution as the background solar wind in terms of
instantaneous physical properties. A potential area for future research is to
detect flux ropes from multiple-spacecraft observations, eliminating the need
for the threshold, since one can be more confident of a detected flux rope not
being a mere fluctuation.

Our comparison between SMFRs shorter or longer than an hour in duration
(corresponding to small and large: SMFRs tend to be Parker spiral aligned (Hu
et al., 2018), so duration is a good proxy for size) suggests that there is
almost no different between their physical properties. Even though, as noted
by Hu et al. (2018), larger SMFRs have lower proton temperature, proton beta,
and gas pressure, the difference is small compared to how low the temperature
and beta of MCs is. Larger SMFRs have more in common with smaller SMFRs and
the background solar wind than with MCs. The main difference in physical
properties between larger SMFRs and smaller SMFRs is that larger SMFRs have
slightly lower gas pressure, temperature, and proton beta. However, the
distributions of these parameters have the same shape as the distributions of
smaller events, whereas both are very different from MCs. Additionally, the
magnetic field fluctuates more in smaller flux ropes than larger ones,
although that is due to their structure’s size, not their plasma properties.

Based on our results, we hypothesize that most SMFRs at 1 au share the same
plasma properties as the ambient solar wind. This is in contrast to the
typical description of SMFRs as transient phenomena like MCs, which are
distinct plasma regimes in the solar wind. Considering that the SMFR events in
the current catalog with the $5\text{\,}\mathrm{nT}$ threshold applied already
span 26% of the full time, i.e. SMFRs are present at 1 au over one quarter of
the time, this makes sense. In the context of the flux-tube picture of the
solar wind by Borovsky (2008), the above discussion suggests that the flux
tubes filling the solar wind are often twisted, which may have important
ramifications for the interaction between the solar wind and the
magnetosphere. It also provides further support for most SMFRs at 1 au having
a local origin within the solar wind, e.g. by turbulence (Hu et al., 2018).
However, this study is based on the bulk fluid properties of the solar wind.
Analysis of the microscopic properties such as particle distribution functions
may yield different conclusions regarding the origin of the SMFRs.

Future works may extend this feature ranking methodology to data from other
spacecraft in different parts of the heliosphere and compare the change in
results based on variations in ecliptic latitude and distance from the Sun.
Furthermore, this methodology can be applied to other phenomena with various
in-situ time series measurements of which the relative importances are of
interest.

## Acknowledgments

We thank NASA’s Space Physics Data Facility for providing the _Wind_
spacecraft data, and we thank the teams at the University of Alabama in
Huntsville (UAH) and NASA for the SMFR (available at `https://fluxrope.info`)
and MC (available at
`https://wind.nasa.gov/ICME_catalog/ICME_catalog_viewer.php`) catalogs,
respectively. We acknowledge the support of NASA grant 80NSSC20K1282 and NSF
grants AGS-1927578, AGS-2229064, and OPP-2032421. The figures and results in
this study were generated using Python and Python libraries including
matplotlib (Hunter, 2007), SciKit-Learn (Pedregosa et al., 2011), numpy
(Harris et al., 2020), scipy (Virtanen et al., 2020), and pandas (Wes
McKinney, 2010).

## References

  * Abduallah et al. (2022) Abduallah, Y., Jordanova, V. K., Liu, H., et al. 2022, The Astrophysical Journal Supplement Series, 260, 16.  https://doi.org/10.3847/1538-4365/ac5f56
  * Alpaydin (2020) Alpaydin, E. 2020, Introduction to machine learning (MIT press)
  * Baker et al. (2004) Baker, D. N., Daly, E., Daglis, I., Kappenman, J. G., & Panasyuk, M. 2004, Space Weather, 2.  https://doi.org/10.1029/2003SW000044
  * Borovsky (2008) Borovsky, J. E. 2008, Journal of Geophysical Research: Space Physics, 113.  https://doi.org/10.1029/2007JA012684
  * Breiman et al. (1984) Breiman, L., Friedman, J. H., & Olshen, R. A. 1984, Classification and Regression Trees (Boca Raton, Florida: Chapman and Hall/CRC).  https://doi.org/10.1201/9781315139470
  * Burlaga et al. (1981) Burlaga, L., Sittler, E., Mariani, F., & Schwenn, a. R. 1981, Journal of Geophysical Research: Space Physics, 86, 6673.  https://doi.org/10.1029/JA086iA08p06673
  * Camporeale et al. (2017) Camporeale, E., Carè, A., & Borovsky, J. E. 2017, Journal of Geophysical Research: Space Physics, 122, 10, doi: https://doi.org/10.1002/2017JA024383
  * Chen (2017) Chen, J. 2017, Physics of Plasmas, 24, 090501.  https://doi.org/10.1063/1.4993929
  * Chen & Hu (2020) Chen, Y., & Hu, Q. 2020, The Astrophysical Journal, 894, 25.  https://doi.org/10.3847/1538-4357/ab8294
  * Chen & Hu (2022) —. 2022, The Astrophysical Journal, 924, 43.  https://doi.org/10.3847/1538-4357/ac3487
  * Chen et al. (2019) Chen, Y., Hu, Q., & le Roux, J. A. 2019, The Astrophysical Journal, 881, 58.  https://doi.org/10.3847/1538-4357/ab2ccf
  * Chen et al. (2020) Chen, Y., Hu, Q., Zhao, L., et al. 2020, The Astrophysical Journal, 903, 76.  https://doi.org/10.3847/1538-4357/abb820
  * dos Santos et al. (2020) dos Santos, L. F., Narock, A., Nieves-Chinchilla, T., Nuñez, M., & Kirk, M. 2020, Solar Physics, 295, 1.  https://doi.org/10.1007/s11207-020-01697-x
  * Gibson (2018) Gibson, S. E. 2018, Living reviews in solar physics, 15, 1
  * Gosling & Phan (2013) Gosling, J. T., & Phan, T. D. 2013, The Astrophysical Journal, 763, L39.  https://doi.org/10.1088/2041-8205/763/2/l39
  * Hanley & McNeil (1982) Hanley, J. A., & McNeil, B. J. 1982, Radiology, 143, 29.  https://doi.org/10.1148/radiology.143.1.7063747
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357.  https://doi.org/10.1038/s41586-020-2649-2
  * Howard & Tappin (2009) Howard, T. A., & Tappin, S. J. 2009, Space science reviews, 147, 31.  https://doi.org/10.1007/s11214-009-9542-5
  * Hu et al. (2019) Hu, Q., Chen, Y., & le Roux, J. 2019, Journal of Physics: Conference Series, 1332, 012005.  https://doi.org/10.1088/1742-6596/1332/1/012005
  * Hu et al. (2014) Hu, Q., Qiu, J., Dasgupta, B., Khare, A., & Webb, G. M. 2014, The Astrophysical Journal, 793, 53.  https://doi.org/10.1088/0004-637x/793/1/53
  * Hu et al. (2018) Hu, Q., Zheng, J., Chen, Y., le Roux, J., & Zhao, L. 2018, The Astrophysical Journal Supplement Series, 239, 12.  https://doi.org/10.3847/1538-4365/aae57d
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90.  https://doi.org/10.1109/MCSE.2007.55
  * Kilpua et al. (2017) Kilpua, E., Koskinen, H. E., & Pulkkinen, T. I. 2017, Living Reviews in Solar Physics, 14, 1.  https://doi.org/10.1007/s41116-017-0009-6
  * Klein & Burlaga (1982) Klein, L., & Burlaga, L. 1982, Journal of Geophysical Research: Space Physics, 87, 613, doi: https://doi.org/10.1029/JA087iA02p00613
  * Lavraud et al. (2013) Lavraud, B., Larroque, E., Budnik, E., et al. 2013, Journal of Geophysical Research: Space Physics, 118, 1089.  https://doi.org/10.1002/jgra.50145
  * Lepping et al. (1995) Lepping, R., Acũna, M., Burlaga, L., et al. 1995, Space Science Reviews, 71, 207, doi: https://doi.org/10.1007/BF00751330
  * Li et al. (2020) Li, H., Wang, C., Tu, C., & Xu, F. 2020, Earth and Space Science, 7, e2019EA000997, doi: https://doi.org/10.1029/2019EA000997
  * Liu et al. (2017) Liu, C., Deng, N., Wang, J. T. L., & Wang, H. 2017, The Astrophysical Journal, 843, 104.  https://doi.org/10.3847/1538-4357/aa789b
  * Liu et al. (2020) Liu, H., Liu, C., Wang, J. T. L., & Wang, H. 2020, The Astrophysical Journal, 890, 12.  https://doi.org/10.3847/1538-4357/ab6850
  * Liu (2020) Liu, R. 2020, Research in Astronomy and Astrophysics, 20, 165.  https://doi.org/10.1088/1674-4527/20/10/165
  * Liu et al. (2005) Liu, Y., Richardson, J., & Belcher, J. 2005, Planetary and Space Science, 53, 3.  https://doi.org/10.1016/j.pss.2004.09.023
  * Moldwin et al. (2000) Moldwin, M., Ford, S., Lepping, R., Slavin, J., & Szabo, A. 2000, Geophysical research letters, 27, 57.  https://doi.org/10.1029/1999GL010724
  * Moldwin et al. (1995) Moldwin, M., Phillips, J., Gosling, J., et al. 1995, Journal of Geophysical Research: Space Physics, 100, 19903.  https://doi.org/10.1029/95JA01123
  * Mullan & Smith (2006) Mullan, D., & Smith, C. 2006, Solar Physics, 234, 325.  https://doi.org/10.1007/s11207-006-2077-y
  * Narock et al. (2022) Narock, T., Narock, A., Dos Santos, L. F., & Nieves-Chinchilla, T. 2022, Frontiers in Astronomy and Space Sciences, 9, 27.  https://doi.org/10.3389/fspas.2022.838442
  * Nguyen et al. (2019) Nguyen, G., Aunai, N., Fontaine, D., et al. 2019, The Astrophysical Journal, 874, 145, doi: https://doi.org/10.3847/1538-4357/ab0d24
  * Nieves-Chinchilla et al. (2018) Nieves-Chinchilla, T., Vourlidas, A., Raymond, J., et al. 2018, Solar Physics, 293, 1.  https://doi.org/10.1007/s11207-018-1247-z
  * Ogilvie et al. (1995) Ogilvie, K., Chornay, D., Fritzenreiter, R., et al. 1995, Space Science Reviews, 71, 55, doi: https://doi.org/10.1007/BF00751326
  * Pal et al. (2022) Pal, S., Lynch, B. J., Good, S. W., et al. 2022, Frontiers in Astronomy and Space Sciences, 136.  https://doi.org/10.3389/fspas.2022.903676
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Raheem et al. (2021) Raheem, A.-u., Cavus, H., Coban, G. C., et al. 2021, Monthly Notices of the Royal Astronomical Society, 506, 1916.  https://doi.org/10.1093/mnras/stab1816
  * Reiss et al. (2021) Reiss, M. A., Möstl, C., Bailey, R. L., et al. 2021, Space Weather, 19.  https://doi.org/10.1029/2021SW002859
  * Richardson & Cohen (2021) Richardson, D. K., & Cohen, M. B. 2021, Journal of Geophysical Research: Space Physics, 126.  https://doi.org/10.1029/2021JA029689
  * Richardson & Cane (1995) Richardson, I., & Cane, H. 1995, Journal of Geophysical Research: Space Physics, 100, 23397.  https://doi.org/10.1029/95JA02684
  * Roberts et al. (2020) Roberts, D. A., Karimabadi, H., Sipes, T., Ko, Y.-K., & Lepri, S. 2020, The Astrophysical Journal, 889, 153.  https://doi.org/10.3847/1538-4357/ab5a7a
  * Rouillard et al. (2010) Rouillard, A., Davies, J., Lavraud, B., et al. 2010, Journal of Geophysical Research: Space Physics, 115, doi: https://doi.org/10.1029/2009JA014471
  * Rouillard et al. (2011) Rouillard, A., Sheeley, N., Cooper, T., et al. 2011, The Astrophysical Journal, 734, 7, doi: https://doi.org/10.1088/0004-637X/734/1/7
  * Sanchez-Diaz et al. (2019) Sanchez-Diaz, E., Rouillard, A., Lavraud, B., Kilpua, E., & Davies, J. 2019, The Astrophysical Journal, 882, 51, doi: https://doi.org/10.3847/1538-4357/ab341c
  * Sanchez-Diaz et al. (2017a) Sanchez-Diaz, E., Rouillard, A. P., Davies, J. A., et al. 2017a, The Astrophysical Journal, 851, 32, doi: https://doi.org/10.3847/1538-4357/aa98e2
  * Sanchez-Diaz et al. (2017b) —. 2017b, The Astrophysical Journal Letters, 835, L7, doi: https://doi.org/10.3847/2041-8213/835/1/L7
  * Sundberg et al. (2017) Sundberg, T., Burgess, D., Scholer, M., Masters, A., & Sulaiman, A. H. 2017, The Astrophysical Journal, 836, L4.  https://doi.org/10.3847/2041-8213/836/1/l4
  * Vech & Malaspina (2021) Vech, D., & Malaspina, D. M. 2021, Journal of Geophysical Research: Space Physics, 126.  https://doi.org/10.1029/2021JA029567
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61.  https://doi.org/10.25080/Majora-92bf1922-00a
  * Wilcoxon (1947) Wilcoxon, F. 1947, Biometrics, 3, 119.  https://doi.org/10.2307/3001946
  * Zewdie et al. (2021) Zewdie, G. K., Valladares, C., Cohen, M. B., et al. 2021, Space Weather, 19.  https://doi.org/10.1029/2020SW002639
  * Zhao et al. (2020) Zhao, L.-L., Zank, G. P., Adhikari, L., et al. 2020, The Astrophysical Journal Supplement Series, 246, 26.  https://doi.org/10.3847/1538-4365/ab4ff1
  * Zheng et al. (2017) Zheng, J., Hu, Q., Chen, Y., & le Roux, J. 2017, Journal of Physics: Conference Series, 900, 012024.  https://doi.org/10.1088/1742-6596/900/1/012024


# Stressed and sliding ice surfaces liquefy without much heating

Achraf Atila  Sergey V. Sukhomlinov  Martin H. Müser<EMAIL_ADDRESS>saarland.de Department of Material Science and Engineering, Saarland
University, Saarbrücken, 66123, Germany

###### Abstract

The low kinetic friction observed between ice or snow and numerous
counterbodies is commonly attributed to a thin interfacial water layer [1, 2,
3], which is believed to exist because of pressure melting [4], surface
melting [5, 6], or friction-induced heating [7]. However, even the currently
leading theory of frictional melting keeps being challenged, for example, due
to the lack of detectable warming of snow surfaces under a rotating slider at
$-7^{\circ}$C temperature and $1$ m/s sliding velocity despite high
temporospatial resolution [8]. Here we present molecular simulations of ice
interfaces that reveal that ice surfaces liquefy readily without melting
thermally but rather by displacement-driven amorphization, normal-stress
gradients, and tensile in-plane stress. Yet, friction coefficients below 0.01,
as observed during the sliding of hydrophobic solids over ice [9], appear
possible only when the counterfaces are smooth and allow water to slip past
them. Our findings provide fundamental guidelines on how to optimize ice
friction and challenge experimentalists to measure the surface temperature of
ice and snow at minute scales and with unprecedented speed.

Skidding on ice or snow is a well-known phenomenon, often dreaded, sometimes
loved. Several mechanisms leading to low ice and snow friction have been
proposed and become popular [1, 2, 3, 4, 5, 7]. They are all based on the
assumption of self-lubrication through the melting of surface ice. However, no
single theory is conclusive [8]. Pressure melting would require the true
contact between a ski and the ice below it to be unreasonably small to explain
skiing at $-20^{\circ}$C [7]. Although the molecular mobility of sub-0∘C
surface water correlates with ice friction [10], the variation in friction
coefficients with different counterbodies remains unexplained. Melting by
frictional heating must also be questioned given that neither snow surfaces
while being rubbed [6] nor a contact-induced water film on ice produced at
high sliding velocities of $v_{0}=5$ m/s [11] exceeded 0∘C. Thus, either the
main reason for the presence of water in sliding ice interfaces varies from
case to case, or some crucial ice-liquefaction mechanism has not been
accounted for hitherto. One candidate mechanism could be related to Moras et
al.’s hypothesis [12] that sliding ice undergoes layer-by-layer amorphization
as do diamond as well as silicon in incommensurate, self-mated contacts. Their
conjecture is supported by simulations of sliding grain boundaries at
$-120^{\circ}$C [13], which revealed substantial disordering at small
velocity. Furthermore, ice surfaces would lose their load-bearing ability and
highly dissipative plowing would set in if the interface was heated to 0∘C
[14].

Determining the applicability of different mechanisms that produce interfacial
water requires simultaneous analysis of interfacial stresses, temperature, and
structure for diverse initial and boundary conditions in moving, buried
interfaces. This can be achieved by means of molecular dynamics (MD)
simulations when using force fields that accurately reproduce the relevant
thermodynamic, dynamic, and mechanical properties, as the TIP4P/Ice potential
does for water [15, 16, 17, 18, 19, 20]. It allows us to model sliding ice
interfaces effectively and to investigate how sliding induces or maintains
sub-0∘C interfacial water. We start by simulating flat, incommensurate ice-ice
interfaces. Such calculations provide a lower bound for the friction between
ice crystals as they disregard roughness-induced plasticity, ploughing, snow
compression, capillaries, and other phenomena that enhance energy dissipation,
though grain-size distribution and macroscopic roughness have been shown to
scarcely affect kinetic ice-ice friction [21]. To ascertain its load-bearing
ability and the role of hydrophobicity, we also study ice sliding past
corrugated counterfaces of different hydrophobicity.

Figure 1:  Structural evolution of a sliding ice interface. Snapshots of an
incommensurate interface at $T=10$ K (a) before the onset of sliding and (d)
after sliding a distance of $d=20$ nm at a velocity of $v_{0}=1$ m/s. Colors
indicate the potential energy per molecule. (b) Zooms of the interface before
and during initial sliding. (c) Magnitude of the electrostatic field in a unit
cell located one interlayer spacing above a free surface of an ideal crystal
structure. Top-layer atomic positions are included. (e) Negative change of
total system height, $-\Delta h$, as a function of $d$. Full lines reflect an
initial harmonic response w.r.t. $d$ and shear-induced amorphization through a
$\sqrt{d-d_{0}}$ dependence, where $d_{0}$ is the slid distance at the onset
of plasticity.

In the first set of simulations, two ice crystals are brought into contact at
a temperature of $T=10$ K with a small approach velocity $v_{\perp}$ until the
normal force between their misaligned [0001] surfaces vanishes. Localized
zones of roughening, a few Angstroms wide, appear where the potential energy
per molecule is lower than in the crystal, as evidenced in Fig. 1 (a,b). The
low-energy zones arise when the dipoles of surface molecules are aligned with
and thus attracted by the counterbody’s electrostatic field, which is shown in
Fig. 1(c).

Once sliding starts, these low-energy zones act like cold-welded spots causing
plastic deformation in their vicinity, strengthening these surrounding areas
but weakening the originally cold-welded site. In this way, old low-energy
zones disappear, while new ones arise as sliding progresses. Due to the open
structure and low packing of hexagonal ice, dislocations are not needed for
plastic deformation to occur during this process. The instabilities, which
destroy crystalline order, are local, whereby energy releases and associated
temperature bursts are small. Similar dynamics were observed for other
interfaces and another popular water potential, namely SPC/E [22] (see
examples shown in Fig. S1). Since the [0001] surface is the most densely
packed ice surface and the maximum misorientation of $30^{\circ}$ provides the
best possible condition for structural lubricity (SL), we can conclude that SL
will also not occur at other ice-ice interfaces. SL [23, 24, 25] refers to a
state of small, Stokes-like friction, which is caused by the systematic
cancellation of lateral forces when two atomically flat, incommensurate
crystal surfaces slide past each other without invoking instabilities.

Due to the impossibility of structurally lubric ice-ice contacts, the shear
stress can only be small given the presence of sufficiently thick pre-existing
or stress-disordered layers, whose structure (see e.g., Fig. 1(d)) resembles
that of sheared supercooled bulk water. The claim of resemblance is supported
by comparisons of the pair distribution function (Fig. S2), a three-body,
mixed radial, angular distribution function [26] (Fig. S3), and the
concentration of five-coordinated water molecules (Fig. S4). The disordered
zone has many five-coordinated molecules as regular and supercooled water,
causing the liquid to be denser than the crystal. Thus, reductions in the
separation between the two outermost layers $\Delta h$ allow the thickness or
width $w$ of the amorphized zone to be determined.

Figure 2:  Interfacial evolution during sliding. Width $w_{\rm p}$ versus slid
distance $d_{\textrm{p}}$ at different temperatures $T$ both since onset of
plasticity for sliding velocity (a) $v_{0}=5$ m/s and (b) $20$ m/s. Simulation
data are represented by circles ($T=10$ K), downward triangles ($T=100$ K),
crosses ($T=173$ K), squares ($T=233$ K), diamonds ($T=253$ K), and upward
triangles ($T=263$ K), while the dashed lines correspond to the square-root
fits to the data. (c) Increase of amorphous-layer width $\Delta w$ with time
$t$ at $T=263$ K for $v_{0}=5$ m/s (orange diamonds) and $v_{0}=20$ m/s
(yellow squares) as well as at rest and $T=278$ K: strain-free (blue circles)
and under isotropic, in-plane strain (magenta crosses). Dashed lines are drawn
to guide the eye. (d) Shear stress $\tau$ as a function of slid distance $d$.

Ice turns out to have the same proportionality of $-\Delta h_{\textrm{p}}$,
see Fig. 1(e), or $w_{\textrm{p}}$, see Fig. 2(a,b) and movie M1, to the
square-root of the slid distance $\sqrt{d_{\textrm{p}}}$ as diamond and
silicon [27, 12]. Here, the subscript p indicates that quantities are given
relative to the point at which (substantial) amorphization sets in, which is
also where the stiction peak is located, e.g., at $d\approx 3.5$ nm in Fig.
1(e). The $w_{\textrm{p}}\propto\sqrt{d_{\textrm{p}}}$ relation indicates that
the probability for a surface molecule to abandon its crystallographic
position is linear in a distance increment $\Delta d_{\textrm{p}}$ but
inversely proportional to $w_{\textrm{p}}$ and thus that amorphization is
displacement-driven. Further evidence against thermal melting is provided in
Fig. 2(c), which shows that shear-disordering the first nanometer of thickness
of $-10^{\circ}$C cold ice with $v_{0}=5$ m/s takes roughly as long as it
would take at rest to melt the same amount of ice of a fully thermostatted
crystal superheated to $+5^{\circ}$C. Yet, temperature in the sliding system
only rose to at most $-5^{\circ}$C, as shown in Fig. S4(b). Tensile strains,
which tend to be high at the trailing edge of sliding contacts [28], can be
more important than heating since a 1.5% isotropic, in-plane strain at
$+5^{\circ}$C almost doubles the melting rate compared to the unstrained case,
see again Fig. 2(c).

While the amorphization coefficient $\alpha\equiv
w_{\textrm{p}}/\sqrt{d_{\textrm{p}}}$ changes non-monotonically with
temperature—a relative minimum occurs near $T=233$ K, i.e., slightly above the
temperature separating the low- and high-density regimes of supercooled water
[29]—the shear stress decreases continuously with increasing temperature, as
revealed in Fig. 2(c). This is because the effective viscosity, which is
defined as the ratio of shear stress to shear rate, has a strong temperature
dependence [3, 19] as opposed to $\alpha$. Thus, the reduction in viscosity
due to frictional heating appears to be the primary reason why Bowden and
Hughes [7] observed insulating skis to have lower friction than heat-
conducting skis.

In order to relate our results to experiments, shear stresses must be
converted into friction coefficients $\mu$, defined as the ratio of shear
force to normal load. This requires local water-film heights and normal
contact pressures to be estimated, which is a non-trivial task, because they
depend on contact-patch geometries, the squeeze-out dynamics of water, and on
the rate- and potentially scale-dependent ice plasticity or creep [30, 31,
32]. However, once $v_{\perp}$ exceeds $0.1$ mm/s locally, normal pressures
much greater than the quasi-static penetration hardness $p_{\textrm{H}}$ can
occur, which has a value close to 10 MPa at temperatures where skiing and
skating are most readily possible [14], that is, close but not too close to
melting, e.g., at $-10^{\circ}$C, which is the temperature we focus on in the
following.

Figure 3:  Amorphization and friction of nanoindented ice. (a) Molecular
configuration of a slice of ice subjected to nanoindentation at $-10^{\circ}$C
during the initial stage of sliding with $v_{0}=5$ m/s. (b) Zoom of the region
highlighted by the box in (a). Further zooms of other regions during sliding
for adhesive (c) and nonadhesive (d) indenters. (e) Friction coefficient $\mu$
as a function of sliding distance $d$ normalized to the length $L_{x}=46$ nm
of the simulation cell along the sliding direction. Stars indicate the moment
in time at which the snapshots shown in (a–d) were taken. (f) Penetration
hardness ($P_{\text{h}}$) as a function of the indentation velocity using two
normal loads. In all snapshots only O-O bonds are shown for clarity.

To generate realistic local stresses occurring during ice-asperity
interactions without relying on rough continuum-mechanics calculations, we
simulate a rigid, single-sinusoidal corrugated plate that first indents and
then scratches over an initially flat ice surface at $-10^{\circ}$C. Parts of
the final configuration of an indentation process with an adhesive, rigid tip
are depicted in Fig. 3(a,b). It was produced by first applying a nominal
normal pressure of $p_{\textrm{n}}=100$ MPa for 0.6 ns, which was followed by
a 40 ps relaxation at $p_{\textrm{n}}=20$ MPa. The high initial pressure
creates an indentation mark of similar depth to that obtained with lower
pressure over a time span far beyond the scope of molecular simulations.

Corrugated indenters will sink into ice while amorphizing it even when the
contact pressure has fallen below 200 MPa, which is where ice undergoes a
phase transition to form water at $-10^{\circ}$C. Nonetheless, simulated ice
withstands indentation pressures of up to 300 MPa over extended time periods,
when the counterface is flat, i.e., its radius of curvature $R_{\textrm{c}}$
is formally infinite. This indicates that stress gradients accelerate
amorphization and explains why Fig. 3(f) reveals a similar
$v_{\perp}$-dependent $p_{\textrm{H}}$ as that observed experimentally [14],
albeit shifted to much larger $v_{\perp}$. The experimental $R_{\textrm{c}}$
were more than four orders of magnitude larger than that of our corrugated
counterface, $R_{\textrm{c}}\approx 12.2$ nm. While a diminished ability to
withstand stress at smaller scales echoes the inverse Hall-Petch effect, it
cannot be attributed to grain-boundary sliding as the simulated ice is single
crystalline. Instead, it could be due to the enhancement of non-affine
deformations when the scale over which stress varies is shortened. Further
evidence for a scale-dependent $p_{\textrm{H}}$ comes from atomic-force
microscopy. Using tips with $R_{\textrm{c}}\approx 50$ nm, Butt et al. [33]
found $p_{\textrm{H}}$ to be 10% of the macroscopic indentation hardness, that
is, before they reinterpreted their data and potential errors to reduce the
gap between results and expectations to a factor of 2.5.

When sliding at $v_{0}=5$ m/s under $p_{\textrm{n}}=20$ MPa, see Fig. 3(e) and
movie M3, friction is lowest shortly before the tip sinks into the indentation
mark and highest when reaching its bottom, where $\mu=0.6$, similar to values
obtained with an atomic-force microscope and a tip radius of 200 nm [34]. The
time dependence is a consequence of capillary rather than mechanical effects:
even a thin water layer attempts to reduce its surface energy when rough
solids slide past each other. Things repeat on the second stroke (periodic-
boundary conditions mimic pin-on-disk tribometers). However, spatial
variations in the friction decrease.

A hydrophobic counterface behaves like the hydrophilic one, albeit both the
stiction peak and the kinetic friction are halved. The friction now assumes
the small values that would usually be associated with slippery ice, i.e.,
well below 0.1 at a safe distance from the indentation mark. Just before the
downhill motion $\mu$ even becomes formally negative.

The smaller friction of curved, hydrophobic surfaces not only originates from
the finite slip length of water [3], but also because the temporospatial
stress fluctuations, and thus dissipation, are diminished. The reduction is
substantial although the structural differences between hydrophilic and
hydrophobic setups are subtle, see Fig. 3 (c,d). In the given case, the
increase of the mean friction coefficient from a flat to a curved hydrophilic
counterbody is 0.175 (from 0.103 to 0.278) versus 0.088 (from 0.015 to 0.103)
in the hydrophobic case. This leaves a missing 50% difference of 0.175-0.088 =
0.087, which can only originate from the adhesion-enhanced viscoelastic
dissipation caused by the hydrophilic surface near the leading and the
trailing contact edges. Thus, in addition to interfacial water, counterbodies
must be smooth and hydrophobic for ice to have very low friction coefficients,
i.e., capillary effects have to be small.

In conclusion, our simulation results reveal novel explanations for the
(thermo-)dynamics of sliding ice. Firstly, thin, pre-molten ice layers reduce
shear stresses noticeably only during the first few nanometers of sliding,
because sliding-induced amorphous zones quickly become substantially thicker
than pre-existing equilibrium layers. Their indisputable friction-reducing
effect established at high humidity even for solids as hard as steel [35],
could be indirect in the case of ice, possibly by mitigating maximum stresses
of nano-scale asperities so that fewer apserities break off and turn into
abrasives [8]. Secondly, disordering during sliding is primarily driven by
displacement. Thermally assisted liquefaction may become relevant very close
to the melting point. Thirdly, pressure-induced melting, though frequently
dismissed as irrelevant outside of glaciology, matters whenever local
roughness induces large stress gradients, which enhances amorphization.
Fourthly, friction coefficients well below 0.01, which can arise in speed
skating or car tires on black-ice at say $-10^{\circ}$C seem implausible
without substantial water slip, as on carbonaceous surfaces [36] or without
excess water. This is because squeeze-out and recrystallization [37] will
eventually counterbalance shear-induced amorphization, which is estimated to
happen near $\tau=2$ MPa at $-10^{\circ}$C in the supplementay materials (SM).
This would yield $\mu=0.01$, assuming a high, normal stress of 200 MPa and
contact patches whose linear size exceeds 1 $\mu$m. Characteristic patch sizes
are estimated in the SM using prior experimental [38, 39], theoretical [40],
and numerical [41, 42, 43] work. However, the large normal stress rests on the
paradoxical assumption that surfaces must be simultaneously rough and smooth;
rough to yield large $v_{\perp}$ required for a large $p_{\textrm{H}}$ at
sliding velocities of the order of 1 m/s and smooth to avoid stress gradients
inducing amorphization and dissipation. Thus, when ordinary stick conditions
apply, $\mu$ should clearly exceed the lower bound of
$\tau/P_{\textrm{max}}=0.01$.

One shortcoming of prior models lies in equating frictional and thermal
heating. In contrast, rubbing ice produces more defects but fewer vibrations
than thermal random collisions doing the same amount of work. This and other
insights gained from this work enable more reliable predictions of ice
friction than previously feasible and may aid in understanding processes in
related systems, such as the “purely mechanical wear of diamond by
considerably softer lubricated metal cones” [44], through a potentially
overlooked extreme scale dependence of penetration hardness of loosely packed
crystals.

This scale dependence also implies that (hydrophobic) counterfaces that are
locally smooth but mesoscopically rough should have particularly small kinetic
friction. Such a design avoids small-scale plasticity while allowing for
sufficiently large normal velocities during sliding to boost the rate-
dependent hardness.

##  I Methods

The TIP4P/Ice potential [15] is used. It accurately reproduces the properties
relevant to our study, including the pressure dependence of the melting point
[17], water’s density anomaly [16], the pressure and temperature dependence
[19] of the viscosity, the thickness of the water surface layer [18], and,
most important for our purposes, the structure of ball-milled ice [20].

The starting configurations consist of two crystals facing each other with
perfectly parallel [0001] surfaces. They are misaligned by $30^{\circ}$, which
is realized through a $90^{\circ}$ rotation. To approximate the $\sqrt{2}$
ratio of $a$ and $b$ lattice constants while using a square simulation cell
and periodic boundary conditions, fourteen unit cells are used in the $x$ and
ten unit cells in the $y$ direction (before rotation) and stretched in the
$x$-direction by approximately 0.5%, while compressing the $y$-direction by a
similar amount. This makes one water layer, which contains $1344$ water
molecules, have an area of $10.7\times 10.7$ nm2. Each crystal is made up of
$14$ unit cells containing two planes and a rigid outer plane yielding a total
number of 77,952 H2O molecules in the system. A unit cell has the height
$7.35$ Å, which results in a total initial height of $21.8$ nm, which includes
a separation region of $\approx 4.5$ Å .

When sliding ice against sinusoidal indenters, we use $60\times 13\times 22$
unit cells, resulting in initial system dimensions of $46.0\times 11.4\times
15.5$ nm3 at $T=-10^{\circ}$C, and a total number of 287,040 H2O molecules.
The counterbody has a height profile of $h(x,y)=h_{0}\cos(2\pi x/\lambda$),
where $\lambda=46$ nm is the length of the simulation cell parallel to the
sliding direction and $h_{0}=4$ nm. Lennard-Jones interaction sites assigned
to the surface of the indenter are placed in a hexagonal lattice with a
nearest-neighbor spacing of $a_{0}=2.6$ Å. Interactions only take place with
the oxygen atoms with a length of $\sigma=3.5$ Å and $\varepsilon=0.15$
kcal/mol. Hydrophobicity was modeled by cutting off the interaction at the
minimum of the potential at $r_{\textrm{c}}=\sqrt[6]{2}\sigma$ and at
$r_{\textrm{c}}=2.5~{}\sigma$ for the hydrophilic surfaces, which results in
the surface energy gained upon making contact be similar to that of ice-ice
interfaces.

Sliding is imposed by moving the two (rigid) outermost layers at $v_{1,2}=\pm
v/2$, resulting in a relative in-plane velocity of $v$, chosen within a 1–50
m/s range, which is relevant for many situations, in particular winter sports.
In the ice-ice simulations, normal stresses are usually set to zero and
sometimes to 20 MPa, which is close to a recent estimate for the quasi-static
penetration hardness of ice at -10∘C [14], though differences between these
two sets of simulations turn out to be minor. To be minimally invasive, only
the second-outermost layers are thermostatted, which is done relative to the
respective outermost, rigid layer. Simulations are conducted with LAMMPS [45],
the visualization with OVITO [46], and the post-analysis with in-house codes.

The amorphization width $w$ is determined as the average of an upper and a
lower bound. Both are estimated from the number density of amorphous and
hexagonal-ice molecules as identified using the Chill+ algorithm [47] with a
cutoff of $3.4$ Å. The upper bound is located when positioning the interface
between ice and water at the value of $z$ where the number density of
“amorphous molecules” is 0.08 Å-3. As the sample was cut into slices along
$z$-axis with a bin size of $0.5$ Å, the error in the upper bound of $w$
should not exceed $1.0$ Å. The lower bound is located in the same fashion,
however, using “hexagonal-ice molecules” instead of “amorphous molecules.” At
low temperature, this method gives close results as when correlating the
height difference of the two outermost layers with the width $w$. However, the
latter method gives poor results at high temperature and fails to locate a
crystal surface in the absence of a counterbody.

## Data and materials availability:

All data are available in the main text or the supplementary materials

## References

  * [1] R. Rosenberg.  Why is ice slippery?  Phys. Today, 58:50–54, 2005.
  * [2] J. G. Dash, A. W. Rempel, and J. S. Wettlaufer.  The physics of premelted ice and its geophysical consequences.  Rev. Mod. Phys., 78:695–741, 2006.
  * [3] Ł. Baran, P. Llombart, W. R., and L. G. MacDowell.  Ice friction at the nanoscale.  Proc. Natl. Acad. Sci. U.S.A., 119:e2209545119, 2022.
  * [4] J. Thompson.  I. on recent theories and experiments regarding ice at or near its melting-point.  Proc. R. Soc. Lond., 10:151–160, 1860.
  * [5] Faraday.  XXIV. on regelation, and on the conservation of force.  Lond. Edinb. Dublin Philos. Mag., 17:162–169, 1859.
  * [6] B. Slater and A. Michaelides.  Surface premelting of water ice.  Nat. Rev. Chem., 3:172–188, 2019.
  * [7] F. P. Bowden and T. P. Hughes.  The mechanism of sliding on ice and snow.  Proc. R. Soc. A: Math. Phys. Sci., 172:280–298, 1939.
  * [8] J. H. Lever, S. Taylor, A. J. Song, Z. R. Courville, R. Lieblappen, and J. C. Weale.  The mechanics of snow friction as revealed by micro-scale interface observations.  J. of Glaciol., 64:27–36, 2017.
  * [9] Markus Ovaska and Ari J. Tuononen.  Multiscale imaging of wear tracks in ice skate friction.  Tribol. Int., 121:280–286, 2018.
  * [10] B. Weber, Y. Nagata, S. Ketzetzi, F. Tang, W. J. Smit, H. J. Bakker, E. H. G. Backus, M. Bonn, and D. Bonn.  Molecular insight into the slipperiness of ice.  J. Phys. Chem. Lett., 9:2838–2842, 2018.
  * [11] L Bäurle, D Szabó, M Fauve, H Rhyner, and N D Spencer.  Sliding friction of polyethylene on ice: tribometer measurements.  Tribol. Lett., 24:77–84, 2006.
  * [12] G. Moras, A. Klemenz, T. Reichenbach, A. Gola, H. Uetsuka, M. Moseler, and L. Pastewka.  Shear melting of silicon and diamond and the disappearance of the polyamorphic transition under shear.  Phys. Rev. Mater., 2:083601, 2018.
  * [13] I. de Almeida Ribeiro and M. de Koning.  Grain-boundary sliding in ice i${}_{\textrm{h}}$: tribology and rheology at the nanoscale.  J. Phys. Chem. C, 125:627–634, 2021.
  * [14] R. W. Liefferink, F.-C. Hsia, B. Weber, and D. Bonn.  Friction on ice: How temperature, pressure, and speed control the slipperiness of ice.  Phys. Rev. X, 11:011025, 2021.
  * [15] J. L. F. Abascal, E. Sanz, R. García Fernández, and C. Vega.  A potential model for the study of ices and amorphous water: TIP4P/Ice.  J. Chem. Phys., 122:234511, 2005.
  * [16] C. Vega and J. L. F. Abascal.  Relation between the melting temperature and the temperature of maximum density for the most common models of water.  J. Chem. Phys., 123:2056539, 2005.
  * [17] E. G. Noya, C. Menduiña, J. L. Aragones, and C. Vega.  Equation of state, thermal expansion coefficient, and isothermal compressibility for ices ih, II, III, v, and VI, as obtained from computer simulation.  J. Phys. Chem. C, 111:15877–15888, 2007.
  * [18] M. M. Conde, C. Vega, and A. Patrykiejew.  The thickness of a liquid layer on the free surface of ice as obtained from computer simulation.  J. Chem. Phys., 129:014702, 2008.
  * [19] Ł. Baran, W. Rżysko, and L. G. MacDowell.  Self-diffusion and shear viscosity for the TIP4p/ice water model.  J. Chem. Phys., 158:0134932, 2023.
  * [20] A. Rosu-Finsen, M. B. Davies, A. Amon, H. Wu, A. Sella, A. Michaelides, and C. G. Salzmann.  Medium-density amorphous ice.  Science, 379:474–478, 2023.
  * [21] E. M. Schulson and A. L. Fortt.  Friction of ice on ice.  J. Geophys. Res. Solid Earth, 117:B12204, 2012.
  * [22] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma.  The missing term in effective pair potentials.  J. Phys. Chem., 91:6269–6271, 1987.
  * [23] M. Hirano and K. Shinjo.  Atomistic locking and friction.  Phys. Rev. B, 41:11837–11851, 1990.
  * [24] M. H Müser.  Structural lubricity: Role of dimension and symmetry.  Europhys. Lett., 66:97–103, 2004.
  * [25] D. Dietzel, M. Feldmann, U. D. Schwarz, H. Fuchs, and A. Schirmeisen.  Scaling laws of structural lubricity.  Phys. Rev. Lett., 111:235502, 2013.
  * [26] S. V. Sukhomlinov and M. H. Müser.  A mixed radial, angular, three-body distribution function as a tool for local structure characterization: Application to single-component structures.  J. Chem. Phys., 152:0007964, 2020.
  * [27] L. Pastewka, S. Moser, P. Gumbsch, and M. Moseler.  Anisotropic mechanical amorphization drives wear in diamond.  Nat. Mater., 10:34–38, 2010.
  * [28] C. Müller, M. H. Müser, G. Carbone, and N. Menga.  Significance of elastic coupling for stresses and leakage in frictional contacts.  Phys. Rev. Lett., 131:156201, 2023.
  * [29] P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Panagiotopoulos, J. Russo, J. A. Sellberg, H. E. Stanley, H. Tanaka, C. Vega, L. Xu, and L. G. M. Pettersson.  Water: A tale of two liquids.  Chem. Rev., 116:7463–7500, 2016.
  * [30] O. Lahayne, B. Pichler, R. Reihsner, J. Eberhardsteiner, J. Suh, D. Kim, S. Nam, H. Paek, B. Lorenz, and B. N. J. Persson.  Rubber friction on ice: Experiments and modeling.  Tribol. Lett., 62, 2016.
  * [31] T. Tada, S. Kawasaki, R. Shimizu, and B. N. J. Persson.  Rubber-ice friction.  Friction, 11:1534–1543, 2023.
  * [32] B. N. J. Persson and E. C. Tyrode.  Ice breakloose friction.  J. Chem. Phys., 158:0155545, 2023.
  * [33] H.-J. Butt, A. Döppenschmidt, G. Hüttl, E. Müller, and O. I. Vinogradova.  Analysis of plastic deformation in atomic force microscopy: Application to ice.  J. Chem. Phys., 113:1194–1203, 2000.
  * [34] H. Bluhm, D. F. Ogletree, C. S. Fadley, Z. Hussain, and M. Salmeron.  The premelting of ice studied with photoelectron spectroscopy.  J. Phys.: Cond. Matt., 14:L227–L233, 2002.
  * [35] W. C. Hong, K. Fukuda, and S. Liza.  Influences of atmospheric humidity on sliding speed characteristics of dry sliding phenomena.  Tribol. Online, 18(6):339–352, 2023.
  * [36] E. Secchi, S. Marbach, A. Niguès, D. Stein, A. Siria, and L. Bocquet.  Massive radius-dependent flow slippage in carbon nanotubes.  Nature, 537(7619):210–213, 2016.
  * [37] Y. Xu, N. G. Petrik, R. S. Smith, B. D. Kay, and G. A Kimmel.  Growth rate of crystalline ice and the diffusivity of supercooled water from 126 to 262 K.  Proc. Natl. Acad. Sci. U. S. A., 113:14921–14925, 2016.
  * [38] J. H. Dieterich and B. D. Kilgore.  Direct observation of frictional contacts: New insights for state-dependent properties.  Pure Appl. Geophys., 143:283–302, 1994.
  * [39] L. Bäurle, Th.U. Kaempfer, D. Szabó, and N.D. Spencer.  Sliding friction of polyethylene on snow and ice: Contact area and modeling.  Cold Reg. Sci. Technol., 47:276–289, 2007.
  * [40] B. N. J. Persson.  Elastoplastic contact between randomly rough surfaces.  Phys. Rev. Lett., 87:116101, 2001.
  * [41] L. Pei, S. Hyun, J. Molinari, and M. O. Robbins.  Finite element modeling of elasto-plastic contact between rough surfaces.  J. Mech. Phys. Sol., 53:2385–2409, 2005.
  * [42] M. Müser and A. Wang.  Contact-patch-size distribution and limits of self-affinity in contacts between randomly rough surfaces.  Lubricants, 6:85, 2018.
  * [43] J. M. Monti, L. Pastewka, and M. O. Robbins.  Fractal geometry of contacting patches in rough elastic contacts.  J. Mech. and Phys. Sol., 160:104797, 2022.
  * [44] C. A. Brookes, M. P. Shaw, and P. E. Tanner.  Non-metallic crystals undergoing cumulative work-hardening and wear due to softer lubricated metal sliding surfaces.  Proc. Math. Phys. Eng. Sci., 409:141–159, 1987.
  * [45] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Comput. Phys. Commun., 271:108171, 2022.
  * [46] A. Stukowski.  Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool.  Model. Simul. Mat. Sci. Eng., 18:015012, 2009.
  * [47] Andrew H. Nguyen and Valeria Molinero.  Identification of clathrate hydrates, hexagonal ice, cubic ice, and liquid water in simulations: the chill+ algorithm.  J. Phys. Chem. B, 119(29):9369–9376, 2014.

Acknowledgement The authors gratefully acknowledge the Gauss Centre for
Supercomputing e.V. (www.gauss-centre.eu) for funding this project by
providing computing time through the John von Neumann Institute for Computing
(NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).
MHM acknowledges helpful discussions with Lars Pastewka, Bo Persson, and Chris
May.
Funding: This work was supported by the DFG through grant MU 1694/5-2.
Author contributions:
Conceptualization: M.M. Methodology: A.A, M.M, S.S. Theory: M.M. Data
curation: A.A, S.S. Visualization: A.A, M.M, S.S. Computing-time acquisition:
A.A, M.M, S.S. Funding acquisition: M.M. Writing – original draft: M.M.
Writing – review & editing: A.A, M.M, S.S.
Competing interests: The authors declare that they have no competing
interests.

.

Supplementary Materials:
Supplementary Notes 1–5
Supplementary Figs. S1–S4
References ([45, 46, 47])
Movies M1 to M3