# XAI-BayesHAR: A novel Framework for Human Activity Recognition with
Integrated Uncertainty and Shapely Values

1st Anand Dubey  Infineon Technologies
Munich, Germany
<EMAIL_ADDRESS>2nd Niall Lyons  Infineon Technologies
Dublin, Ireland
<EMAIL_ADDRESS>3rd Avik Santra  Infineon Technologies
Irvine, CA
<EMAIL_ADDRESS>4th Ashutosh Pandey  Infineon Technologies
Irvine, CA
<EMAIL_ADDRESS>
###### Abstract

Human activity recognition (HAR) using IMU sensors, namely accelerometer and
gyroscope, has several applications in smart homes, healthcare and human-
machine interface systems. In practice, the IMU-based HAR system is expected
to encounter variations in measurement due to sensor degradation, alien
environment or sensor noise and will be subjected to unknown activities. In
view of practical deployment of the solution, analysis of statistical
confidence over the activity class score are important metrics. In this paper,
we therefore propose _XAI-BayesHAR_ , an integrated Bayesian framework, that
improves the overall activity classification accuracy of IMU-based HAR
solutions by recursively tracking the feature embedding vector and its
associated uncertainty via Kalman filter. Additionally, _XAI-BayesHAR_ acts as
an out of data distribution (OOD) detector using the predictive uncertainty
which help to evaluate and detect alien input data distribution. Furthermore,
Shapley value-based performance of the proposed framework is also evaluated to
understand the importance of the feature embedding vector and accordingly used
for model compression.

###### Index Terms:

human activity classification, representation learning, accelerometer-
gyroscope, XAI-BayesHAR inference, Explainable methods

##  I Introduction

Human activity recognition can assist to automatically save energy in smart
homes, such as heating, ventilation, air conditioning and lighting by
understanding user’s intentions [1, 2]. Several sensor modalities for human
activity recognition are investigated in literature such as camera based,
radar [3, 4, 5, 6, 7], Inertial Measurement Unit (IMU) [8, 9, 10], infrared,
thermal imaging sensors, etc. The machine learning models used for HAR are
optimized using categorical distribution based loss functions such as softmax
[3, 11, 12].

Recently, [13] evaluated algorithms using CNN, LSTM, BiLSTM, MLP and SVM based
architectures for HAR and demonstrated high performance on activity
classification by allowing separability of the human activity in the feature
space. Although the performance metrics are high for these methods, they fail
to provide discrimination between different activity classes in the feature
space. This is a critical requirement for HAR system to work in an open world
classification setting under variations in the input examples due to sensor
noise, or alien environments. Deep metric learning in HAR problems, enables
the ability for both separability and discriminating learning by providing
distinct class clusters. Contrastive loss [14, 15], triplet loss [16] and
quadruplet loss [17] employ this principle for attaining robust metric based
representation learning.

However, real-time continuous human sensing and activity classification has
challenges arising due to missed missed signals or transitioning between known
and unknown activities. To address this, in [18, 19], solutions using deep
variational embedding models for classification utilizing a tracker based on a
Kalman Filter (KF) is introduced. Here, a metric learning based triplet and
quadruplet loss function for the optimization of the encoder-decoder
architecture is implemented. This helped the network to learn the projection
of time series input data into an embedding vector, where data from similar
activities are grouped together while dis-similar activities are far apart.
During inference, the embedding vector from the learned model is fed into a
Kalman filter to track the embedding vector over time, while handling spurious
miss-classifications and smoothening the embedding vectors towards an activity
cluster centroid. Additionally, the tracked feature embedding is classified
into the desired activity class using a K-NN classifier.

Encoder
---
($q_{\phi}$)
Classification
---
Gating
Tracker
---
(KF)
Model
---
Compression
IMU
---
Signal
Filtering
---
Normali -
---
zation
K-NN
---
Classifier
${\mu}$
---
${\sigma}$
---
Shapley
---
Values
FC-BNN
---

${E_{\mu}}$
---

${E_{\sigma}}$
---

(a)
---
(b)
---
Figure 1: Comparison of (a) State-of-the-art and (b) _XAI-BayesHAR_ framework
for HAR. In contrast to SOTA, _XAI-BayesHAR_ enables the framework to follow
complete Bayesian inference over feature extraction and prediction stage.
Additionally, _XAI-BayesHAR_ uses Shapley value for feature understanding and
model compression during training (illustrated in dotted lines).

While the integrated framework leads to an improved performance in the HAR
domain over SOTA [19], the framework uses K-NN classifier which fails to
estimate uncertainty over HAR classification scores. As a result, this leads
to further developments in this paper, where we proposed an end to end
integrated _XAI-BayesHAR_ framework for HAR. For this purpose, an encoder-
decoder architecture is trained, similar to [19], to estimate both mean
embedding vector and its associated variance. This helped the integrated
tracker to estimate feature vector following Bayesian characteristic. Taking
advantage of Bayesian estimate, a 4-layered fully-connected Bayesian neural
network (FC-BNN) is proposed to classify the human activity into desired
class. FC-BNN leverages the estimate from the tracker by taking both mean and
variance as input and thus predicting both classification score with
associated uncertainty. Furthermore, Shapley values based concept of feature
understanding and selection is applied over learned feature embedding. As a
result, this helped to reduce the dimension feature embedding and thus, the
size of FC-BNN using iterative algorithm.

This paper is organized in 5 section where first section gives introduction to
the HAR system and state-of-the-art (SOTA) methods. Later, section 2 gives
mathmatical formulation and understanding on _XAI-BayesHAR_ followed by
training setup in section 3. At the end, both quantitative and qualitative
analysis is done in detail. The contributions of our paper are:

  * •

The uncertainty associated with activity classification during both training
and inference is evaluated to show the reliability and robustness of the
framework.

  * •

The generalization of the proposed framework is demonstrated by showing its
ability to successfully discriminate new unseen highly-correlated targets.

  * •

The importance and advantage of interpretation and explainability to reduce
the dimension of FC-BNN model.

##  II Proposed Framework

Fig. 1(a) illustrates the state-of-the-art (SOTA) pipeline for HAR where an
encoder-decoder architecture is optimized in combination with variational
approximation together with metric learning at latent space and reconstruction
loss at decoder. The learned mean feature embedding from the latent space is
then integrated into a K-NN classifier for activity classification. In
contrast to this, our framework (termed as _XAI-BayesHAR_), shown in Fig.
1(b), follows a encoder-decoder architecture similar to [18, 19]. This aids
the _XAI-BayesHAR_ to follow Bayesian characteristic over latent feature
embeddings. The _XAI-BayesHAR_ feature embedding from the variational encoder
is integrated inside the Kalman filter, followed by our proposed FC-BNN
network. Due to the inherent nature of Kalman filters and state vector being
Bayesian in nature, the tracked and estimated feature embedding also follow
Bayesian characteristics. As a result, FC-BNN is optimized within the _XAI-
BayesHAR_ framework and outputs class scores and its associated uncertainty.
The uncertainty quantification not only helps to trust the prediction over a
specific target class, but also helps to reject the alien sample treated as
out-of-distribution. Finally, a Shapley value based feature understanding is
done during the training of the FC-BNN. This helped the network to weight
important feature embeddings and thus reduce the size of the FC-BNN to half of
original size. The details on _XAI-BayesHAR_ integration and formulation are
described below.

###  II-A XAI-BayesHAR Integration

An encoder is trained to follow variational inference by mapping input data to
a distribution over a plausible latent feature embedding. Thus, it returns
both mean (confidence $E_{\mu}^{m}$) and variance (uncertainty
$E_{\sigma}^{m}$) over the feature embedding. The extracted Bayesian feature
embedding is augmented inside the Bayesian KF tracker. The KF assumes the
state vector as a Gaussian random variable distribution. Thus, the integration
of the latent embedding distribution into the tracker facilitates the
processing in obtaining not only the value of the current state of the human
activity, but also the uncertainty associated with it. This enables complete
Bayesian inference, the estimation of the feature embedding associated with a
human activity class. Here, the variance over the embedding vector is used for
updating the state uncertainty corresponding to the activity class in the KF.
Additionally, due to the nature of the state vector, which carries a
probabilistic distribution, a Mahalanobis distance as the association metric
is used for the association inside the tracker. This acts as a multivariate
Euclidean norm which is a function of both the mean and covariance of the
predicted state vector.

###  II-B FC-BNN Architecture

The estimated Bayesian embedding distribution from the tracker is used for
HAR. In contrast to the SOTA where K-NN based classifier is applied over mean
embedding distribution for activity classification while ignoring variance
over it, our proposed FC-BNN takes advantage of Bayesian estimates from the
tracker and incorporates input uncertainty in the form of mean and variance.
For this purpose, both input and hidden units are designed to follows a
Bayesian formulation to predict the categorical distribution, in contrast to
the point-estimate neural networks (NNs).

Although NNs can act as universal approximation function for complex and non-
linear functions between inputs and outputs, one of the main limitations of
deterministic NNs is that they are fundamentally frequentist in nature. This
can be understood from the basic formulation of a cost function during
training, i.e., negative log likelihood 111This is equivalent to cross-entropy
for categorical and mean square error for Gaussian distribution. The network
is optimized by maximizing likelihood estimates (MLE) over training data
$\mathcal{D}=\left\\{\left(x_{i},y_{i}\right)\right\\}$ given the network
parameters $\boldsymbol{w}.$ Since the network is trying to maximize the
probability of data itself, with minimum data, deterministic NNs will often
over-fit the data and fail to generalize. As an alternative to this is,
instead of calculating the MLE, the maximum a posteriori (MAP) point estimates
can be calculated.

$\displaystyle\boldsymbol{w}^{\mathrm{MLE}}=\operatorname{argmax}_{\boldsymbol{w}}\log
P(\mathcal{D}\mid\boldsymbol{w})$ (1)
$\displaystyle\boldsymbol{w}^{\mathrm{MAP}}=\operatorname{argmax}_{\boldsymbol{w}}\log
P(\boldsymbol{w}\mid\mathcal{D})$

Whereas, both MLE and MAP give point estimates of parameters, they still bring
limitation to quantify uncertainty over NN estimates. Thus, this is covered by
the posterior predictive distribution, $p(y\mid\mathbf{x},\mathcal{D})=\int
p(y\mid\mathbf{x},\mathbf{w})p(\mathbf{w}\mid\mathcal{D})d\mathbf{w}$ in which
the model parameters follow probability distributions and have been
marginalized out over parameter $\mathbf{w}$. This helps to reject unseen data
which is OOD by quantifying associated uncertainty over its prediction.

The FC-BNN formulation of the hidden layers and prediction layer helps to
propagate both aleatoric and epistemic uncertainty caused due to randomness or
error in true estimation from the tracker and lack of model knowledge due to
the limited data set [20]. Due to the distribution nature of the model
parameters, instead of typical direct backpropagation, these weight
distribution parameters are learned through variational inference. This is
done by minimizing the Kullback-Leibler (KL) divergence between
$q(\mathbf{w}\mid\boldsymbol{\theta})$ and the true posterior $p(\mathbf{w})$
with respect to $\theta$ following below formation.

$\displaystyle\mathcal{F}(\mathcal{D},\boldsymbol{\theta})=$
$\displaystyle\mathrm{KL}(q(\mathbf{w}\mid\boldsymbol{\theta})\|p(\mathbf{w}))-$
(2) $\displaystyle\mathbb{E}_{q(\mathbf{w}\mid\boldsymbol{\theta})}\log
p(\mathcal{D}\mid\mathbf{w})$ $\displaystyle=$
$\displaystyle\mathbb{E}_{q(\mathbf{w}\mid\boldsymbol{\theta})}\log
q(\mathbf{w}\mid\boldsymbol{\theta})-\mathbb{E}_{q(\mathbf{w}\mid\boldsymbol{\theta})}\log
p(\mathbf{w})-$ $\displaystyle\mathbb{E}_{q(\mathbf{w}\mid\theta)}\log
p(\mathcal{D}\mid\mathbf{w})$

As it can be seen from Eq. 2, all three formulations are expectation terms
with respect to the variational distribution
$q(\mathbf{w}\mid\boldsymbol{\theta})$. While the first two terms are data-
independent and can be evaluated layer-wise, the last term is data-dependent
and is evaluated at the end of the forward-pass. Due to the multi-variant
probability distribution nature of the model, it’s not possible to compute
gradient during backpropagation. Thus, taking advantage of stochastic sampling
during forward pass and re-parameterization trick during backward pass,
optimization of the network is done. In our case, we initialize model
parameters with a Gaussian distribution parameterized by,
$\boldsymbol{\theta}=(\boldsymbol{\mu},\boldsymbol{\sigma})$ where
$\boldsymbol{\mu}$ is the mean vector of the distribution and $\sigma$ is the
standard deviation vector. For better understanding and generalization of FC-
BNN, the network is optimized and analyzed with respect to the state-of-the-
art (SOTA) i.e., class feature embedding before integration of KF and our
proposed framework i.e., class embedding after temporal smoothening from
integrated KF.

###  II-C Global Interpretability

Furthermore, ethical issues surrounding the transparency of data and lack of
understanding and trust in machine learning frameworks, creates environments
where explainability plays a key role. The ability to interpret and explain
why computationally complex models make particular decisions based on given
input variables, is inherent to a robust, faithful and trustworthy machine
learning platform [21]. Our proposed framework utilizes the SHAP (SHapley
Additive exPlanations) tool, [22] based on the game theoretically optimal
Shapely Values [23]. The outcome of this approach illustrates the contribution
of each feature on the predicted output. In particular, KernelShap was used
throughout to estimate the contributions of each feature value to a
prediction. KernelShap improves the sample efficiency of model-agnostic
estimations of SHAP values, thus improving the interpretability and
explainability of our framework. Additionally, the SHAP values are used in a
closed loop with FC-BNN, as illustrated in Fig. 1, to downsize both input-
dimension and FC-BNN model size.

Require: $D^{m},D^{v},D^{t}$

Baseline Training :

1: Initialize model $\theta^{0}$

2: $\theta^{m}=\operatorname{train}\left(D^{m},\theta^{0},D^{v}\right)$

3: $Accuracy^{m}=\operatorname{evaluate}\left(\theta^{m},D^{t}\right)$

4: Generate SHAP Value

5: maxSHAP = max(SHAP)

6: minSHAP = min(SHAP)

7: while
_$\operatorname{test}\left(D^{m},\theta^{0_{new}},D^{v}\right)\approx\operatorname{test}\left(D^{m},\theta^{0},D^{v}\right)$_
do

8: for ${shap}_{i}$ in $[{shap}_{1},,{shap}_{2},....{shap}_{n}]$ do

9: $\quad drop(I_{i})$ if $shap^{I_{i}}<SHAP_{th}$

10: end for

11: $\theta^{0}$ = $\theta^{0_{new}}$

end while

Algorithm 1 Retraining Procedure using Shapely Values.

Algorithm 1 gives an empirical understanding on the FC-BNN model optimization
using SHAP value in a closed loop. Here, $D^{m},D^{v},D^{t}$ represents the
training set, validation set, and the test set used for the optimization of
the FC-BNN parameter $\theta^{m}$. After the baseline training of the FC-BNN,
SHAP vales are estimated for input feature embeddings. Therefore, based on
global maximum-minimum SHAP values, input features ($I_{i}$) having SHAP vales
greater than a threshold (SHAPth) are selected. In parallel, the hidden units
of FC-BNN are also reduced. This is performed in an iterative way until SHAP
model accuracy decreases.

##  III Training Setup

###  III-A Dataset Preparation

The datasets used throughout this paper are considered from [19]. Six Known
human activities and one unknown activity are represented using IMU
accelerometer and gyroscope sensors, collected by a PSoC 6 Wi-Fi BT
CY8CKIT-062 Pioneer Kit and the CY8CKIT-028-TFT shield [24]. Sensor data was
collected at 100Hz, in an indoor setting. Known training data consisted of
$21875$ samples, a validation size of $9874$ and a testing size $10197$
samples respectively. Unknown activity samples equated to, $9923$. Time series
data were preprocessed using Z-score normalization [25] and a third order
Butterworth filter [26] with a corner frequency of 0.3 Hz to filter noise.

###  III-B Loss Function

In the paradigm of metric learning, selection of both triplet and quadruplet
pairs are essential. These pairs consist of anchor sample ($x_{a}$), i.e., any
random sample, positive sample ($x_{p}$), which is from the same class as the
anchor, and a negative samples ($x_{n}$), which is a sample from any different
class in comparison to the anchor class. The loss function is computed over
feature embedding in the latent space ($z$), which is output from the encoder
$q(.)$. Both triplet and quadruplet based optimization is done using online
hard and semi-hard pairs following a min-max distance learning between
selected pairs, similar to [18].

Eq. 3 gives mathematical understanding on formulation of triplet loss, where
the distance between the anchor and positive samples is minimized forcing
$d(q_{\phi}(x_{a}),q_{\phi}(x_{p}))$ to $0$ and the distance between the
anchor and negative samples is maximized by making
$d(q_{\phi}(x_{a}),q_{\phi}(x_{p}))+\alpha_{margin}$ less than
$d(q_{\phi}(x_{a}),q_{\phi}(x_{n}))$.

$\displaystyle\mathcal{L}_{\text{triplet}}=$
$\displaystyle\max(\|q_{\phi}(x_{a})-q_{\phi}(x_{p})\|^{2}-\|q_{\phi}(x_{a})-q_{\phi}(x_{n})\|^{2}+$
(3) $\displaystyle\alpha_{\text{margin}},0),$

Here, $\|.\|$ is the euclidean distance function and $\alpha_{\text{margin}}$
is a hyperparameter, which defines the boundary condition between the similar
and dissimilar pairs. In contrast to triplet loss, quadruplet loss includes
another negative sample at the cost of another hyperparameter $\alpha_{2}$,
summarized by Eq. 4. This helps the network to have a better inter and intra-
class distance by adding an extra parameter optimization to separate the
negative class from each other. The resulting new loss function is termed as
Quadruplet loss ($\mathcal{L}_{\text{quadruplet}}$) and can be summarized by
Eq. 4. Here, sample $x_{i}$ and $x_{j}$ belong to the same class and represent
an anchor and positive sample, $x_{k}$ and $x_{l}$ belong to two different
classes, which are also not an anchor class.

$\displaystyle\mathcal{L}_{\text{quadruplet}}=$
$\displaystyle\sum_{i,j,k}^{N}\left\|q\left(x_{i},x_{j}\right)^{2}-q\left(x_{i},x_{k}\right)^{2}+\alpha_{1}\right\|{+}$
(4)
$\displaystyle\dot{\sum}_{i,j,k,l}^{N}\left\|q\left(x_{i},x_{j}\right)^{2}-q\left(x_{l},x_{k}\right)^{2}+\alpha_{2}\right\|,$
$\displaystyle s_{i}=s_{j},s_{l}\neq s_{k},s_{i}\neq s_{l},s_{i}\neq s_{k}$

Furthermore, to enable probabilistic inference for the model over learned
feature embeddings, the concept of variational inference is adapted which
performs an approximate Bayesian inference efficiently by having continuous
feature information. The network is optimized by minimizing the upper-bound on
the expected negative log-likelihood of the data, together with the desired
loss function ($\mathcal{L}_{\text{metric}}$). As a result, in this paper, the
overall loss function for both triplet and quadruplet learning can be
summarized by Eq. 5.

$\displaystyle\mathcal{L}_{\text{ XAI-BayesHAR}}$
$\displaystyle=0.7*\mathcal{L}_{\text{reconstruction}}+0.3*(\mathcal{L}_{\text{KL}}+\mathcal{L}_{\text{metric}}),$
(5)

where $L_{\text{KL}}$ is the KL divergence loss at the latent vector to
minimize its deviation from Gaussianity with 0 mean and unit variance.
$L_{\text{MSE}}$ is the mean-squared error of reconstructing the denoised time
series data at each of the decoder, thus termed as reconstruction loss
function. And, $\mathcal{L}_{\text{metric}}$ denotes metric loss function,
which can be either triplet $\mathcal{L}_{\text{triplet}}$ or quadruplet loss
$\mathcal{L}_{\text{quadruplet}}$.

##  IV Results & Discussion

This section investigates on the reliability towards known and unknown target
classes using both qualitative and quantitative methods for _XAI-BayesHAR_ in
comparison to the SOTA. Additionally, the concept of explainability is used to
interpret latent embedding vectors and further used for model compression.

###  IV-A Original Feature Embedding

The variability of pre-trained model parameter distribution is illustrated in
Fig. 2(a) and Fig. 2(b) for triplet and quadruplet based proposed framework,
respectively.

Figure 2: Visualization of model parameter variability distribution for (a)
triplet and (b) quadruplet based framework. The left column represents SOTA
framework and right column corresponds to _XAI-BayesHAR_ where different color
represents different depth of the FC-BNN. Figure 3: A qualitative analysis of
activity classification for class $0$ using FC-BNN for (a) triplet and (b)
quadruplet based framework, where left column corresponds to SOTA and right
one is for _XAI-BayesHAR_. While both SOTA and _XAI-BayesHAR_ classifies
activity into true target class, the uncertainty over their prediction for
SOTA is much higher in comparison to _XAI-BayesHAR_.

The classification accuracy of SOTA using triplet and quadruplet based
optimization is $0.62$ and $0.93$, respectively. The proposed _XAI-BayesHAR_
improves classification accuracy to $0.85$ for triplet based framework and
$0.97$ for quadruplet based framework. To further understand performance of
the different methodologies, variability of learned parameter distribution
from FC-BNN is visualized and analyzed in Fig. 2. The variance over model
parameter reflects optimization of the function with varied independent
multimode distribution. Thus, a uniform variance characteristic could lead to
higher confusion between closely related known or unknown class distribution.
Interestingly, the variability of parameter distribution for triplet based
SOTA (T-SOTA) in the left column of Fig. 2(a) is almost flat in contrast to
the triplet based proposed framework (T-XAI-BayesHAR), illustrated in the
right column in Fig. 2(a). This gives understanding on poor generalization of
the T-SOTA framework in contrast to the T-XAI-BayesHAR for similar or alien
human activity target class. In contract to this, SOTA optimized using
quadruplet loss (Q-SOTA) shows better variability in parameter distribution in
comparison to the T-SOTA, thus results in better performance and
generalization. Furthermore, quadruplet based proposed framework (Q-XAI-
BayesHAR) enhance the performance of Q-SOTA. This hypothesis is analyzed in
the next section using both qualitative and quantitative approach for known
activity class data distribution and unknown activity class data distribution.

Figure 4: Illustration of qualitative analysis of classification accuracy
using pre-trained FC-BNN over unseen activity data-distribution for (a)
triplet and (b) quadruplet based framework, where left column corresponds to
SOTA and right one is for _XAI-BayesHAR_.

####  IV-A1 Known-Data Distribution

Fig. 3 gives a visual understanding of qualitative analysis of activity
classification using FC-BNN applied over feature embedding estimated from SOTA
(left column) and _XAI-BayesHAR_ (right column). While Fig. 3(a) corresponds
to triplet based optimized encoder, Fig. 3(b) represents analysis for
quadruplet based optimized encoder. Additionally, in every subplot, the first
sub-plot from the left represents class prediction variability (uncertainty)
by posterior sample followed by classification score. The x-axis represents
human activity (“Idle”: 0, “Jump”: 1, “Sit”: 2, “Squat”: 3, “Stairs”: 4,
“Stand”: 5, “Walk”: 6) and the y-axis represents confidence scores in the case
of predictive probability and uncertainty in posterior samples.

It is evident that both SOTA and _XAI-BayesHAR_ predicts the true activity
class for both the triplet and quadruplet based framework, the T-SOTA shows
much higher uncertainty for class $2$ in comparison to the Q-SOTA. This is due
to the reason that variability on learned weight parameters for T-SOTA are
uniform, as shown in Fig. 2(a)-left, in contrast to Q-SOTA from Fig.
2(b)-left. Further, the uncertainty over the prediction estimate is reduced to
zero with our proposed integrated framework for T-XAI-BayesHAR and Q-XAI-
BayesHAR. Thus, both Fig. 3 and 2 demonstrates that our framework helps to
learn variability among target class distribution while preserving the target
class score. In addition to it, this also shows that FC-BNN successfully
learns the confusion metric by quantifying uncertainty together with class
prediction scores.

####  IV-A2 Out of Distribution

Furthermore, _XAI-BayesHAR_ is evaluated for a new activity class (kicking),
labelled as $7$, which was not seen by both pre-trained encoder and FC-BNN
during training. Prior to evaluation of classification scores and associated
uncertainty, Pearson coefficient for class $7$ is being calculated against
feature embedding estimated from SOTA and _XAI-BayesHAR_ for all other human
activity, as mentioned in Table I.

TABLE I: Pearson coefficient to measure class similarity between known class distribution and out of class distribution for features from SOTA and _XAI-BayesHAR_. Unknown Class | Known Class Distribution
---|---
Distribution (7) |  $0$ |  $1$ |  $2$ |  $3$ |  $4$ |  $5$ |  $6$
Q-XAI-BayesHAR |  $-0.03$ |  0.52 |  $-0.06$ |  $0.47$ |  $0.33$ |  $0.44$ |  $0.43$
Q-SOTA |  0.49 |  0.99 |  0.47 |  0.78 |  -0.55 |  0.51 |  0.49
T-XAI-BayesHAR |  0.24 |  0.47 |  0.29 |  -0.46 |  0.43 |  -0.18 |  0.14
T-SOTA |  $0.37$ |  $0.98$ |  $0.29$ |  $-0.72$ |  $-0.69$ |  $-0.31$ |  $-0.40$

While kicking class $7$ is treated as OOD, the Pearson coefficient for kicking
class shows strong relation with jumping class $1$ for both triple and
quadruplet based SOTA and _XAI-BayesHAR_. In addition to it, Table I also
shows the decreased correlation between pre-trained class distribution and an
alien class for the proposed framework in comparison to SOTA. This further
demonstrates better separability of the learned feature embedding using _XAI-
BayesHAR_ in contrast to the SOTA.

Similar to Fig. 3, Fig. 4 evaluates the performance of the proposed framework
as an OOD detector using uncertainty estimates from FC-BNN. For this purpose,
the feature embedding of class $7$ activity is extracted from both SOTA and
_XAI-BayesHAR_ and passed to pre-trained categorical FC-BNN. Due to the
inherent nature of categorical based optimization, the network is bound to
predict class label with in bounded class distribution (i.e., in our case
between $0$ to $6$). By taking advantage of FC-BNN, estimates can be rejected
if uncertainty is high.

As discussed before, due to low separability between pre-trained human
activity target class and new alien class, the predictive uncertainty from FC-
BNN shows a uniform uncertainty behavior between class $3,4,5,6$ and making it
not rely on predictive estimates in the case of SOTA. Furthermore, the feature
embedding from class 7 is highly related to other target class, as seen in
Table I. In contrast to this, despite _XAI-BayesHAR_ shows better separability
between pre-trained target class’s, it successfully rejects alien class’s with
high uncertainty between class $1$ and class $3$ for T-XAI-BayesHAR and class
$1$ and class $0$ for Q-XAI-BayesHAR.

####  IV-A3 Local Interpretability

_XAI-BayesHAR_ utilizes SHAP to interpret the contribution of each feature to
the predicted output.

Figure 5: Proposed framework SHAP Force Plot of individual classes (a)
Jumping, (b) Kicking

As outlined in Table I, class $7$ (kicking) is treated as out-of-distribution,
however the Pearson coefficient shows strong relation with class $1$
(jumping). SHAP force plots are used to see how features contributed to the
model’s prediction for a specific observation, namely kicking and jumping.
Interestingly, Fig. 5 (a) jumping and (b) kicking share attributes that
contribute to the output, mean_b and mean_d. However, it can be seen that the
out-of-distribution sample attributes push the prediction value lower, while
mean_h in (a) push the prediction value higher. Therefore, although a strong
Pearson coefficient, through model interpretability, the _XAI-BayesHAR_
successfully demonstrates that the framework can handle classes with a low
inter-class difference.

###  IV-B Reduced Feature Embedding

In the second stage of our experiment, we further analyze the contribution of
each feature embedding and compress the FC-BNN model in a close loop, as
mentioned earlier in Algorithm. 1. Additionally, uncertainty analysis for both
inlier and outlier data distribution is analyzed similar to Fig. 3 and Fig. 4.

Fig. 6 illustrates the contributions made by each feature influence the
outcome of the models in ascending order. Both SOTA and _XAI-BayesHAR_
frameworks are evaluated for triplet and quadruplet. Each sample in the
dataset is represented for each attribute, indicating the importance of that
feature on the prediction. The Shap value impact on the horizontal axis of the
plot indicates the importance of the prediction. Thus, how each class
embedding after temporal smoothening from integrated KF will affect the
overall output. The TVAE and QVAE mean output vector, for illustration
purposes, is labelled mean_a-h.

Figure 6: (a) TVAE, and (b) QVAE SHAP values feature summary of SOTA (left)
and proposed framework (right)

Interestingly, it can be seen in Fig. 6 (a1) and (a2) that mean_e, is the most
significant feature for the T-SOTA and T-XAI-BayesHAR. However, SHAP value
impact in the SOTA is negatively correlated to target classes, while it has a
high and positive impact on the target class in the proposed framework. In
comparing the QVAE, SOTA and proposed in Fig. 6 (b1) and (b2) respectively,
the SHAP value feature influences changes with the introduction of temporal
smoothening. It can be seen that mean_b, mean_h, mean_c has high correlations
with the prediction, while mean_e and mean_g have little relevance to the
output prediction. Through the addition of exploring models using
interpretability and explainable techniques, models are re-trained using
knowledge gained from SHAP value plots. Thus, reducing space and time
complexities, and creating a noise reduced more interpretable model that is
integrated into our closed loop proposed framework.

####  IV-B1 Known-Data Distribution

Leveraging the SHAP value corresponding to each feature value, the input
dimension is reduced to $3$ to that of the original dimension of $16$. This
result in the compression of FC-CNN model by a factor of $95\%$ with $402$
trainable parameters in contrast to the original FC-BNN with $7548$
parameters. Following similar procedure as below, the SOTA and _XAI-BayesHAR_
give classification accuracy of $0.80$ and $0.95$ for triplet and $0.70$ and
$0.82$ for quadruplet based embedding. Interestingly, the triplet based SOTA
and _XAI-BayesHAR_ shows better classification accuracy in comparison to
quadruplet based SOTA and _XAI-BayesHAR_. To better understand this, the tiny
FC-BNN classification accuracy is further evaluated using _XAI-BayesHAR_
uncertainty. Prior to looking into uncertainty, learned variability inside
tiny FC-BNN is evaluated and illustrated in Fig. 7, similar to Fig. 2. The
Fig. 7 shows uniform learned parameter variability for the case of both T-SOTA
and Q-SOTA. This gives an indication for over-fitting of tiny FC-BNN model
using reduced triplet embedding and poor separability of different activity
class.

Figure 7: Visualization of model parameter variability distribution for (a)
triplet and (b) quadruplet based framework. The left column represents SOTA
framework and the right column corresponds to _XAI-BayesHAR_ where different
color represents different depth of the tiny FC-BNN. Figure 8: A qualitative
analysis of pre-trained FC-BNN over unseen activity data-distribution
estimated from (a) SOTA and (b) our _XAI-BayesHAR_ framework. The uncertainty
from _XAI-BayesHAR_ is much higher in comparison to SOTA and help to discard
the prediction.

Further, Fig. 8 shows the behavior of tiny FC-BNN over reduced feature
embeddings for known activity class’s. Clearly, SOTA framework shows much
higher uncertainty between multiple target class. This is due to reason of
uniform variability in the learned model parameters inside the tiny FC-BNN.
Additionally, T-SOTA shows higher confidence over the prediction class in
comparison to the Q-SOTA. The similar trend is followed within _XAI-BayesHAR_
whereas, advantage of our _XAI-BayesHAR_ can be reflected on increased
uncertainty over target class 3 in case of false prediction.

####  IV-B2 Unknown-Data Distribution

In addition to uncertainty over reduced embedding for learned activity class
embedding, both SOTA and _XAI-BayesHAR_ is evaluated for reduced dimension
feature embedding for alien class $7$. Table II shows that the correlation
between the alien class and known activity class for reduced feature embedding
follows the similar pattern as to the original dimension of feature embedding
described in Table I.

TABLE II: Pearson coefficient to measure class similarity between known class distribution and out of class distribution for features before integration of KF and after temporal smoothening with KF Unknown Class | Known Class Distribution
---|---
Distribution (7) |  $0$ |  $1$ |  $2$ |  $3$ |  $4$ |  $5$ |  $6$
Q-_XAI-BayesHAR_ |  0.13 |  0.86 |  0.32 |  0.26 |  -0.28 |  -0.58 |  -0.68
Q-SOTA |  $0.6$ |  $0.86$ |  $0.58$ |  $-0.67$ |  $-0.64$ |  $-0.65$ |  $-0.68$
T-XAI-BayesHAR |  0.86 |  -0.01 |  0.25 |  -0.26 |  -0.28 |  0.14 |  0.5
T-SOTA |  $0.99$ |  $0.22$ |  $0.67$ |  $-0.56$ |  $-0.62$ |  $-0.17$ |  $0.60$
Figure 9: A qualitative analysis of pre-trained FC-BNN over unseen activity
data-distribution estimated from (a) SOTA and (b) our _XAI-BayesHAR_
framework. The uncertainty from _XAI-BayesHAR_ is much higher in comparison to
SOTA and helps to discard the prediction.

Fig. 9 shows the behavior of tiny FC-BNN over reduced feature embeddings for
alien data. Clearly, both T-SOTA and Q-SOTA fail to generalize over new unseen
data in reduced feature embedding and leads to failing to identify outlier
class activity. Whereas, both T-XAI-BayesHAR and Q-XAI-BayesHAR shows high
uncertainty over target class $0$, $2$ and $5$ for T-XAI-BayesHAR, and $0$ and
$1$ for Q-XAI-BayesHAR. This indicates that despite quadruplet based
optimization for both SOTA and _XAI-BayesHAR_ over reduced feature embeddings,
it shows reduced accuracy in comparison to the triplet based framework. In
contrast to the quadruplet based framework, which demonstrates better
separability and robustness against outliers while reducing model dimension.

##  V Conclusion

Our proposed framework, _XAI-BayesHAR_ , enables a complete Bayesian
formulation for HAR using IMU sensors. For this purpose, _XAI-BayesHAR_ takes
advantage of Bayesian formulation over both feature learning and activity
classification. As a result, the framework shows robustness against alien
target class (kicking in our case) and rejects it with high uncertainty over
classification scores. Additionally, the metric based learning and temporal
smoothening used inside _XAI-BayesHAR_ further improves the classification in
comparison to the SOTA by $\approx 35$% and $\approx 5$% for triplet and
quadruplet based _XAI-BayesHAR_ , respectively. Furthermore, the concept of
SHAPLEY value based model explainability is used in a closed loop to compress
the FC-BNN by a factor of $95\%$. The tiny FC-BNN successfully retains the
separability of both known and unknown activity class applied over _XAI-
BayesHAR_ in comparison to the SOTA.

## References

  * [1] V. Garg and N. K. Bansal, “Smart occupancy sensors to reduce energy consumption,” _Energy and Buildings_ , vol. 32, no. 1, pp. 81–87, 2000.
  * [2] A. Santra, R. V. Ulaganathan, and T. Finke, “Short-range millimetric-wave radar system for occupancy sensing application,” _IEEE sensors letters_ , vol. 2, no. 3, pp. 1–4, 2018.
  * [3] Y. Kim and T. Moon, “Human detection and activity classification based on micro-doppler signatures using deep convolutional neural networks,” _IEEE geoscience and remote sensing letters_ , vol. 13, no. 1, pp. 8–12, 2015\.
  * [4] R. J. Javier and Y. Kim, “Application of linear predictive coding for human activity classification based on micro-doppler signatures,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 11, no. 10, pp. 1831–1834, 2014\.
  * [5] M. Zenaldin and R. M. Narayanan, “Radar micro-doppler based human activity classification for indoor and outdoor environments,” in _Radar Sensor Technology XX_ , vol. 9829. International Society for Optics and Photonics, 2016, p. 98291B.
  * [6] R. Hernangómez, A. Santra, and S. Stańczak, “Human activity classification with frequency modulated continuous wave radar using deep convolutional neural networks,” in _2019 International Radar Conference (RADAR)_. IEEE, 2019, pp. 1–6.
  * [7] P. Vaishnav and A. Santra, “Continuous human activity classification with unscented kalman filter tracking using fmcw radar,” _IEEE Sensors Letters_ , 2020.
  * [8] A. E. Minarno, W. A. Kusuma, H. Wibowo, D. R. Akbi, and N. Jawas, “Single triaxial accelerometer-gyroscope classification for human activity recognition,” in _2020 8th International Conference on Information and Communication Technology (ICoICT)_. IEEE, 2020, pp. 1–5.
  * [9] A. K. M. Masum, M. E. Hossain, A. Humayra, S. Islam, A. Barua, and G. R. Alam, “A statistical and deep learning approach for human activity recognition,” in _2019 3rd International Conference on Trends in Electronics and Informatics (ICOEI)_. IEEE, 2019, pp. 1332–1337.
  * [10] M. A. I. Anik, M. Hassan, H. Mahmud, and M. K. Hasan, “Activity recognition of a badminton game through accelerometer and gyroscope,” in _2016 19th International Conference on Computer and Information Technology (ICCIT)_. IEEE, 2016, pp. 213–217.
  * [11] S. Ha and S. Choi, “Convolutional neural networks for human activity recognition using multiple accelerometer and gyroscope sensors,” in _2016 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2016, pp. 381–388.
  * [12] T. Stadelmayer, M. Stadelmayer, A. Santra, R. Weigel, and F. Lurz, “Human activity classification using mm-wave fmcw radar by improved representation learning,” in _Proceedings of the 4th ACM Workshop on Millimeter-Wave Networks and Sensing Systems_ , 2020, pp. 1–6.
  * [13] S. Wan, L. Qi, X. Xu, C. Tong, and Z. Gu, “Deep learning models for real-time human activity recognition with smartphones,” _Mobile Networks and Applications_ , vol. 25, no. 2, pp. 743–755, 2020.
  * [14] Y. Taigman, M. Yang, M. Ranzato, and L. Wolf, “Deepface: Closing the gap to human-level performance in face verification,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2014, pp. 1701–1708.
  * [15] J. Weiß and A. Santra, “One-shot learning for robust material classification using millimeter-wave radar system,” _IEEE sensors letters_ , vol. 2, no. 4, pp. 1–4, 2018.
  * [16] A. Santra and S. Hazra, _Deep learning applications of short-range radars_. Artech House, 2020.
  * [17] W. Chen, X. Chen, J. Zhang, and K. Huang, “Beyond triplet loss: a deep quadruplet network for person re-identification,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 403–412.
  * [18] A. Dubey, A. Santra, J. Fuchs, M. Lübke, R. Weigel, and F. Lurz, “A bayesian framework for integrated deep metric learning and tracking of vulnerable road users using automotive radars,” _IEEE Access_ , vol. 9, pp. 68 758–68 777, 2021.
  * [19] N. Lyons, A. Santra, and A. Pandey, “Improved deep representation learning for human activity recognition using imu sensors,” in _2021 20th IEEE International Conference on Machine Learning and Applications (ICMLA)_. IEEE, 2021, pp. 326–332.
  * [20] A. Kendall and Y. Gal, “What uncertainties do we need in bayesian deep learning for computer vision?” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [21] N. Agarwal and S. Das, “Interpretable machine learning tools: A survey,” in _2020 IEEE Symposium Series on Computational Intelligence (SSCI)_. IEEE, 2020, pp. 1528–1534.
  * [22] S. Lundberg, “Shap (shapley additive explanations),” 2020.
  * [23] B. Rozemberczki, L. Watson, P. Bayer, H.-T. Yang, O. Kiss, S. Nilsson, and R. Sarkar, “The shapley value in machine learning,” _arXiv preprint arXiv:2202.05594_ , 2022.
  * [24] S. Singha, A. S. Singh, S. Prasad, and A. Alam, “A study on power optimization techniques in psoc,” in _2019 International Conference on Communication and Signal Processing (ICCSP)_. IEEE, 2019, pp. 0225–0229.
  * [25] S. Patro and K. K. Sahu, “Normalization: A preprocessing stage,” _arXiv preprint arXiv:1503.06462_ , 2015.
  * [26] S. Soi, “Design and characterization of a 3rd order low-pass butterworth filter,” _International Journal of Science and Research (IJSR)_ , vol. 3, no. 7, pp. 1–3, 2014.


∎

11institutetext:  H. Gao 22institutetext:  Dept. Mat. Sci. & Eng., Saarland
University, Campus C6 3, 66123 Saarbrücken, Germany 33institutetext: M. H.
Müser 44institutetext:  Dept. Mat. Sci. & Eng., Saarland University, Campus C6
3, 66123 Saarbrücken, Germany 44email<EMAIL_ADDRESS>
# On the shear-thinning of alkanes

Hongyu Gao and Martin H. Müser

(Received: date / Accepted: date)

###### Abstract

The approximate power-law dependence of the apparent viscosity of liquids on
shear rate is often argued to arise from a distribution of energy barriers.
However, recent work on the Prandtl model, which consists of a point mass
being dragged by a damped, harmonic spring past a sinusoidal potential,
revealed a similar dependence of the friction on velocity as that of many
liquids. Here, we demonstrate that this correlation is not only qualitative
but can also be made quantitative over a broad temperature range using merely
three dimensionless parameters, at least for alkanes, in particular
hexadecane, at elevated pressure $p$. These and other observations made on our
all-atom alkane simulations at elevated pressure point to the existence of an
elementary instability causing shear thinning. In addition, the equilibrium
viscosity shows power law dependence on $p$ near the cavitation pressure but
an exponential dependence at large $p$, while the additional parameter(s) in
the Carreau-Yasuda equation compared to other rheological models turn out
justifiable.

##  1 Introduction

Research on how liquids oppose shear as a function of shear stress or shear
rate, temperature, and also pressure supposedly started almost exactly a
century ago Barnes1997JNNFM , when Schalek and Szegvari Schalek1923KZ
discovered that the apparent viscosity $\eta$ of iron-oxide colloidal
suspension, defined as the ratio of shear stress $\tau$ and shear rate
$\dot{\gamma}$, became continuously smaller with increasing shear rate.
Ostwald Ostwald1925KZ  soon found that this shear thinning, termed
“Strukturviskosität” at the time (structure or intrinsic viscosity in
English), can be approximated rather well with an
$\eta~{}\propto\dot{\gamma}^{n-1}$ power law, where $n$ is called the shear-
thinning exponent. It is now well established that any liquid, even the
prototypical Newtonian fluid water when supercooled deAlmeidaRibeiro2020PRR ,
has a shear-rate dependent viscosity if the shear rate is only sufficiently
high so that labeling fluids as (approximately) Newtonian or not is not a
matter of chemistry but mostly of shear rate. Liquids obeying Ostwald’s
empirical law with $n>1$, $n=1$, and $n<1$ can be labeled shear-thickening,
Newtonian, and shear-thinning, respectively. Under certain situations, $n<0$
might be observable, in which case the friction force decreases with velocity,
e.g., in computer simulations of poly-alpha-olefins and
2,4-dicyclohexyl-2-methylpentane at shear rated exceeding 1 GHz Ewen2019PCCP .
However, macroscopic flow could not be Couette like in such situations,
because unavoidable fluctuations of the laminar flow profile would instantly
grow, thereby resulting in a dynamical instability producing narrow zones of
large slip.

The equilibrium viscosity $\eta_{0}(T,p)$ generally increases with pressure
$p$ and decreases with temperature $T$ Bair2003TT ; Roland2006JCP . One
central reason for this behavior certainly is that elevating pressure makes it
more difficult for molecules to move past each other, mainly because of
increased steric repulsion, while higher temperature assists molecules to
overcome energy barriers, whereby shear flow is facilitated. The functional
form of $\eta(\dot{\gamma})$, e.g., the exponent $n$ in (post-) Ostwald shear-
thinning laws, evolves rather smoothly with pressure and temperature with
potential exceptions occurring in the vicinity of a liquid-liquid phase
transformation. This has enticed researchers to estimate effective viscosities
at temperatures, pressures, and shear rates that are difficult to access
experimentally by using time-temperature and pressure superposition principles
Bair2002PRL , similar to those used for the construction of master curves
describing the linear viscoelasticity of rubbers or glass-forming liquids
Li2000MSE .

It remains actively debated what shear-thinning model describe real-laboratory
or in-silico experimental data the best Ewen2021TL . For example, Bair et al.
Bair2015TL  objected to the validity of the Eyring model Eyring2004JCP  to
describe rheological data, which had been advocated by Spikes and Jie
Spikes2014TL ; Spikes2015TL . In a possible retaliation, Spikes Spikes2017TL
dismissed work by Voeltzel, Vergne et al. Voeltzel2016TL ; Voeltzel2017TL  on
the rheology ionic liquids with variable shear thinning exponent, by favoring
the Eyring model over the originally used Carreau model Carreau1972TSR ].
(These and related models are introduced in Sect. 2). Also the question what
(simple) theory provides the best microscopic explanation for shear thinning
in fluids remains discussed, e.g., whether Eyring’s Eyring2004JCP  or
Prandtl’s Prandtl1928ZAMM  approach is more appropriate Spikes2015TL-R ,
although it seems often unappreciated that Prandtl had already worked out
Eyring theory as a limiting case of his model and that Eyring admitted his
results to be similar to Prandtl’s and related prior work.

Another debate evolves around the question what time-temperature, or even
time-temperature-pressure superposition principle is best suited for
lubricants. For example, Bair Bair2017PNAS  objected to the cross-over from a
non-Arrhenius to an Arrhenius dependence of the viscosity upon cooling that
Jadhao and Robbins had claimed to have identified for squalane Jadhao2017PNAS
; Jadhao2017PNAS-R , though it seems that the transition from non-Arrhenius to
Arrhenius dependencies of viscosity and diffusivity, which in some cases is
called a fragile-to-strong transition, is a generic phenomenon of glass-
forming melts Lucas2021JNCS .

In this work, we revisit the dependence of the viscosity of simple liquids
depends on shear-rate, temperature, pressure and molecular weight. While
numerous simulation studies Bair2002PRL ; Jadhao2017PNAS ; Daivis1994JCP ;
Sivebaek2012PRL ; Lemarchand2015JCP  have computed $\eta(\dot{\gamma})$
relations for realistic lubricant models in the past, often with the attempt
to correlate molecular rearrangement with shear thinning Lemarchand2015JCP ,
see Ref. Ewen2021TL  for a recent review, we focus here on the question if the
similarity of shear-thinning in realistic model liquids and that observed in
the Prandtl model Prandtl1928ZAMM  are merely qualitative in nature or if they
could be quantitative. To this end, we generate a reference relation of
$\eta(\dot{\gamma})$ for hexadecane at $T=500$ K and a hydrostatic pressure of
$p=300$ MPa and identify a Prandtl model that reproduces the
$\eta(\dot{\gamma})$ relationship satisfactorily for that reference. Next,
either temperature or pressure or the molecular weight of the alkane is varied
and the changes in equilibrium viscosity studied and discussed in the
framework of the Prandtl model, which had been introduced as a simple model
for dislocation motion and argued to also explain the velocity dependence of
solid friction or the shear thinning of fluids, around similar times when
shear thinning was reported for the first time Prandtl1928ZAMM ; Popov2012ZAMM
.

For our investigation, in which we mean to understand trends rather than to
produce accurate numbers for specific molecules, we see computer simulation as
the ideal tool, despite essentially unavoidable flaws in interaction
potentials and the difficulty to predict reliable viscosity data for shear
rates below 107 Hz. However, the advantage of simulations is that high shear
rates can be reached while conditions such as shear rate, pressure, and
temperature can be set to a target precision. Moreover, undesired effects like
frictional heating can be suppressed or at least be quantified. Last but not
least, the pressure dependence of the viscosity can be investigated even under
a well-defined isotropic, tensile stress.

We continue this paper with a background section on rheological models in
Sect. 2. Theory and numerical methods are introduced in Sect. 3. Results are
presented in Sect. 4, while conclusions are drawn in Sect. 5.

##  2 Background on rheological models

The rheology of a broad class of different liquids including standard base
oils, water, or blood is captured quite well by the Carreau-Yasuda (CY)
equation,

$\eta(\dot{\gamma})=\eta_{\infty}+(\eta_{0}-\eta_{\infty})\left\\{1+\left(\dot{\gamma}/\dot{\gamma}_{0}\right)^{a}\right\\}^{(n-1)/a},$
(1)

where $\eta_{0}$ is the equilibrium viscosity, $\dot{\gamma}_{0}$ is a
reference shear rate, near which shear thinning starts or has started to set
in, while $\eta_{\infty}$ is a parameter discussed in detail in the next
paragraph. Two other popular models arise as limiting cases of Carreau-Yasuda,
namely the Carreau model Carreau1972TSR  for $a=2$ and the Cross model
Crook1963MPS  for $a=1-n$.

An arguably tricky term in Eq. (1) is the large-shear-rate limit of the
viscosity, $\eta_{\infty}$, even if it is often useful to introduce it in
practice so that fitting functions extend their validity to large shear rates.
However, at the point at which $\eta(\dot{\gamma})$ appears to level off, the
way in which stress and temperature are controlled is likely to matter. For
example, we find that fitting viscosity data when reported at fixed pressure
does not require $\eta_{\infty}$ to be introduced, while it does at constant
density. These observations are in line with classical results by Daivis and
Evans Daivis1994JCP , who found decane to show a finite value of
$\eta_{\infty}$ or even to shear thicken at constant-density simulations but
not at constant pressure. For these reasons, $\eta_{\infty}$ will be ignored
in the following, all the more extreme shear-thinning at very large
$\dot{\gamma}$ is readily built in and thus rationalized in the Prandtl model,
see Sect. 3.1.

Coming back to the discussion of rheological models: Carreau Carreau1972TSR
had originally attempted to simultaneously address the linear-response
viscoelasticity and the out-of-equilibrium shear thinning by assuming the
rates of creation and annihilation of molecular junctions to be functions of
what now would be called the equivalent von Mises strain-rate. Carreau
suggested two ad-hoc functional dependencies for those functions, one of which
he dismissed as unrealistic, which, ironically is the one leading to the
equation now carrying his name, i.e., Eq. (1) for $a=2$.

A beneficial aspect of the Carreau model is that the leading-order corrections
to the linear-response viscosity are, as they should be, of order
$\dot{\gamma}^{2}$ Muser2020L . However, it seems that expanding the Carreau
equation into (even) powers of the shear rate does not accurately reflect the
initial but only the later deviations of the viscosity from the equilibrium
viscosity. We base this conclusion on the observation that almost any high
resolution data that we are aware of, be it experimental or simulations, lies
below the fits at the point where $\eta$ has decreased to say 80% its
equilibrium value.

Another shear-rate dependence of the viscosity, which can be motivated from
Eyring theory, is

$\frac{\eta}{\eta_{0}}=\frac{\dot{\gamma}_{0}}{\dot{\gamma}}\mathrm{asinh}\left(\frac{\dot{\gamma}}{\dot{\gamma}_{0}}\right),$
(2)

where the ideal value for $\dot{\gamma}_{0}$ may differ from that of another
rheological model. This relation makes the shear stress increase
logarithmically with $\dot{\gamma}$ at large $\dot{\gamma}$ while converging
to the equilibrium viscosity $\eta_{0}$ at $\dot{\gamma}\to 0$, as is typical
for boundary lubrication. This is arguably the main reason for the use of this
and related relations in state-rate models of frictional contacts Ruina1983JGR
. However, what appears to be often overlooked is that Eq. (2) can be seen as
a limiting case of Eq. (1), because an increase of shear stress or friction
with a power-law $\dot{\gamma}^{\varepsilon}$ with $0<\varepsilon\ll 1$
behaves essentially logarithmically. Thus, as long as $0<n\ll 1$,
discrepancies between Eyring and the other rheological models can only be
minor, in which case the only advantage of Eyring would be to depend on one
less adjustable parameter.

When evaluating the functional form of a given model, mostly the viscosity and
sometimes the shear stress or friction are often represented in a double
logarithmic plot as a function of shear rate or velocity. This way differences
between different models near the cross-over from Stokes friction to “Ostwald
scaling” have a low optical resolution. It can be improved, in particular for
shear-thinning exponents near “Eyring scaling”, i.e., $0<n\ll 1$, when
multiplying the viscosity with the square-root of the shear rate Muser2011PRB
; Muser2020L . An example of such a what-we-take-the-liberty-to-call-for-
reasons-of-simplicity-and-not-vanity Müser plot is shown in Fig. 1. It
contrasts the regular $\eta(\dot{\gamma})$ representation for the exponent
$n=0.2$ with a Müser plot. In this example, $n$ was chosen roughly half way
between the values identified in this study and $n=1-\varepsilon$,
$0<\varepsilon\ll 1$, at which point the algebraic dependence of shear stress
on shear rate is logarithmic like.

Figure 1: Shear-rate dependence of the viscosity $\eta$ as produced by various
model, (a) in a regular representation and (b) in an alternative fasion. All
model functions assume the same linear-response viscosity $\eta_{0}$ and a
similar $\eta(\dot{\gamma}=10^{3})$. Carreau-Yasuda (blue line), Carreau
(red), and Cross (green) employ the same value for $n$ and $\dot{\gamma}_{0}$.

To quantify the goodness of fits to molecular dynamics (MD) results, we
defined a $\chi^{2}$-error as a logarithmic deviation of the viscosity
function from the reference data, i.e.,

$\chi^{2}=\frac{1}{N}\sum_{n=1}^{N}\left\\{\ln\frac{\eta_{\textrm{model}}(\dot{\gamma}_{n})}{\eta_{\textrm{MD}}(\dot{\gamma}_{n})}\right\\}^{2},$
(3)

where $N$ enumerates the data points. This way, errors do not depend on
whether stress or viscosity are targeted. To not bias the fits to large or low
shear rates, we furthermore used $\dot{\gamma}_{n+1}/\dot{\gamma}_{n}$
constant.

##  3 Theory and methods

###  3.1 Inertia-free Prandtl model

In the Prandtl model, a mass point is pulled with a spring over a corrugated
substrate. In addition to the conservative forces, Stokes’ friction terms
arise, which are argued to reflect lattice vibrations in the substrate.
However, this explicit damping can also be assumed to occur in the pulling
spring, in which case the mass point’s equation of motion in the presence of
thermal noise $\Gamma_{\textrm{th}}(t)$ reads

$\gamma_{0}(\dot{x}-v_{0})=qV_{0}\sin(qx)+k(x-v_{0}t)+\Gamma_{\textrm{th}}(t)$
(4)

Here $x$ is the position of the “Prandtl atom”, $\dot{x}$ its velocity,
$v_{0}$ the velocity of the pulling spring, $\gamma_{0}$ is a damping term of
unit kg/s, $V$ the amplitude of the sinusoidal potential through which the
atom is dragged, $2\pi/q$ is the period of the potential, $k$ the stiffness of
the pulling spring, and $\Gamma_{\textrm{G}}$ is thermal noise with the
expectation values

$\displaystyle\langle\Gamma_{\textrm{th}}(t)\rangle$ $\displaystyle=$
$\displaystyle 0$ (5)
$\displaystyle\langle\Gamma_{\textrm{th}}(t)\Gamma_{\textrm{th}}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\gamma_{0}k_{\textrm{B}}T\delta(t-t^{\prime}),$ (6)

where $k_{\textrm{B}}T$ the thermal energy and $\delta(...)$ is the Dirac
delta function. We have the liberty to use a unit system of our choice, i.e.,
one in which $q$, $V_{0}$, and $\gamma$ define the units used for inverse
length, energy, and weight per time. When doing so, it becomes clear that this
form of the Prandtl model has two parameters defining the shape of a $F(v)$
relation, since only the reduced stiffness $\tilde{k}\equiv k/(q^{2}V)$ and
the reduced temperature $\tilde{T}\equiv k_{\textrm{B}}T/V$ remain, in which
case the dynamical equation to be solved, introducing
$\langle\tilde{\Gamma}_{\textrm{th}}(t)\tilde{\Gamma}_{\textrm{th}}(t^{\prime})\rangle=\delta(t-t^{\prime})$,
simplifies to

$\dot{x}-v_{0}=\sin
x+\tilde{k}(x-v_{0}t)+\sqrt{2\tilde{T}}\tilde{\Gamma}_{\textrm{th}}(t),$ (7)

In this model, only two dimensionless numbers, $\tilde{k}$ and $\tilde{T}$,
determine the shape of a $\eta(\dot{\gamma})$ relation, e.g., the exponents
$n$ and $a$ that would be obtained when fitting simulation data to a Carreau-
Yasuda relation. Absolute values, such as the equilibrium viscosity $\eta_{0}$
or the cross-over shear rate $\dot{\gamma}_{0}$ are then determined by the
model parameters defining the units in the Prandtl model. For example, given a
value of $\tilde{k}$ and assuming $\tilde{k}$ and $V_{0}$ to be temperature-
independent, the small-temperature dependency of the equilibrium-damping
constant of the Prandtl model, which is proportional to the viscosity, reveals
Arrhenius type dependence Muser2020L  with the energy activation barrier
$\Delta E$ whose construction is shown in Fig. 2. Assuming the temperature to
be below the liquid’s fragile-to-strong transition temperature and $\Delta E$
to be know, the unit for $V_{0}$ can then be fixed.

Figure 2: Substrate and total potential energy for $\tilde{k}=0.7$, which is
used at the reference state (300 MPa, 500 K, which corresponds to
$\tilde{t}=0.076$). The energy barrier $\Delta E=0.14$ is defined at the
moment of time when the position of the driving spring coincides with the
maximum of the substrate potential.

When varying temperature of a large range, it seems clear that both $V_{0}$
and $k$ depend on temperature. Specifically, at large $T$, the all-atom liquid
expands, which will reduce $V_{0}$ supposedly more than $k$. To reflect such
effects of thermal dilation in the model, we allowed $\tilde{k}$ to obey

$\tilde{k}=\tilde{k}_{\infty}\exp(-T_{k}/T),$ (8)

where $T_{k}$ is an empirical parameter, which would need to get readjusted as
a function of pressure and the degree of polymerization. The exponential
dependence is motivated by the Boltzmann distribution without any further
justification. However, we note that the value used for $T_{k}$ was only
$T=90.5$ K, is clearly less than the smallest investigated temperature of 300
K.

Before discussing some important asymptotic limits of the Prandtl model, we
wish to justify why we damp the spring motion rather than the velocity of an
atom or rather a material point relative to the substrate. First, the kinetic
energy of monolayers sliding across a substrate was shown to be related to the
(anharmonic) coupling of vibrations in the monolayer Smith1996PRB  and we see
a Prandtl spring to reflect the coupling of a central atom to its in-plane
neighbors of a lamella during laminar flow. Second, damping relative to the
substrate would yield a high-velocity damping, which would produce a finite
$\eta_{\infty}$. However, as discussed in Sect. 2, we see no evidence for the
existence of $\eta_{\infty}$, when evaluating $\eta(\dot{\gamma})$ at a
constant pressure.

A few properties of the Prandtl model are worth reciting, most of which can be
found in Fisher1985PRB . First, instabilities and thus (athermal) static
friction and Coulomb-like friction at low temperatures can only occur when
$\tilde{k}<1$. Second, the athermal low-velocity limit of the friction-
velocity relation reads $F_{\textrm{k}}-F_{\textrm{s}}\propto v^{2/3}$. Thus,
for $\tilde{k}$ being in the immediate vicinity of unity, a shear-thinning
exponent of $n=2/3$ is obtained. This value thereby constitutes an upper bound
for what shear-thinning exponent $n$ can be reproduced with the Prandtl model.
Third, in the limit $0<\tilde{k}\ll 1$, the material point is effectively
moved with a constant force, which is the limit at which the Prandtl model
corresponds to situation described by himself and later again by Eyring. In
that limit, the effective damping at small but finite temperature can be
described with Eq. (2). Thus, the Prandtl model can reproduce any shear-
thinning exponent between $0<n\leq 2/3$, where a small exponent would imply
small $\tilde{k}$ and/or finite temperature, while $n=2/3$ would suggest
thermal fluctuations to be minor and/or $\tilde{k}$ to be relatively close to
unity. In fact, the Prandtl model reproduces the general trend of real
lubricants Ewen2021TL  that $n$ decreases with decreasing temperature
Muser2020L .

The fourth asymptotic limit to be discussed is the high-velocity limit. In
this case, the atom moves at constant velocity of $v_{0}$ plus small
fluctuations, which can be described using perturbation theory, The first-
order perturbation is a time-dependent force $qV\sin(\omega t)$ with
$\omega=qv_{0}$. For $\omega\to\infty$, the response function
$\tilde{x}_{1}(\omega)~{}\propto 1/\omega$ so that the dissipated power
$P\propto|\omega\tilde{x}_{1}|^{2}$ is independent of $\omega$, which in turn
means that the friction force is proportional to $1/\omega$ or $1/v$, which,
when applied to laminar flow would translate to $n=-1$. Couette flow would no
longer be possible, as discussed in the introduction. Adding inertia to the
model would further impede Couette flow at large $\dot{\gamma}$ as
$\tilde{x}_{1}$ would then be proportional to $1/\omega^{2}$ at large
$\omega$, which decreases $n$ to $n=-3$. Irrespective of the choice of the
mass, extreme shear thinning at very high shear rate occurs in the Prandtl
model, when the spring head moves so fast that the energy released during an
instability can no longer be transferred into heat before the atom undergoes
the next instability. We would argue that similar statements can be made for
any other systems and degrees of freedom showing decreasing friction with
increasing velocity.

###  3.2 Molecular-dynamics alkane simulations

Non-equilibrium molecular dynamics (NEMD) simulations are carried out on a
bulk liquid consisting of approximately 5,000 atoms using the open-source
code–LAMMPS Plimpton1995JCP . Interatomic interactions are modeled with the
all-atom L-OPLS-AA force field Price2001JCC ; Siu2012JCTC , which describes
bond-stretching, bending, and torsional in addition to non-bonded
interactions. The shear viscosity is calculated through Couette-flow
simulations in which shear flow is implemented via Lees-Edwards Lees1972JPC
periodic boundary conditions and the SLLOD algorithm Todd2017 . The system
temperature, $320~{}\textrm{K}\leq T\leq 550$ K was maintained with a Nosé-
Hoover thermostat. To this end, the LAMMPS command temp/deform was used that
removes the effect of a mean velocity profile from the determination of
kinetic energy and thus temperature.

Results are reported for a fixed normal pressure, which was achieved as
follows: In a first set of short simulations at fixed volumes and shear rates,
an equation of state for the give shear rate was ascertained and the volume
determined by interpolation (never extrapolation!) at which the target
pressure was expected to occur. This volume was then used in a longer run
during which the effective viscosity was calculated. The average normal stress
in those long runs always turned out within 1% at finite stress and $\pm 1$
MPa at zero stress, where small relative errors are somewhat difficult to
achieve. In the case of moderate, mean tensile stresses, the same average
value $\sigma_{33}$ value can be reached with two different densities.

##  4 Results

###  4.1 Temperature dependence

We first investigate the dependence of the viscosity of hexadecane as a
function of temperature at a pressure of $p_{33}=300$ MPa, where
$p_{33}\equiv-\sigma_{33}$. The latter value was chosen such that it is high
enough to almost represent boundary-lubrication conditions but small enough so
that the equilibrium viscosity could be calculated down to a temperature of
$T=320$ K, which would be at the cold end of engine oil after warming up. Fig.
3 shows the raw data obtained from MD simulations and fits based on the
Carreau-Yasuda equation.

Figure 3: Shear viscosity $\eta$ of $n$-hexadecane at various temperatures
(320-550 K) and a constant normal pressure of 300 MPa. Symbols show NEMD data
and dashed lines fits to the data using the Carreau-Yasuda equation.

The range of shear rates presented in Fig. 3 and others presented in this work
is certainly beyond any practical relevance, which should be roughly
$\dot{\gamma}=1$ GHz, obtained when two surfaces moving at a relative velocity
of 2 m/s have a local separation of 2 nm. However, owing to time-temperature
and time-pressure superposition principles, similar rheology (except for
prefactors) can occur at “reasonable” shear rates when the temperature is
lowered and/or the pressure increased. Nonetheless, using shear-rates clearly
above $10^{12}$ Hz would not be meaningful, because these rates exceed local
librational frequencies.

Although it has been a matter of ongoing debates what functional form
describes data as that presented in Fig. 3 best, we feel that a “Müser plot”
allows differences between fit and data to be much better resolved than
traditional representations of $\eta(\dot{\gamma})$. For the data in question,
our results for viscosity spanning two decades in a regular
$\eta(\dot{\gamma})$ representation as in Fig. 3 vary merely by a factor of
about three in the cross-over domain, as can be seen in Fig. 4. This improves
the resolution of the graph effectively by a factor of approximately 60,
whereby both small deficiencies of rheological model or stochastic errors are
relentlessly revealed.

Figure 4:  A comparison of phenomenological models fitted to a representative
case of $n$-hexadecane at $T$ = 350, 400, and 500 K and $p_{33}$ = 300 MPa
with normalized
$\hat{\dot{\gamma}}=\dot{\gamma}/\dot{\gamma}_{0}^{\textrm{CY}}$ and
$\hat{\eta}=\eta/\eta_{0}^{\textrm{CY}}$, where the upper index CY indicates
that the Carreau-Yasuda values at the given values were taken as units for
viscosity and frequency. Results from $T=400$ K and 500 K are shifted in $y$
direction by 0.5 and 1.0, respectively.

Fig. 4 compares three rheological models for the viscosity of hexadecane at
three different temperatures. It can be seen by the naked eye that the Eyring
model produces the data in the least satisfactory manner. Carreau is a
substantial improvement to Eyring, however, Carreau overestimates the
viscosity at shear rates slightly lower than the maximum. We invite the reader
to check for this deviation to occur in other published data and invite him or
her to remember our assessment that this short-coming of Carreau shows up in
the early shear thinning regime the more clearly the less noisy the real-
laboratory or in-silico data.

Relative standard-deviations converted in percentages, $\chi^{\%}\equiv
100\cdot\chi$, are listed in Table 1 for all rheological models introduced in
this paper so far and all MD-based $\eta(\dot{\gamma})$ reference data. The
latter includes data obtained by varying the pressure or the molecular weight
rather than the temperature from their default values.

3 $P$ |  $T$/K |  $p_{3}$/MPa |  $\eta_{0}$/cP |  $\dot{\gamma}_{0}$/GHz |  $n$ |  $a$ |  $\chi_{\textrm{CY}}^{\%}$ |  $\chi_{\textrm{C}}^{\%}$ |  $\chi_{\textrm{Cr}}^{\%}$ |  $\chi_{\textrm{E}}^{\%}$
---|---|---|---|---|---|---|---|---|---|---
16 |  320 | 300 | 26.76 | 0.71 |  0.310 |  1.28 |  4.12 |  5.45 |  6.46 |  17.77
|  350 | 13.25 | 1.67 |  0.327 |  1.17 |  1.93 |  4.11 |  4.56 |  15.78
|  400 | 5.24 | 4.50 |  0.376 |  1.42 |  1.86 |  2.79 |  5.41 |  13.73
|  450 | 2.96 | 9.65 |  0.391 |  1.34 |  1.25 |  2.67 |  4.53 |  11.53
|  500 | 2.00 | 16.37 |  0.406 |  1.25 |  1.65 |  3.00 |  3.71 |  10.64
|  550 | 1.46 | 22.35 |  0.441 |  1.36 |  0.97 |  2.06 |  3.51 |  9.96
16 |  500 | 100 | 0.66 | 59.25 |  0.300 |  1.16 |  1.61 |  3.29 |  2.41 |  4.66
|  | 200 | 1.20 | 26.37 |  0.404 |  1.35 |  1.67 |  2.58 |  3.69 |  8.03
|  | 300 | 2.00 | 16.37 |  0.406 |  1.25 |  1.65 |  3.00 |  3.71 |  10.64
|  | 400 | 2.97 | 9.66 |  0.425 |  1.45 |  1.74 |  2.44 |  5.00 |  12.63
|  | 500 | 4.36 | 7.61 |  0.398 |  1.26 |  1.84 |  3.16 |  4.38 |  13.47
|  | 600 | 5.90 | 5.31 |  0.404 |  1.45 |  1.07 |  2.14 |  5.45 |  14.46
16 |  500 | 300 | 2.00 | 16.37 |  0.406 |  1.25 |  1.65 |  3.00 |  3.71 |  10.64
20 |  | 2.81 | 8.54 |  0.421 |  1.32 |  1.77 |  2.88 |  4.61 |  13.62
24 |  | 3.83 | 5.26 |  0.418 |  1.19 |  1.97 |  3.53 |  4.51 |  16.26
3 |  |  |  |  |  |  |  |  |  |
Table 1: Results for Carreau-Yasuda (CY) fits to MD data. $P$ stands for the
number of backbone carbons. Relative standard deviations are also included for
fits to Carreau (C), Cross (Cr), and Eyring (E).

Since Carreau (Ca) and Cross (Cr) are exact limits of Carreau-Yasuda (CY),
they cannot outperform CY. It is similarly clear that Eyring (E) should be the
least accurate of the four models, because it arises as an approximate
asymptotic limit from the other models. However, the interesting question is
to what degree increasing the number of fit parameters from 2 (E) to 3 (Ca,
Cr) to 4 (CY) improves the goodness of the fits. Our personal view is that the
dividing line between justified and unjustified additional parameter is a
reduction of $\chi$ by a factor that should be clearly smaller than $k/(k+1)$
when $k$ is the number of originally adjustable coefficients.

Table 1 reveals quite clear trends. Errors associated with Eyring exceed 10%
in all but two out of 12 sets of simulations. The Cross model reduces the
errors compared to Eyring by a factor always exceeding 1.9 and in the clear
majority of cases by a factor greater than 2.5. Thus, the additional
adjustable parameter seems justified. Carreau further reduces the error
compared to Cross quite substantially, except at the highest temperature and
the lowest pressure reported in Table 1. Carreau-Yasuda is again a clear
improvement on Carreau. Improvements scatter between reduction factors of 1.3
and 2.2, however, the median improvement in the reported data lies at about
1.8. For this reason, we judge the extra parameter $a$ in Carreau-Yasuda as
significant, in particular as the Carreau-Yasuda errors are often only
marginally larger than the statistical uncertainty of the produced data.
Unfortunately, it is beyond our current resources to substantially reduce the
thermal/statistical errors, because halving statistical errors implies
quadrupling the numerical effort.

While all investigated model functions, even that based on Eyring, matches the
MD data in the representation of Fig. 4 substantially better than an
“arbitrary fitting function”—such as a skewed Gaussian, which, like Carreau-
Yasuda, has four adjustable parameters—Cross, Carreau, and Carreau-Yasuda are
purely heuristic functions with built-in linear response at small shear rates
and power-law dependence at large shear rates. However, unlike Eyring and even
more so Prandtl, they lack a clear microscopic justification. The Prandtl
model is a physics-based model, in which linear response at small velocities
arises naturally while quasi-power-law behavior occurs at large shear rates
Popov2012ZAMM ; Muser2020L .

The just-mentioned feature of the Prandtl model motivated us to explore to
what extent the two dimensionless parameters of the Prandtl model can be tuned
to match experimental or realistic, in-silico data. The interesting aspect of
such an analysis is that the model parameters can be fixed at one temperature
and changes in $\eta(\dot{\gamma})$ be estimated without requiring adjustment
of the parameters beyond including a smooth temperature dependence of
$\tilde{k}$.

Fig. 5 reveals that three dimensionless parameters and thus a total of five
parameters is sufficient to reproduce the full-scale NEMD results at different
temperatures. Interestingly, it also reproduces the observation that the
exponent $n$ decreases thereby approaching Eyring-like behavior with
decreasing temperature. Specifically, the full simulations show $n=0.391$ at
$T=450$ K and $n=0.327$ at $T=350$ K, while the Prandtl model finds $n=0.366$
and $0.321$ at those two temperatures, respectively. Thus, the Prandtl model
does not perfectly reproduce the reduction in the shear-thinning exponent but
supposedly much better than expected from a simple model.

Figure 5:  Comparison of results obtained from the Prandtl model to Carreau-
Yasuda fits to the hexadecane MD results at $p_{33}=300$ MPa (dashed curves).
The Prandtl model was chosen with $\tilde{k}=0.84\cdot(1-90.5~{}K/T)$ and
$\tilde{t}=1.52\cdot 10^{-4}~{}T$/K.

Fig. 6 further reveals the large degree of correlation between the results of
all-atom simulations and those of the Prandtl model. Both data sets have the
same inverse correlation of $\eta_{0}$ and $\dot{\gamma}_{0}$, as would be
expected from the equilibrium dynamics of an activated process.

Figure 6:  Relation between equilibrium viscosity $\eta_{0}$ and cross-over
shear rate $\dot{\gamma}_{0}$ as obtained in full NEMD simulations and the
Prandtl model. The dashed line reflects a slightly sublinear, inverse power
law, i.e., $\eta_{0}\propto 1/\dot{\gamma}^{0.84}$, which is drawn to guide
the eye.

The final analysis in this section is concerned with the temperature
dependence of hexadecane’s viscosity at the relatively large, constant
compressive stress of $p_{33}=300$ MPa. Fig. 7 hints to a cross-over from a
non-Arrhenius to an Arrhenius-like dependence in the all-atom MD data.
Unfortunately, some doubt on the robustness of the Arrhenius-type dependence
at small temperature remains, as we have not yet managed to compute an
equilibrium viscosity below 320 K within our statistical error margin of about
1% from the full MD simulations. However, the Prandtl model, which accounts
for the cross-over reasonably well, while clearly having a constant energy
barrier supports an extension of an Arrhenius dependence to smaller
temperature. The model thereby also reveals that the deviation from Arrhenius
at large $T$ might not be caused by temperature-dependent energy barriers but
could be owing to the ratio $\Delta E/k_{\textrm{B}}T$ not being a large
number, at which point the exponential factor $\exp(\Delta E/k_{\textrm{B}}T)$
no longer dominates other temperature dependencies.

Figure 7:  Newtonian viscosity $\eta_{0}$ from all-atom simulations (full, red
circles) and Prandtl model (open, black circles) and inverse diffusion
coefficient $1/D$ of hexadecane (blue crosses) as a function of inverse
temperature, for $n$-hexadecane, at 300 MPa. Both ordinates, i.e., $\eta_{0}$
on the left and $1/D$ on the righ, span the same factor of 100.

To further illuminate the reason for the FST, we computed the diffusion
coefficient $D$ of carbon atoms from their time-dependent root-mean-square
displacement (see also Sect. 4.4) to explore if the temperature dependence of
the diffusion coefficient starts to decouple from that of the viscosity, as it
happens at the FST of glass-forming melts. Fig. 7 provides some evidence for
this to happen, which, however is relatively weak, because accurate diffusion
coefficients at small temperatures are difficult to determine. Thus, more
research, which is currently outside the scope of our computational
feasibility, is required to see if (a) the viscosity is really Arrhenius-like
at small temperature and (b) if the cross-over can already be rationalized
within the Prandtl model or is in need of a more elaborate explanation like
the FST. The difficulty is that statistical error bars are small at relatively
small times, at which point the functional form used to fit the time-dependent
root-mean-square displacement affect the value for $D$, while stochastic
errors are large at large times. Thus, values for $1/D$ reported at the lowest
temperatures must be taken with a grain of salt and might increase by roughly
20% with different fitting procedures than those that we assumed.

###  4.2 Pressure dependence

Besides the temperature dependence of the viscosity $\eta$ at fixed pressure
$p_{33}$, we also analyzed $\eta(p)$ at fixed temperature. This time, however,
we present the results in the form of shear stress $\tau$ as a function of
shear rate $\dot{\gamma}$, see the main panel in Fig. 8. Results are also
reported for simulations at a negative normal pressure. They reflect meta-
stable conditions, because a thermally or shear-induced void in the liquid
would unavoidably grow indefinitely under a tensile load once the void
surpasses a critical nucleation size. However, including tensile stresses into
the analysis allows us to explore the pressure-dependence on viscosity at
values of $p_{33}$, where this dependence is particularly steep.

Figure 8:  Shear stress $\tau$ of $n$-hexadecane at various pressures (0-600
MPa) and a constant temperature of 500 K. Symbols show NEMD data and dashed
lines fits to it using the Carreau-Yasuda equation. At high shear rates, the
system can be either in the liquid or in the gas phase. Open symbols refer to
the gas phase as well as data for the liquid phase, if some configurations
revealed shear localization.

Analysis of the equilibrium viscosity, see Fig. 9, reveals that $\eta_{0}$
might disappear as a power law at a negative pressure of $p^{*}=-10.5$ MPa.
This value is close to the pressure at which the volume of a meta-stable
liquid, having been “stress-quenched” to a negative pressure, is expected to
diverge. However, at large pressure, $\eta(p)$ crosses over to an exponential
dependence of pressure. A simple functional form reflecting this behavior is

$\eta(p)=\eta_{1}\left(p/p^{*}-1\right)^{\mu}+\eta_{2}\left(e^{\beta
pV^{*}}-e^{\beta p^{*}V^{*}}\right),$ (9)

where the parameters $\eta_{1}$, $\eta_{2}$, $\mu$, and $V^{*}$ are
phenomenological adjustable parameter and $\beta=1/k_{\textrm{B}}T$ $V^{*}$
can be interpreted as $V^{*}\approx\partial\Delta E/\partial p$ at
intermediate pressures and be called “free volume” or Peter, both names
probably being similarly meaningful, also because $V^{*}$ can take a different
value in non-equilibrium situations when replacing the hydrostatic pressure
$p$ in Eq. (9) with the normal pressure $p_{33}$. The inset of Fig. 9 shows
the viscosity as a function of volume. No extended domain of a linear
$\ln\eta_{0}$ versus volume can be detected.

Figure 9: Main graph: Newtonian viscosity $\eta_{0}$ as a function of normal
pressure $p$, $n$-hexadecane, 500 K. Symbols show numerical results, while the
line is a fit to Eq. (9). Inset: $\eta_{0}$ as a function of the volume $V$
per atom. The line is a fourth-order polynomial of $\ln\eta_{0}$ in $V$.

A few more words on the decrease of viscosity at $p\gtrsim p^{*}$ may be in
place. The line represented in Fig. 9 is based on an exponent $\mu\approx
2/9$. However, the uncertainty of the exponent is rather large, as the value
changes quite substantially depending on what data are included during the
fitting. (We excluded the point at the lowest pressure.) Here, we would only
want to state a rough upper bound of 1/2 for the exponent. Moreover, we would
argue that the decrease in viscosity arises from a much reduced steric
repulsion as molecules slide past each other under a tensile hydrostatic
stress.

We do not intend to claim that the functional form of $\eta(p)$ proposed in
Eq. (9) is superior to other descriptions, in particular generalizations of
the McEwen equation McEwen1952JIP , for example by Bair Bair2015HTHP  or other
others discussed by recently by Kondratyuk, Pisarev, and Ewen
Kondratyuk2020JCP . However, our analysis suggests that a quasi-exponential
pressure dependence at large compression might be an appropriate description,
at least at low shear rates.

In fact, we find the McEwen $\eta(p)$ dependence, which can be cast as the
$\eta_{2}=0$ limiting case of Eq. (9) to hold extremely well at compressive
stress and fixed, large shear rates, as is demonstrated in the main part of
Fig. 10. Here $\eta_{\textrm{ref}}$ is the effective viscosity at a high shear
rate and zero normal stress, the exponent $\mu$ is set to $\mu=1/2$, and
$\sigma_{0}$ a hypothetical tensile stress at which the McEwen formula would
predict the viscosity to disappear. (Using $\mu=1/2$ substantially worsened
the agreement between fit and data.) However, without further modification,
the McEwen equation makes poor predictions once the stress becomes tensile, as
can be seen in the inset of Fig. 10.

Figure 10: Effective viscosity at a high shear rate of $\dot{\nu}=10^{11}$ Hz
at $T=500$ K as a function of pressure $p_{33}$. The inset highlights result
for tensile stresses.

###  4.3 Chain-length dependence

In addition to changing pressure and temperature, we also altered the degree
of polymerization, $P$, however, only in a relatively narrow range in which
molecules are much too short to entangle. From the results reported in Table
1, it can be seen that $\eta_{0}$ and $\dot{\gamma}_{0}$ scale roughly with
$P^{2}$ and $1/P^{2}$, respectively. Changes in the exponents $n$ and $a$ are
rather small, which points once more to the existence of an elementary
instability. This result does not necessarily contradict findings by Sivabeak
and Persson (SP) Sivebaek2012PRL , who suggested a more noticeable $n(P)$
dependence, where their symbol $n$ corresponded to our $1-n$. Further, more
substantial difference between their and our study are that SP explored a
relatively small normal pressure of 10 MPa while altering $P$ by a factor of
70, which covers the range from oligomer to strongly entangled polymers.

An interesting feature born out in our data is that $\eta(\dot{\gamma})$
overlaps at very high shear rates for different degree of polymerization $P$,
i.e., not only the exponent $n$ but also the prefactors are similar at high
shear rates for different alkanes.

Figure 11: Scaled shear stress vs. shear rate for $n$-hexadecane at 500 K and
300 MPa.

###  4.4 Microscopic analysis

Since the early days of shear-thinning studies, attempts were made to
rationalize phenomenological relations in terms of microscopic models and to
relate shear thinning to structural changes Lacks2001PRL . Early theories of
polymers in solution Lodge1956TFS , “derive” the Carreau equation by
considering the rate of formation and breaking of contact points. Even if the
explanation were valid for such systems, it can scarcely matter for dense
polymer mixtures, when the resistance to flow is largely affected by steric
repulsion. Yet, non-local configurational changes like polymer alignment have
been held responsible for changes in viscosity Jadhao2019TL . Specifically,
polymers oriented along the shear-flow direction were expected to counteract
shear in that direction than those with quasi-spherical equilibrium shapes.

It might have came as a surprise to many when Jadhao and Robbins Jadhao2019TL
reported the alignment of polymers to have saturated when the viscosity had
merely changed by a factor of two, which, in their master plot was at the very
beginning of the shear-thinning regime. We can confirm their finding in our
work and extend on it by noticing that the molecules tend to coil up again at
constant pressure, rather than constant volume, at very high shear rates.
Thus, Ostwald scaling may occur even if the overall molecular shape develops
in a non-monotonic fashion.

Figure 12: Different cartesian components to $R_{\textrm{g}}^{2}$ vs shear
rate, for $n$-hexadecane at $p=300$ MPa and $T=500$ K (full symbols). Open
symbols refer to gas-phase data at $p=3.4$ MPa.

For reason of completeness, we also present data that we produced to determine
the diffusion coefficient from equilibrium simulations in Fig. 13, which shows
the atom averaged mean-square displacement as a function of time.
Unfortunately, the value for the diffusion coefficient is sensitive to what
functional form is assumed. Assuming
$\left\langle\left\\{x(t)-x(0)\right\\}^{2}\right\rangle$ to be a constant
plus $Dt$ for the large times yields slightly different results than when
using the sum of a sub-diffusive and a diffusive function. Nonetheless, we
note that the time, where sub-diffusive motion crosses over to diffusive
motion is of similar order of magnitude as the inverse of the shear rate, at
which shear-thinning becomes noticeable.

Figure 13: Mean-squared displacement (MSD) of $n$-hexadecane from equilibrium
MD simulation at 300 MPa. Lines are drawn to guide the eye.

##  5 Discussion and Conclusions

This work demonstrates that the shear-thinning of alkanes at high pressure can
be mapped quite well onto the Prandtl model. Fixing its two dimensional
parameters at a given temperature and pressure even allows changes in the
equilibrium viscosity and the functional form of $\eta(\dot{\gamma})$ to be
reproduced over a reasonably broad temperature range without readjusting
parameters. To a large degree, the parametrization is even transferable
between different alkanes. Moreover, the obtained activation barrier of
$\Delta E\approx 0.08$ eV that we obtain at a pressure of $p=300$ MPa is
greater but of similar order of magnitude as experimental result for the
relaxational dynamics in polyethylen melts at zero pressure, which are close
to 0.05 eV at ambient pressure Qiu1999M . Despite this success, we do not
claim our parametrization to be optimal but rather minimal in the sense that
the smallest number of adjustable parameters needed to describe shear-thinning
was used. In principle, it could be attempted to also target the linear-
response frequency dependence $\eta(\omega)$.

Our simulation data on all-atom alkanes, in particular $n$-hexadecane also
indicates that the additional parameter(s) used in a Carreau-Yasuda fit are
justified, at least at elevated pressure. We find a fit to have a median error
of roughly 12%, which reduces to about 6% when using an additional parameter
in the Cross model. This error is cut once more into one half when using
Carreau instaed of Cross, and by almost another factor of two with Carreau-
Yasuda, at which point the systematic error is on par with our stochastic
error.

## Acknowledgments

We thank Jörg Baschnagel and Daniele Dini for useful discussions.

## Declarations

### Funding

This research was supported by the German Research Foundation (DFG) under
grant number GA 3059/2-1.

### Conflict of interest

No.

### Authors’ contributions

MHM designed and supervised the project, wrote the manuscript. HG carried out
all the simulations, analyzed the data and prepared the figures. MHM and HG
reviewed and modified the manuscript.

## References

  * [1] S. Bair.  Purported fragile-to-arrhenius crossover in squalane.  Proc. Natl. Acad. Sci. U.S.A., 114(42), October 2017.
  * [2] S. Bair and P. Kottke.  Pressure-viscosity relationships for elastohydrodynamics.  Tribol. Trans., 46(3):289–295, January 2003.
  * [3] S. Bair, P. Vergne, P. Kumar, G. Poll, I. Krupka, M. Hartl, W. Habchi, and R. Larsson.  Comment on “history, origins and prediction of elastohydrodynamic friction” by spikes and jie.  Tribol. Lett., 58(1), March 2015.
  * [4] Scott Bair.  Choosing pressure-viscosity relations.  High Temperatures-High Pressures, 44:415–428, 2015.
  * [5] Scott Bair, Clare McCabe, and Peter T. Cummings.  Comparison of nonequilibrium molecular dynamics with experimental measurements in the nonlinear shear-thinning regime.  Phys. Rev. Lett., 88:058302, Jan 2002.
  * [6] Howard A Barnes.  Thixotropy—a review.  Journal of Non-Newtonian Fluid Mechanics, 70(1-2):1–33, May 1997\.
  * [7] Pierre J. Carreau.  Rheological equations from molecular network theories.  Transactions of the Society of Rheology, 16(1):99–127, March 1972\.
  * [8] A. W. Crook and Thomas Edward Allibone.  The lubrication of rollers iv. measurements of friction and effective viscosity.  Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 255(1056):281–312, 1963.
  * [9] Peter J. Daivis and Denis J. Evans.  Comparison of constant pressure and constant volume nonequilibrium simulations of sheared model decane.  The Journal of Chemical Physics, 100(1):541–547, 01 1994.
  * [10] Ingrid de Almeida Ribeiro and Maurice de Koning.  Non-newtonian flow effects in supercooled water.  Phys. Rev. Res., 2(2), April 2020.
  * [11] James P. Ewen, Hongyu Gao, Martin H. Müser, and Daniele Dini.  Shear heating, flow, and friction of confined molecular fluids at high pressure.  Phys. Chem. Chem. Phys., 21:5813–5823, 2019.
  * [12] James P. Ewen, Hugh A. Spikes, and Daniele Dini.  Contributions of molecular dynamics simulations to elastohydrodynamic lubrication.  Tribology Letters, 69:24, 2021.
  * [13] Henry Eyring.  Viscosity, Plasticity, and Diffusion as Examples of Absolute Reaction Rates.  The Journal of Chemical Physics, 4(4):283–291, 12 1936.
  * [14] Daniel S. Fisher.  Sliding charge-density waves as a dynamic critical phenomenon.  Physical Review B, 31(3):1396–1427, February 1985.
  * [15] V. Jadhao and M. O. Robbins.  Reply to bair: Crossover to arrhenius behavior at high viscosities in squalane.  Proc. Natl. Acad. Sci. U.S.A., 114(42), 2017.
  * [16] Vikram Jadhao and Mark O. Robbins.  Probing large viscosities in glass-formers with nonequilibrium simulations.  Proc. Natl. Acad. Sci. U.S.A., 114(30):7952–7957, 2017.
  * [17] Vikram Jadhao and Mark O. Robbins.  Rheological properties of liquids under conditions of elastohydrodynamic lubrication.  Tribology Letters, 67:66, 2019.
  * [18] Nikolay D. Kondratyuk, Vasily V. Pisarev, and James P. Ewen.  Probing the high-pressure viscosity of hydrocarbon mixtures using molecular dynamics simulations.  The Journal of Chemical Physics, 153(15):154502, 10 2020.
  * [19] Daniel J. Lacks.  Energy landscapes and the non-newtonian viscosity of liquids and glasses.  Phys. Rev. Lett., 87:225502, Nov 2001.
  * [20] A W Lees and S F Edwards.  The computer study of transport processes under extreme conditions.  Journal of Physics C: Solid State Physics, 5(15):1921, aug 1972\.
  * [21] Claire A. Lemarchand, Nicholas P. Bailey, Billy D. Todd, Peter J. Daivis, and Jesper S. Hansen.  Non-Newtonian behavior and molecular structure of Cooee bitumen under shear flow: A non-equilibrium molecular dynamics study.  The Journal of Chemical Physics, 142(24):244501, 06 2015.
  * [22] R Li.  Time-temperature superposition method for glass transition temperature of plastic materials.  Mater. Sci. Eng., A278:36–45, 2000.
  * [23] A. S. Lodge.  A network theory of flow birefringence and stress in concentrated polymer solutions.  Trans. Faraday Soc., 52:120–130, 1956.
  * [24] Pierre Lucas.  Fragile-to-strong transitions in glass forming liquids.  J. Non-Cryst. Solids, 557:119367, April 2021.
  * [25] E. McEwen.  The effect of variation of viscosity with pressure on the load-carrying capacity of the oil film between gear-teeth.  J. Inst. Pet., 38(1):646–672, 1952.
  * [26] Martin H. Müser.  Velocity dependence of kinetic friction in the prandtl-tomlinson model.  Physical Review B, 84(12), September 2011.
  * [27] Martin H. Müser.  Shear thinning in the Prandtl model and its relation to generalized Newtonian fluids.  Lubricants, 8(4):38, March 2020.
  * [28] Wolfgang Ostwald.  Ueber die geschwindigkeitsfunktion der viskosität disperser systeme. i.  Kolloid-Zeitschrift, 36(2):99–117, February 1925.
  * [29] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of Computational Physics, 117(1):1–19, 1995.
  * [30] V.L. Popov and J.A.T. Gray.  Prandtl-tomlinson model: History and applications in friction, plasticity, and nanotechnologies.  ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 92(9):683–708, July 2012\.
  * [31] L. Prandtl.  Ein gedankenmodell zur kinetischen theorie der festen körper.  J. Appl. Math. Mech., 8(2):85–106, 1928.
  * [32] Melissa L. P. Price, Dennis Ostrovsky, and William L. Jorgensen.  Gas-phase and liquid-state properties of esters, nitriles, and nitro compounds with the opls-aa force field.  Journal of Computational Chemistry, 22(13):1340–1352, 2001.
  * [33] XiaoHua Qiu and M. D. Ediger.  Local and global dynamics of unentangled polyethylene melts by 13c nmr.  Macromolecules, 33(2):490–498, December 1999.
  * [34] C. M. Roland, S. Bair, and R. Casalini.  Thermodynamic scaling of the viscosity of van der waals, h-bonded, and ionic liquids.  J. Chem. Phys., 125(12):124508, September 2006.
  * [35] Andy Ruina.  Slip instability and state variable friction laws.  Journal of Geophysical Research: Solid Earth, 88(B12):10359–10370, December 1983.
  * [36] E. Schalek and A. Szegvari.  Die langsame koagulation konzentrierter eisenoxydsole zu reversiblen gallerten.  Kolloid-Zeitschrift, 33(6):326–334, December 1923.
  * [37] Shirley W. I. Siu, Kristyna Pluhackova, and Rainer A. Böckmann.  Optimization of the opls-aa force field for long hydrocarbons.  Journal of Chemical Theory and Computation, 8(4):1459–1470, 2012\.  PMID: 26596756.
  * [38] I. M. Sivebaek, V. N. Samoilov, and B. N. J. Persson.  Effective viscosity of confined hydrocarbons.  Physical Review Letters, 108(3):036102, January 2012.
  * [39] Elizabeth D. Smith, Mark O. Robbins, and Marek Cieplak.  Friction on adsorbed monolayers.  Physical Review B, 54(11):8252–8260, September 1996.
  * [40] H. Spikes and Z. Jie.  History, origins and prediction of elastohydrodynamic friction.  Tribol. Lett., 56(1):1–25, August 2014.
  * [41] H. Spikes and W. Tysoe.  On the commonality between theoretical models for fluid and solid friction, wear and tribochemistry.  Tribol. Lett., 59(1), June 2015.
  * [42] H. Spikes and J. Zhang.  Reply to the comment by scott bair, philippe vergne, punit kumar, gerhard poll, ivan krupka, martin hartl, wassim habchi, roland larson on “history, origins and prediction of elastohydrodynamic friction” by spikes and jie in tribology letters.  Tribol. Lett., 58(1), March 2015.
  * [43] H. A. Spikes.  Comment on: Rheology of an ionic liquid with variable carreau exponent: A full picture by molecular simulation with experimental contribution, by nicolas voeltzel, philippe vergne, nicolas fillot, nathalie bouscharain, laurent joly, tribology letters (2016) 64:25.  Tribol. Lett., 65(2), April 2017.
  * [44] Billy D. Todd and Peter J. Daivis.  Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications.  Cambridge University Press, 2017.
  * [45] N. Voeltzel, P. Vergne, N. Fillot, N. Bouscharain, and L. Joly.  Rheology of an ionic liquid with variable carreau exponent: A full picture by molecular simulation with experimental contribution.  Tribol. Lett., 64(2), October 2016.
  * [46] Nicolas Voeltzel, Philippe Vergne, Nicolas Fillot, Nathalie Bouscharain, and Laurent Joly.  Reply to the “comment on ‘rheology of an ionic liquid with variable carreau exponent: A full picture by molecular simulation with experimental contribution, ’ by n. voeltzel, p. vergne, n. fillot, n. bouscharain, l. joly, tribology letters (2016) 64:25” by h. a. spikes.  Tribol. Lett., 65(2), April 2017.