# ADL-ID: Adversarial Disentanglement Learning Network for RF Fingerprinting
Temporal Domain Adaptation

Abdurrahman Elmaghbub, Bechir Hamdaoui and Weng-Keen Wong  School of
Electrical Engineering and Computer Science, Oregon State University
{elmaghba, hamdaoui<EMAIL_ADDRESS>
# ADL-ID: Adversarial Disentanglement Learning for Wireless Device
Fingerprinting Temporal Domain Adaptation

Abdurrahman Elmaghbub, Bechir Hamdaoui and Weng-Keen Wong  School of
Electrical Engineering and Computer Science, Oregon State University
{elmaghba, hamdaoui<EMAIL_ADDRESS>
###### Abstract

As the journey of 5G standardization is coming to an end, academia and
industry have already begun to consider the sixth-generation (6G) wireless
networks, with an aim to meet the service demands for the next decade. Deep
learning-based RF fingerprinting (DL-RFFP) has recently been recognized as a
potential solution for enabling key wireless network applications and
services, such as spectrum policy enforcement and network access control. The
state-of-the-art DL-RFFP frameworks suffer from a significant performance drop
when tested with data drawn from a domain that is different from that used for
training data. In this paper, we propose ADL-ID, an unsupervised domain
adaption framework that is based on adversarial disentanglement representation
to address the temporal domain adaptation for the RFFP task. Our framework has
been evaluated on real LoRa and WiFi datasets and showed about $24\%$
improvement in accuracy when compared to the baseline CNN network on short-
term temporal adaptation. It also improves the classification accuracy by up
to $9\%$ on long-term temporal adaptation. Furthermore, we release a $5$-day,
$2.1$TB, large-scale WiFi 802.11b dataset collected from $50$ Pycom devices to
support the research community efforts in developing and validating robust
RFFP methods.

###### Index Terms:

RF Fingerprinting, Domain Adaptation, Disentanglement Learning, WiFi Dataset,
IoT Testbed.

##  I Introduction

An important highlight of evolving next-generation wireless networks is
security and privacy, which have been somehow overlooked, to some extent, in
previous generations [1]. While traditional security relies heavily on the
upper-layer mechanisms, it is essential to complement these mechanisms with
unclonable, physical-layer security technologies to increase the security
robustness of such systems. As an efficient physical-layer security mechanism,
RF fingerprinting (RFFP) technology allows for device identification based on
unique fingerprints inherited in the RF signal. Due to the imperfection of the
manufacturing process, circuit components in the signal path introduce a
variety of random hardware impairments such as I/Q mismatch, phase noise, and
non-linearity [2] that collectively impact the transmitted signal,
constituting a unique RF Fingerprint (RFF). Deep Learning (DL) has reduced the
need for preprocessing and expert-defined feature extraction techniques,
building state-of-the-art RFFP frameworks [3, 4]. However, most of these rely
on the assumption that the training and testing sets are drawn from the same
distribution, which falls short of the conditions of realistic RF scenarios.

We define domain portability as the ability of a model to maintain its
training domain performance when applied to new domains. In the RFFP context,
the domain refers to the settings in which data has been collected. This
includes time, channel, receiver, protocol, and environment. Hence, any
considerable change in the training settings yields a different domain. In [5,
6], we conducted an experimental study on LoRa devices, which disclosed the
sensitivity of DL-RFFP to domain changes. In that work, the models (i) fail to
maintain their high accuracy when channel conditions change and (ii)
completely lose their classification ability when the protocol configuration
or receiver changes. This lack of portability represents a substantial
impediment to the implementation of DL-based RFFP methods in real-world
applications.

The complexity of the RFFP problem comes from the fact that we are dealing
with two layers of obscurity: the fingerprint interactions in both circuit and
channel levels and the feature vector learned by the deep learning network.
The lack of behavioral comprehension of RFFs and deep learning networks opens
up the floor to several hypotheses to explain the exposed sensitivity of deep
learning models when the domain changes. It could be that the fingerprints
themselves change as the domain changes; hence, identifying devices based on
their old/learned fingerprints yields poor accuracy in cross-domain
validation. This assumption disqualifies these features as fingerprints since
authentication schemes require them to be robust and permanent. More
reasonably, the feature vector learned by the network includes both RFFs and a
domain-specific component which leads to confusion when the latter component
changes as the domain changes. Most of the proposed solutions adapted the
latter assumption in their endeavors. Some works have focused on eliminating
the background effects from the raw signal to improve its domain
generalization. One such example is the use of channel equalization to improve
the performance by eliminating channel influences [7], but this is a partial
solution since they eliminate many beneficial spurious emissions in the
passband, which reduces the user capacity of the RFFs [8]. On the other hand,
DeepLoRa [9], a data augmentation technique based on the ITU-R channel models,
was proposed to counter the degradation introduced by the wireless channel.
The best performance reported by DeepLoRa was $36\%$ with only $5$ devices
being classified, achieved when the model is trained on one day and tested on
another day. This accuracy drops to $8\%$ when the number of devices goes to
$50$.

More recently, the portability problem of the RFFP has been formulated as a
domain adaptation problem [10, 11], allowing us to leverage the advancements
of this branch of transfer learning in the RF context. This formulation
assumes that the source and target domains have slightly different
distributions. In [12], a multi-discrepancy domain adaptation method has been
proposed to minimize the intraclass discrepancy across the source and target
domains while maximizing the inter-class discrepancy, leading to $15\%$
improvement over the ResNet network. Adversarial domain adaptation [13] has
also been proposed to reduce the influence of domain features on the RFFs by
maximizing the domain discrimination loss, resulting in an enhancement in the
portability of the network. However, the target domains on the two mentioned
domain adaptation methods have been created via simulation by introducing
white Gaussian noise to the source domain dataset to create several target
domains with different signal-to-noise ratio conditions. These simulation
settings are far from realistic RF scenarios where multiple settings change at
the same time. Addressing multiple varying settings, Tweak [11] combined the
metric learning with a lightweight calibration method to address the
portability problem in hardware, channel, and configuration dimensions,
separately. Tweak shows a considerable performance improvement in each of
these dimensions, but needs to be recalibrated each time the target domain
changes.

Another exciting domain adaptation category is disentangled representation
learning, an unsupervised learning technique that separates (disentangles)
features into disjoint parts of the underlying representation[14]. Motivated
by the intuition that the RFFs are entangled with domain-specific features and
inspired by the disentanglement concept [15], we propose ADL-ID, a domain
adaptation framework that integrates disentanglement representation learning
with adversarial learning to address the temporal portability problem in RFFP.
Our framework factors the feature vector into a device-specific component
(RFFs) and a domain-specific component. Then, it uses only the device-specific
component, which contains domain-invariant features, as input to the
classifier. Our proposed framework, tested using real WiFi dataset (detailed
in Section II) and LoRa datasets (described in [16]), provides a high-
performance gain, $20\%-24\%$, in short-term-varied domain adaptation. We are
also excited to release our large-scale WiFi 802.11b dataset of $50$ Pycom
devices to the research community to support the efforts toward enabling RFFP
portability. This is the first public WiFi 802.11b RFFP dataset, which is
common among low-rate IoT devices. Our testbed includes the same devices used
before to build and release the LoRa dataset [16], which extends the
opportunity to analyze other domain dimensions, such as protocol
configurations and hardware receivers. Our main contributions are as follows:

  * •

Proposed an adversarial-based disentanglement learning network to address the
RFFP temporal portability issue.

  * •

Released the first large-scale WiFi 802.11b RFFP dataset of $50$ Pycom devices
over $5$ consecutive days for both indoor and outdoor scenarios.

  * •

Provided novel insights into the performance of RFFP methods due to short-term
and long-term temporal domain variations using both WiFi and LoRa datasets.

The rest of the paper is organized as follows. Section II describes the device
testbed and WiFi dataset. Following that, Section III reveals the
investigation on short-term variations performance. Then, our framework is
explained in Section IV while our results are presented and discussed in
Section V. The paper is concluded in Section VI.

##  II Testbed and WiFi Dataset Description

Figure 1: 50-device Pycom WiFi testbed

In this section, we begin by describing the hardware, software, and protocol
components of the tested used for collecting our dataset. Then, we describe
the collected WiFi dataset. The description and download information of the
LoRa dataset that we also used in our evaluation can be found in [16, 5].

###  II-A Testbed Description

####  II-A1 Hardware

Our testbed, shown in Fig. 1, consists of $50$ Pycom devices ($25$ FiPy & $25$
LoPy boards on top of $50$ PySense extension boards). These boards mainly
include ESP$32$-D$0$WDQ, Semtech SX$1276$, and Sequans Monarch chips to
support WiFi, Bluetooth, LoRa, Sigfox, and LTE networks, respectively. On the
reception side, we used an Ettus USRP (Universal Software Radio Peripheral)
B$210$ with a VERT$900$ antenna for the data acquisition. The RF-front of the
USRP was configured with a center frequency of $2.412$GHz, a $20$MHz
Bandwidth, and a $25$MSps sample rate. We used lipo batteries to power the
Pycom devices. Each Pycom device was connected to a dedicated WiFi antenna and
configured to transmit IEEE802.11b packets in the $2$GHz spectrum.

####  II-A2 Software

Pycom devices are based on FREERTOS operating system. To perform transmission,
we programmed and configured our Pycom boards using MicroPython, an efficient
implementation of Python3 composed of a subset of standard Python libraries
and optimized to run on microcontrollers and constrained environments. Also,
we used Pymakr plugin as a REPL console that connects to Pycom boards to run
codes or upload files. To perform reception, we used the GNURadio software, a
real-time signal processing graphical tool, along with USRP Hardware Driver
(UHD) to set up and configure the USRP receiver to capture WiFi transmissions,
plot their time and frequency domains, and store the collected samples in
their files.

####  II-A3 Protocol

We used the WiFi 802.11b protocol with the high-rate direct-sequence spread-
spectrum (HR/DSSS) physical layer. Transmitters were configured to transmit at
$1$Mbps, with a carrier frequency of $2.412$GHz and $20$MHz bandwidth.

###  II-B Dataset Description

Our WiFi dataset contains $2.1$TB of WiFi transmissions of $50$ Pycom devices
captured over $5$ consecutive days in both indoor (laboratory) and outdoor
environments. For each day, we capture $5$ transmissions from each device in a
round-robin fashion, where each capture consists of $50$M complex-valued I/Q
samples. The time-gap between two consecutive captures of the same device is
$5$mins. Transmitters were located about $5$ meters away from the access point
and the receiver, with a clear line of sight. The collected signals are down-
converted by the USRP to the baseband and then stored in the form of I/Q data
in both time and frequency domain representations. For each capture,
associated metadata compliant with the latest SigMF specifications and
additional field extensions for usability are recorded.

More details and use cases of the WiFi dataset can be found in [5]. The
datasets can be downloaded from NetSTAR Lab at
https://research.engr.oregonstate.edu/hamdaoui/datasets/http://research.engr.oregonstate.edu/hamdaoui/datasets.

##  III The Temporal Portability Problem

Several works [5, 17, 9] consider the days as the change unit in the temporal
dimension of the domain portability problem. Hence, they analyze the
performance of the models when the training and testing data are from
different days. However, there is nothing special about the 24 hours as far as
the RF fingerprints are concerned. We postulate that the data distribution of
fingerprints changes dramatically even within the same day. To validate this
hypothesis, we used our LoRa and WiFi datasets in which we have several
captures of the same device within the same day. Starting with the LoRa
dataset [16], the time gap between two consecutive captures is $5$ minutes
with $5$m distance between transmitters and the receiver. The captures were
collected in indoor and outdoor scenarios, and the devices have not been
rebooted between the captures. We trained a 6-Layer CNN network [5] using
capture $1$ and tested it on data from captures $1$, $2$, $3$, and $4$ of the
$4$ days. Surprisingly, the model’s performance suffers from a huge
degradation even when the time gap between the training and testing sets is as
short as $5$ minutes. Furthermore, as we can see in Fig 2—depicting accuracy
results when using the LoRa data, the performance continues to degrade as we
move further in time from the training domain. This trend manifests itself in
the four days and both scenarios.

(a) LoRa Indoor Scenario. (b) LoRa Outdoor Scenario.

Figure 2: Performance in short-term variations of LoRa captures where the time
gap between two consecutive captures is $5$mins. The CNN model was trained on
captures $1$ of each day and tested on captures $1$, $2$, $3$, and $4$.

(a) WiFi Indoor Scenario. (b) WiFi Outdoor Scenario.

Figure 3: Performance in short-term variations of WiFi captures where the time
gap between two consecutive captures is $5$mins. The CNN model was trained on
captures $1$ of each day and tested on captures $1$, $2$, $3$, and $4$. Figure
4: Performance in long-term variations on LoRa (1st column) and WiFi (2nd
column) captures where the captures are collected on $4$ different days. The
CNN model was trained on capture $1$ and tested on captures $1$, $2$, $3$, and
$4$.

To confirm our observation, we carried out the same investigation in our
50-device WiFi dataset, which also includes $5$ captures of the same devices
within the same day for $5$ consecutive days, where the time gap is similar to
that in the LoRa dataset. The difference is that the devices here have been
rebooted between the captures. Results in Fig. 3 show that the WiFi dataset is
no different, and the model also suffers from short-term variations on the
WiFi dataset. These results consolidate our hypothesis that even a slight
departure from the network’s training domain can introduce confusion to the
classifier and hurt its performance severely.

We also notice that on most of the days, the most significant drop rate in
performance occurs in the transition from capture $1$ to capture $2$, and that
rate decreases as we move away from the training domain. We claim that the
drop rate decreases until it reaches a point when the time gap does not matter
anymore; i.e., aging has no to little effect. Hence, with respect to domain
adaptation performance, $2$-day or $10$-day time separation between the source
and target domains has a similar effect. Some of that can be seen in Fig. 4,
which indicates that the effect on the different days is similar. We refer to
this effect as the effect of long-term variations. As expected, Fig. 4 shows
that the drop in the network’s performance due to the long-term variations is
more severe than the previous scenario as the accuracy drops from around
$90\%$ to $20\%$, when the model is trained on capture $1$ and tested on
captures $2$, $3$, and $4$, which have been collected on day $2$, $3$, and
$4$, respectively.

##  IV The Proposed Framework: ADL-ID

Figure 5: Architecture of the proposed unsupervised adversarial-based
disentanglement learning network for the RFFP task.

The proposed framework, shown in Fig. 5, aims to train an encoder, called
fingerprint encoder, $E_{F}$, to extract only the device-specific feature
vector (RFFs) and feed it to the cross-domain device classification task. The
intuition here is that the RFFs are shared across domains while the background
and domain artifacts are different from one domain to another. This is
achieved by having three encoders: fingerprinting encoder, $E_{F}$, source
encoder, $E_{S}$, and target encoder, $E_{T}$ that are trained to model
different components of the domain representation. Using labeled source domain
data, $X_{S}=\\{(x_{s}^{i},y_{s}^{i})\\}^{N_{s}}_{i=0}$, and unlabeled target
domain data, $X_{T}=\\{(x_{t}^{i})\\}^{N_{t}}_{i=0}$, the fingerprint encoder
learns to extract discriminative and shared features across the source and
target domains while the source and target encoders learn the features that
are specific to the source and target domains, respectively.
At each training iteration, the output of the fingerprinting encoder will be
fed to the domain discriminator, $D(x)$, to eliminate domain-indicative
features to make the RFFs similar to each others regardless of their original
domains. We minimize the domain classification loss, $L_{Domain}$, while
having the gradient reversal layer (GRL) [18] in the data path. GRL acts as an
identity transformation in the forward path while reversing the gradient’s
sign in the backpropagation. Hence, it maximizes the confusion of the domain
discriminator as the updated weights of $E_{F}$ result in RFFs with less
domain-indicative information. To enhance the separation between the RFF and
the domain-specific representations of an input frame as the output of $E_{F}$
and $E_{S}$ or $E_{T}$ encoders, we apply a difference loss, $L_{difference}$
that encourages the orthogonality of the two representations and pushes them
further away from each other by minimizing the squared Frobenius norm of the
multiplication of two matrices $H_{f}$ and $H_{s}$ or $H_{t}$ whose rows are
$h_{f}$ and $h_{s}$ or $h_{t}$, respectively [15]. To avoid trivial solutions
from the previous losses, we use a signal decoder, $G(x)$, that decodes the
combination of the RFF and domain-specific features to get back the original
input frame. The reconstruction process is governed by the reconstruction
loss, $L_{Reconstruct}$, which is a scale–invariant mean squared error (SI-
MSE) loss that experimentally proves to be more effective for disentanglement
representations [19]. Finally, $h_{f}^{s}$, the RFF representation of $X_{s}$
is passed to the main-task classifier with a cross-entropy loss, $L_{Class}$,
that encourages the RFF representation to not just be domain-invariant and
separable from the domain-specific representation, but also to be
discriminative with regard to the device classification task. The total loss
of our network can thus be written as

$L_{total}=L_{Class}+\alpha L_{Reconstruct}+\beta L_{Difference}+\gamma
L_{Domain}$

where the hyperparameters, $\alpha$, $\beta$, and $\gamma$, control the
strength of these losses in the overall training process.

Each encoder in our framework consists of $6$ convolution blocks followed by a
fully-connected layer. Each convolution block contains 2D convolution, batch
normalization, leakyReLu, and max pooling layers. The signal decoder has also
the same components in an opposite order. On the other hand, the device
classifier and domain discriminator has three and two fully-connected layers,
respectively, followed by a Softmax layer. We start the training with a warm
up phase of $10000$ steps where only the supervised task and the
reconstruction losses are used before switching to the other losses. This can
help in speeding up the training process as well as leading to a better
generalization [20].

##  V Performance Evaluation and Discussion

(a) LoRa Dataset. (b) WiFi Dataset.

Figure 6: ADL-ID vs. CNN on captures from Day $1$, Indoor scenario for both
(a) LoRa and (b) WiFi. ADL-ID was trained on labeled capture $1$ data (source)
and unlabeled capture $2$ data (target) and tested on captures $1$, $2$, $3$,
and $4$.

(a) LoRa Dataset. (b) WiFi Dataset.

Figure 7: ADL-ID vs. CNN on captures from 4 different days, Indoor scenario
for both (a) LoRa and (b) WiFi. ADL-ID was trained on labeled capture $1$ data
(source) and unlabeled capture $2$ data (target) and tested on captures $1$,
$2$, $3$, and $4$.

In this section, we evaluate our proposed framework for the indoor setup of
our LoRa and WiFi datasets. For each device, we used $20,000$ I/Q frames each
of size $2$x$1024$, partitioned into training, validation, and testing sets.
Our network has been trained for $50$ epochs on an Nvidia DGX2 node. We use
two testing methodologies to run our evaluation: (i) testing on same-day
captures to evaluate the short-term temporal adaptation; and (ii) testing on
different-days captures to evaluate the long-term temporal adaptation.

First, we start with the short-term temporal adaption scenario. In the LoRa
dataset, the time gap between the source domain and target domain is $5$mins,
and the transmitters have not been rebooted between the captures.
Specifically, we train our framework with labeled data from capture 1 (source
domain) and unlabeled data from capture 2 (target domain) and then test it on
captures $1$, $2$, $3$, and $4$, all captured on the same day with $5$mins in
between. Compared to the CNN baseline, Fig. 6a shows that ADL-ID (in red),
when tested on LoRa data, provides a $20\%-24\%$ improvement in the
classification accuracy, pushing the performance from $57\%$ to $76\%$ on
capture $2$, from $41\%$ to $65\%$ on capture $3$, and from $29\%$ to $53\%$
on capture $4$. This improvement can be linked to the efforts of our network
to minimize the presence of domain-specific features from the source and
target domains in the RFF feature vector, which enhances the network’s
generalization. On the WiFi dataset (Fig. 6b), our proposed framework
continues to show an improvement in accuracy with highest gain of $9\%$,
increasing the accuracy from $55\%$ to $64\%$ on capture $2$, and from $40\%$
to $42\%$ on capture $3$, and from $22\%$ to $29\%$ on capture $4$. It can be
noticed that the improvement percentage is higher in the LoRa dataset. We
speculate that one reason could be that, in contrast to the LoRa dataset, the
transmitters in the WiFi scenario were rebooted between the different
captures. Consequently, a partial RF calibration has taken place, which
affects the RFF distribution and stretches the distance between the source and
target domains.

Second, we evaluate our network on long-term temporal adaptation.
Specifically, we consider day $1$ data as a source domain and day $2$ data as
the target domain. Hence, we train the network with labeled day $1$ data and
unlabeled day $2$ data and then test it on data from the $4$ days. While Fig.
7 shows that our network provides an improvement over the baseline CNN network
in the long-term temporal adaptation as well, it must be marked that it did
not succeed in bringing the performance on long-term varying domains closer
enough to the source domain’s performance. A reason for this is that our
framework struggles to adapt to new domains when the difference between the
source and target domains is too severe, which is the case in the long-term
variation scenario. This results also conveys that performance of the
framework is dependent on the similarity between the source and target
domains. It is worth mentioning that similar behavior has been reported when
domain adaptation techniques are applied in vision applications [20]. The
variations in target domains in real-world RF applications are exceptionally
more complicated compared to the vision domain. In the RFFP context, the
variations in the target domain might include changes in the temporal
dimension, channel characteristics, protocol configurations, and receiver
hardware variations. Furthermore, these variations have made the interaction
between the RFFs and domain-specific features more subtle, which results in a
more difficult disentanglement process.

##  VI Conclusion

We presented a large-scale WiFi 802.11b dataset, collected from 50 Pycom
devices, and analyzed the performance of CNN-based RFFP on short- and long-
term variations using LoRa and WiFi data. We also proposed ADL-ID, which
provides up to $24\%$ improvement in short-term temporal adaptations and up to
$10\%$ improvement in long-term variations over the baseline CNN. These
results prove the potential of disentanglement representation techniques in
enabling robust RFFP. However, more work is still needed to fully understand
the RFF behaviors when the domain changes and to develop more robust RFFP
methods.

##  VII Acknowledgment

This research is supported in part by Intel/NSF MLWiNS Award No. 2003273. We
would like to thank Intel researchers, Dr. Kathiravetpillai Sivanesan, Dr.
Lily Yang, Dr. Richard Dorrance, Dr. Vesh Raj Sharma Banjade, and Dr. Hechen
Wang for their constructive discussions and feedback.

## References

  * [1] S. Dang, O. Amin, B. Shihada, and M.-S. Alouini, “What should 6g be?” _Nature Electronics_ , vol. 3, no. 1, pp. 20–29, 2020.
  * [2] A. Elmaghbub, B. Hamdaoui, and A. Natarajan, “Widescan: Exploiting out-of-band distortion for device classification using deep learning,” in _GLOBECOM 2020-2020 IEEE Global Communications Conference_. IEEE, 2020, pp. 1–6.
  * [3] K. Merchant, S. Revay, G. Stantchev, and B. Nousain, “Deep learning for rf device fingerprinting in cognitive communication networks,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 12, no. 1, pp. 160–167, 2018.
  * [4] B. Hamdaoui, A. Elmaghbub, and S. Mejri, “Deep neural network feature designs for RF data-driven wireless device classification,” _IEEE Network_ , vol. 35, no. 3, pp. 191–197, 2020.
  * [5] A. Elmaghbub and B. Hamdaoui, “LoRa device fingerprinting in the wild: Disclosing RF data-driven fingerprint sensitivity to deployment variability,” _IEEE Access_ , 2021.
  * [6] B. Hamdaoui and A. Elmaghbub, “Deep-learning-based device fingerprinting for increased lora-iot security: Sensitivity to network deployment changes,” _IEEE Network_ , vol. 36, no. 3, pp. 204–210, 2022.
  * [7] F. Restuccia, S. D’Oro, A. Al-Shawabka, M. Belgiovine, L. Angioloni, S. Ioannidis, K. Chowdhury, and T. Melodia, “Deepradioid: Real-time channel-resilient optimization of deep learning-based radio fingerprinting algorithms,” in _Proc. of the 20th ACM International Symposium on Mobile Ad Hoc Networking and Computing_ , 2019, pp. 51–60.
  * [8] S. Rajendran and Z. Sun, “Rf impairment model-based iot physical-layer identification for enhanced domain generalization,” _IEEE Transactions on Information Forensics and Security_ , vol. 17, pp. 1285–1299, 2022.
  * [9] A. Al-Shawabka, P. Pietraski, S. B. Pattar, F. Restuccia, and T. Melodia, “Deeplora: Fingerprinting lora devices at scale through deep learning and data augmentation,” in _Proceedings of the Twenty-second International Symposium on Theory, Algorithmic Foundations, and Protocol Design for Mobile Networks and Mobile Computing_ , 2021, pp. 251–260.
  * [10] I. Redko, E. Morvant, A. Habrard, M. Sebban, and Y. Bennani, “A survey on domain adaptation theory: learning bounds and theoretical guarantees,” _arXiv preprint arXiv:2004.11829_ , 2020.
  * [11] J. Gaskin, B. Hamdaoui, and W.-K. Wong, “Tweak: Towards portable deep learning models for domain-agnostic LoRa device authentication,” in _2022 IEEE conference on communications and network security (CNS)_. IEEE, 2022.
  * [12] T. Wang, L. Yu, W. Wang, and J. Wei, “Specific emitter identification based on the multi-discrepancy deep adaptation network,” _IET Radar, Sonar & Navigation_, 2022.
  * [13] J. Wang, B. Zhang, J. Zhang, N. Yang, G. Wei, and D. Guo, “Specific emitter identification based on deep adversarial domain adaptation,” in _2021 4th International Conference on Information Communication and Signal Processing (ICICSP)_. IEEE, 2021, pp. 104–109.
  * [14] I. Higgins, D. Amos, D. Pfau, S. Racaniere, L. Matthey, D. Rezende, and A. Lerchner, “Towards a definition of disentangled representations,” _arXiv preprint arXiv:1812.02230_ , 2018.
  * [15] K. Bousmalis, G. Trigeorgis, N. Silberman, D. Krishnan, and D. Erhan, “Domain separation networks,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [16] A. Elmaghbub and B. Hamdaoui, “Comprehensive rf dataset collection and release: A deep learning-based device fingerprinting use case,” in _2021 IEEE Globecom Workshops (GC Wkshps)_. IEEE, 2021, pp. 1–7.
  * [17] A. Al-Shawabka, F. Restuccia, S. D’Oro, T. Jian, B. C. Rendon, N. Soltani, J. Dy, S. Ioannidis, K. Chowdhury, and T. Melodia, “Exposing the fingerprint: Dissecting the impact of the wireless channel on radio fingerprinting,” in _IEEE INFOCOM 2020-IEEE Conference on Computer Communications_. IEEE, 2020, pp. 646–655.
  * [18] Y. Ganin, E. Ustinova, H. Ajakan, P. Germain, H. Larochelle, F. Laviolette, M. Marchand, and V. Lempitsky, “Domain-adversarial training of neural networks,” _The journal of machine learning research_ , vol. 17, no. 1, pp. 2096–2030, 2016.
  * [19] D. Eigen, C. Puhrsch, and R. Fergus, “Depth map prediction from a single image using a multi-scale deep network,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [20] M. Wulfmeier, A. Bewley, and I. Posner, “Addressing appearance change in outdoor robotics with adversarial domain adaptation,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 1551–1558.


Also at ]Office of Theoretical and Applied
Mechanics, Academy of Athens, 10679, Greece

# A Rheological Analogue for Brownian Motion with Hydrodynamic Memory

Nicos Makris<EMAIL_ADDRESS>[  Dept. of Civil and Environmental Engineering,
Southern Methodist University, Dallas, Texas, 75276

###### Abstract

When the density of the fluid surrounding suspended Brownian particles is
appreciable, in addition to the forces appearing in the traditional Ornstein
and Uhlenbeck theory of Brownian motion, additional forces emerge as the
displaced fluid in the vicinity of the randomly moving Brownian particle acts
back on the particle giving rise to long-range force correlations which
manifest as a “long-time tail” in the decay of the velocity autocorrelation
function known as hydrodynamic memory. In this paper, after recognizing that
for Brownian particles immersed in a Newtonian, viscous fluid, the
hydrodynamic memory term in the generalized Langevin equation is essentially
the $\nicefrac{{\text{1}}}{{\text{2}}}$ fractional derivative of the velocity
of the Brownian particle, we present a rheological analogue for Brownian
motion with hydrodynamic memory which consists of a linear dashpot of a
fractional Scott-Blair element and an inerter. The synthesis of the proposed
mechanical network that is suggested from the structure of the generalized
Langevin equation simplifies appreciably the calculations of the mean-square
displacement and its time-derivatives which can also be expressed in terms of
the two-parameter Mittag–Leffler function.

Brownian Motion; Random process; Unsteady motion; Viscoelasticity;
Hydrodynamic Memory; Time-Response Functions

##  I Introduction

Following Einstein’s [1] ‘long-term’ theory of Brownian motion, Langevin [2]
presented a force–balance dynamic approach where the inertia of a Brownian
microsphere with mass $m$ and radius $R$ balances the friction force that
develops as the particle moves inside the surrounding fluid, together with the
random excitation forces $f_{R}(t)$, that originate from the collisions of the
fluid molecules on the Brownian microsphere

$m\frac{\mathrm{d}v(t)}{\mathrm{d}t}=-\zeta v(t)+f_{R}(t)$ (1)

In Eq. (1), $v(t)=\frac{\mathrm{d}r(t)}{\mathrm{d}t}$ is the particle velocity
and $\zeta v(t)$ is a viscous drag force proportional to the velocity of the
Brownian particle. For a steady motion in a memoryless, viscous fluid with
viscosity $\eta$, the drag coefficient is given by Stokes’ law
$\zeta=\text{6}\pi R\eta$ [3] which has been traditionally accepted when long-
range correlations are not present [4, 5]. According to the Ornstein and
Uhlenbeck theory of Brownian motion [6, 7], the random excitation $f_{R}(t)$
has a zero average value over time $\left\langle f_{R}(t)\right\rangle=$ 0;
while for the memoryless, Newtonian viscous fluid that only dissipates energy,
the correlation function contracts to a Dirac delta function [8]

$\left\langle
f_{R}(t_{\text{1}})f_{R}(t_{\text{2}})\right\rangle=A\delta(t_{\text{1}}-t_{\text{2}})$
(2)

with $A$ being a constant that expresses the strength of the random forces.
Equation (2) indicates that successive collisions of the fluid molecules on
the Brownian particle are uncorrelated.

Given the random nature of the excitation force $f_{R}(t)$, the Langevin Eq.
(1) can be integrated in terms of ensemble averages in association with Eq.
(2) and the velocity autocorrelation function of Brownian particles immersed
in a memoryless viscous, Newtonian fluid with viscosity $\eta$ is [6, 7]

$\displaystyle\left\langle v(\text{0})v(t)\right\rangle$
$\displaystyle=\left\langle v(\xi)v(\xi+t)\right\rangle$ (3)
$\displaystyle=\lim_{T\rightarrow\infty}\frac{\text{1}}{T}\int_{\text{0}}^{T}v(\xi)v(\xi+t)\,\mathrm{d}\xi=\frac{NK_{B}T}{m}e^{-\nicefrac{{t}}{{\tau}}}$

where $N\in\left\\{\text{1}\text{,}\text{2}\text{,}\text{3}\right\\}$ is the
number of spatial dimensions, $K_{B}$ is Boltzmann’s constant, $T$ is the
equilibrium temperature and $\tau=\frac{m}{\text{6}\pi R\eta}$ is the
dissipation time-scale of the perpetual fluctuation–dissipation process.

When the density of the fluid surrounding the Brownian particle is
appreciable, there is a convoluted interplay between the unsteady motion of
the particle and the motion of the displaced fluid. Accordingly, in addition
to the forces appearing in Eq. (1), the Brownian particle also displaces the
fluid in its immediate vicinity and in return the fluid acts back on the
particle and gives rise to long-range correlations which are different than
the delta-correlations of Eq. (2). The interaction of the Brownian particle
with the displaced fluid can be described by the addition of inertia and
memory terms to the basic Langevin Eq. (1) [3, 9, 10, 11, 12, 13].

$\displaystyle m\frac{\mathrm{d}v(t)}{\mathrm{d}t}$
$\displaystyle=-\frac{\text{2}}{\text{3}}\pi
R^{\text{3}}\rho_{f}\frac{\mathrm{d}v(t)}{\mathrm{d}t}-\text{6}\pi R\eta v(t)$
(4)
$\displaystyle-\text{6}R^{\text{2}}\sqrt{\pi\rho_{f}\eta}\int_{\text{0}^{-}}^{t}\frac{\nicefrac{{\mathrm{d}v(\xi)}}{{\mathrm{d}\xi}}}{(t-\xi)^{\nicefrac{{\text{1}}}{{\text{2}}}}}\,\mathrm{d}\xi+f_{R}(t)$

The quantity $\frac{\text{2}}{\text{3}}\pi R^{\text{3}}\rho_{f}$ in the
additional inertia term of Eq. (4) is half the mass of the displaced fluid
$\frac{m_{f}}{\text{2}}$ $\bigg{(}m_{f}=\frac{\text{4}}{\text{3}}\pi
R^{\text{3}}\rho_{f}$ with $\rho_{f}=$ density of the surrounding
fluid$\bigg{)}$; whereas, the following term $\text{6}\pi R\eta v(t)$ is the
ordinary Stokes’ friction force, also present in Eq. (1). The convolution
integral term in Eq. (4) is a hydrodynamic memory force that emerges from the
reaction of the displaced fluid on the moving Brownian particle [3].

By recalling that the Riemann–Liouville fractional integral of a continuous
function $f(t)$ is defined as [14, 15, 16, 17]

$I^{q}\left[f(t)\right]=\frac{\text{1}}{\Gamma(q)}\int_{\text{0}^{-}}^{t}\frac{f(\xi)}{(t-\xi)^{\text{1}-q}}\,\mathrm{d}\xi\text{,}\quad
q\in\mathbb{R}^{+}$ (5)

the hydrodynamic memory integral in Eq. (4) is essentially the fractional
integral of order $q=$ $\nicefrac{{\text{1}}}{{\text{2}}}$ of the acceleration
history $\frac{\mathrm{d}v(t)}{\mathrm{d}t}$ of the Brownian particle.
Accordingly, by virtue of Eq. (5), Eq. (4) is expressed as

$M\frac{\mathrm{d}v(t)}{\mathrm{d}t}+\text{6}\pi R\eta
v(t)+\text{6}R^{\text{2}}\sqrt{\pi\rho_{f}\eta}\,\Gamma\left(\nicefrac{{\text{1}}}{{\text{2}}}\right)I^{\nicefrac{{\text{1}}}{{\text{2}}}}\left[\frac{\mathrm{d}v(t)}{\mathrm{d}t}\right]=f_{R}(t)$
(6)

where $M=m+\frac{\text{1}}{\text{2}}m_{f}=\frac{\text{4}}{\text{3}}\pi
R^{\text{3}}\left(\rho_{p}+\frac{\text{1}}{\text{2}}\rho_{f}\right)$ with
$\rho_{p}=$ mass density of the Brownian particle. By using the properties of
fractional calculus [14, 15, 16, 17] the fractional integral of order
$\nicefrac{{\text{1}}}{{\text{2}}}$ of the derivative of the velocity (first-
order derivative) is the fractional derivative of order
$\nicefrac{{\text{1}}}{{\text{2}}}$ of the velocity history

$I^{\nicefrac{{\text{1}}}{{\text{2}}}}\left[\frac{\mathrm{d}v(t)}{\mathrm{d}t}\right]=\frac{\mathrm{d}^{\nicefrac{{\text{1}}}{{\text{2}}}}v(t)}{\mathrm{d}t^{\nicefrac{{\text{1}}}{{\text{2}}}}}=\frac{\text{1}}{\Gamma\left(-\frac{\text{1}}{\text{2}}\right)}\int_{\text{0}^{-}}^{\text{1}}\frac{v(\xi)}{(t-\xi)^{\nicefrac{{\text{3}}}{{\text{2}}}}}\,\mathrm{d}\xi$
(7)

More generally, the fractional derivative of order $q\in\mathbb{R}^{+}$ of a
continuous function, $f(t)$, is defined within the context of generalized
functions as the convolution of $f(t)$ with the kernel
$\frac{\text{1}}{\Gamma(-q)}\frac{\text{1}}{t^{q+\text{1}}}$, whereas, the
Laplace transform of the fractional derivative of $f(t)$ is
$\mathcal{L}\left\\{\frac{\mathrm{d}^{q}f(t)}{\mathrm{d}t^{q}}\right\\}=$
$\int_{\text{0}^{-}}^{\infty}\frac{\mathrm{d}^{q}f(t)}{\mathrm{d}t^{q}}e^{-st}\,\mathrm{d}t=s^{q}\mathcal{L}\left\\{f(t)\right\\}=s^{q}\mathcal{F}(s)$
[14, 16, 18, 19]. By employing the result of Eq. (7) in association with
$\Gamma\left(\frac{\text{1}}{\text{2}}\right)=\sqrt{\pi}$, Eq. (6) simplifies
to

$M\frac{\mathrm{d}v(t)}{\mathrm{d}t}+\text{6}\pi
R^{\text{2}}\sqrt{\rho_{f}\eta}\frac{\mathrm{d}^{\nicefrac{{\text{1}}}{{\text{2}}}}v(t)}{\mathrm{d}t^{\nicefrac{{\text{1}}}{{\text{2}}}}}+\text{6}\pi
R\eta v(t)=f_{R}(t)$ (8)

Equation (8) offers the remarkable result that fractional differentials do not
only appear when modeling the phenomenological power-law relaxation of a wide
range of viscoelastic materials [20, 21, 22, 23], but also they emerge
naturally from the solution of continuum mechanics equations as they result
from fundamental conservation laws [3, 9, 10].

The Langevin Eq. (4) or (8) which accounts for the hydrodynamic memory was
solved analytically by Widom [10] after solving an integral equation for the
velocity autocorrelation function of the Brownian particle in association with
the appropriate long-range correlation of the random process. Building on
Widom’s (1971) solution the following result for the velocity autocorrelation
function can be obtained [13, 24]

$\displaystyle\left\langle
v(\text{0})v(t)\right\rangle=\frac{NK_{B}T}{M(b-a)}$
$\displaystyle\left[be^{b^{\text{2}}t}\text{Erfc}\left(b\sqrt{t}\right)\right.$
(9)
$\displaystyle\left.-ae^{a^{\text{2}}t}\text{Erfc}\left(a\sqrt{t}\right)\right]$

with

$a=\frac{z+\sqrt{z^{\text{2}}-\text{4}\zeta M}}{\text{2}M}\quad\text{and}\quad
b=\frac{z-\sqrt{z^{\text{2}}-\text{4}\zeta M}}{\text{2}M}$ (10)

In this paper we follow the notation of Clercx and Schram [13] where
$z=\text{6}\pi R^{\text{2}}\sqrt{\rho_{f}\eta}$ with units
[M][L][T]${}^{-\nicefrac{{\text{3}}}{{\text{2}}}}$ (say Pa
s${}^{\nicefrac{{\text{3}}}{{\text{2}}}}$) — that is the coefficient of the
$\nicefrac{{\text{1}}}{{\text{2}}}$ fractional derivative of the velocity of
the Brownian particle appearing in the Langevin eq. (8).

##  II The Dashpot–Inertpot–Inerter Parallel Connection

For the analysis of classical Brownian motion, as was treated by [6] and [7]
(when effects of the hydrodynamic memory are neglected) of microparticles
suspended in any linear, isotropic viscoelastic material, Makris [25]
presented a viscous–viscoelastic correspondence principle which reveals that
the mean-square displacement and its time derivatives [26] of the randomly
moving Brownian microspheres with mass $m$ and radius $R$ are identical to the
deterministic time response functions (creep compliance
$J(t)=\frac{\text{3}\pi R}{NK_{B}T}\left\langle\Delta
r^{\text{2}}(t)\right\rangle$, impulse response function
$h(t)=\frac{\text{3}\pi R}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}$ and impulse strain–rate response
function $\psi(t)=\frac{\text{3}\pi
R}{NK_{B}T}\frac{\mathrm{d}^{\text{2}}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t^{\text{2}}}=\frac{\text{6}\pi
R}{NK_{B}T}\left\langle v(\text{0})v(t)\right\rangle$) of a mechanical network
that is a parallel connection of the linear viscoelastic material (within
which the microspheres are immersed) with an inerter with distributed
inertance $m_{R}=\frac{m}{\text{6}\pi R}$.

In this paper we build upon the viscous–viscoelastic correspondence principle
for Brownian motion [25, 26] to develop a mechanical analogue for Brownian
motion with hydrodynamic memory which simplifies appreciably the calculations
of the velocity autocorrelation function and the mean-square displacement. The
synthesis of the proposed macroscopic mechanical network is suggested from the
left-hand side of the Langevin Eq. (8), which consists of an inertia term, a
$\nicefrac{{\text{1}}}{{\text{2}}}$-order fractional derivative of the
velocity term and a viscous term acting in parallel to balance the random
force $f_{R}(t)$. Accordingly, we examine the frequency and the time-response
functions of the mechanical network shown in Fig. 1 which is a parallel
connection of a dashpot with shear viscosity $\eta$, a fractional Scott–Blair
element with material constant
$\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}=R\sqrt{\rho_{f}\eta}=\frac{z}{\text{6}\pi
R}$ and an inerter with distributed inertance $m_{R}=\frac{M}{\text{6}\pi
R}=\frac{\text{2}}{\text{9}}R^{\text{2}}\rho_{f}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)$.

Figure 1: The rheological analogue for Brownian motion in a Newtonian viscous
fluid with hydrodynamic memory. A dashpot with viscosity $\eta$ is connected
in parallel with a Scott-Blair element with material constant
$\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}=R\sqrt{\rho_{f}\eta}$ and an inerter
with distributed inertance $m_{R}=\frac{M}{\text{6}\pi
R}=\frac{\text{2}}{\text{9}}R^{\text{2}}\rho_{f}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)$
(the dashpot–inertpot–inerter parallel connection).

Following Nutting’s [20] observations that the stress response of several
fluid-like materials when subjected to a step–strain
$(\gamma(t)=U(t-\text{0}))$ excitation decays following a power law:
$\tau(t)=G_{\text{ve}}(t)\sim t^{-q}$ with 0 $\leqslant q\leqslant$ 1 and the
early work of Gemant [21, 27] on fractional differentials, Scott-Blair [22,
23], proposed the springpot element (0 $\leqslant q\leqslant$ 1) which is a
mechanical idealization in–between an elastic spring and a viscous dashpot
with constitutive law

$\tau(t)=\mu_{q}\frac{\mathrm{d}^{q}\gamma(t)}{\mathrm{d}t^{q}}\text{,}\quad
q\in\mathbb{R}^{+}$ (11)

In this study we adopt the name ‘Scott–Blair element’ compared to the
restrictive ‘springpot’ element since the order of differentiation of the
fractional element appearing in Fig. 1 is
$q=\nicefrac{{\text{3}}}{{\text{2}}}$. In view that 1 $\leqslant
q=\nicefrac{{\text{3}}}{{\text{2}}}\leqslant$ 2 and by adopting the
terminology introduced by Scott–Blair (‘springpot’), the fractional element
appearing in Fig. 1 is coined the ‘inertpot’ since it is a mechanical
idealization in–between an inerter ($q=$ 2) and a viscous dashpot ($q=$ 1).

An inerter is a linear mechanical element for which, at the force–displacement
level, the output force is proportional only to the relative acceleration of
its end-nodes (terminals) and its constant of proportionality is the inertance
with units of mass [M] [28, 29, 30, 31]. At the stress–strain level, the
constant of proportionality of the inerter shown schematically in Fig. 1 is
the distributed inertance, $m_{R}$, with units [M][L]${}^{-\text{1}}$ (i.e Pa
s${}^{\text{2}}$). Given the parallel connection of the
dashpot–inertpot–inerter network shown in Fig. 1, its constitutive law is

$\tau(t)=\eta\frac{\mathrm{d}\gamma(t)}{\mathrm{d}t}+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}\frac{\mathrm{d}^{\nicefrac{{\text{3}}}{{\text{2}}}}\gamma(t)}{\mathrm{d}t^{\nicefrac{{\text{3}}}{{\text{2}}}}}+m_{R}\frac{\mathrm{d}^{\text{2}}\gamma(t)}{\mathrm{d}t^{\text{2}}}$
(12)

The Laplace transform of Eq. (12) gives

$\tau(s)=\mathcal{G}(s)\gamma(s)=\left(\eta
s+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+m_{R}s^{\text{2}}\right)\gamma(s)$
(13)

where $\mathcal{G}(s)=\left(\eta
s+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+m_{R}s^{\text{2}}\right)$
is the complex dynamic modulus. The complex dynamic compliance of the
rheological network shown in Fig. 1 is

$\displaystyle\mathcal{J}(s)=\frac{\text{1}}{\mathcal{G}(s)}$
$\displaystyle=\frac{\text{1}}{\eta
s+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+m_{R}s^{\text{2}}}$
(14)
$\displaystyle=\frac{\text{1}}{m_{R}}\frac{\text{1}}{s}\frac{\text{1}}{\left(s+\frac{\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}}{m_{R}}s^{\nicefrac{{\text{1}}}{{\text{2}}}}+\frac{\eta}{m_{R}}\right)}$

The recently published viscous–viscoelastic correspondence principle for
classical Brownian motion (without hydrodynamic memory) [25] uncovered that
ensemble averages of the motion of Brownian particles are proportional to
macroscopic viscoelastic response functions. For instance the Laplace
transform of the velocity autocorrelation function of Brownian particles
satisfies the identity

$\left\langle v(s)v(\text{0})\right\rangle=\frac{NK_{B}T}{\text{6}\pi
R}s\mathcal{J}(s)$ (15)

where
$s\mathcal{J}(s)=\phi(s)=\frac{\accentset{\mbox{\large.}}{\gamma}(s)}{\tau(s)}$
is the complex dynamic fluidity (admittance) of the rheological network that
resists the random forces, $f_{R}(t)$. In this paper we extend this concept
for the case with hydrodynamic memory and is implemented to the Langevin Eq.
(8). Accordingly, substitution of Eq. (14) into Eq. (15) gives

$\displaystyle\left\langle v(s)v(\text{0})\right\rangle$
$\displaystyle=\frac{NK_{B}T}{\text{6}\pi
R}\frac{\text{1}}{m_{R}}\frac{\text{1}}{s+\frac{\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}}{m_{R}}s^{\nicefrac{{\text{1}}}{{\text{2}}}}+\frac{\eta}{m_{R}}}$
(16)
$\displaystyle=\frac{NK_{B}T}{M}\frac{\text{1}}{s+\frac{z}{M}s^{\nicefrac{{\text{1}}}{{\text{2}}}}+\frac{\zeta}{M}}$

where $m_{R}=\frac{M}{\text{6}\pi
R}=\frac{\text{2}}{\text{9}}R^{\text{2}}\rho_{f}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)$,
$\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}=R\sqrt{\rho_{f}\eta}=\frac{z}{\text{6}\pi
R}$ and the last expression of Eq. (16) follows the Clercx and Schram [13]
notation. The inverse Laplace transform of Eq. (16) is evaluated by employing
the fraction expansion

$\frac{\text{1}}{s+(a+b)s^{\nicefrac{{\text{1}}}{{\text{2}}}}+ab}=\frac{\text{1}}{b-a}\left(\frac{\text{1}}{s^{\nicefrac{{\text{1}}}{{\text{2}}}}+a}-\frac{\text{1}}{s^{\nicefrac{{\text{1}}}{{\text{2}}}}+b}\right)$
(17)

together with the result [32]

$\displaystyle\mathcal{L}^{-\text{1}}\left\\{\frac{\text{1}}{s^{\nicefrac{{\text{1}}}{{\text{2}}}}+a}\right\\}$
$\displaystyle=\frac{\text{1}}{\sqrt{\pi
t}}-ae^{a^{\text{2}}t}\text{Erfc}(a\sqrt{t})$ (18)
$\displaystyle=\frac{\text{1}}{\sqrt{t}}E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{1}}}{{\text{2}}}}(-a\sqrt{t})$

where Erfc$(x)$ is the complementary error function and $E_{p,q}(x)$ is the
two parameter Mittag–Leffler function [33, 34]

$E_{p,q}(x)=\sum_{j=\text{0}}^{\infty}\frac{x^{j}}{\Gamma(jp+q)}\text{,}\quad
p,q>\text{0}$ (19)

With reference to the denominators of the left-hand side of Eqs. (17) and the
right-hand side of Eq. (16), the solution of the system of equations
$a+b=\frac{z}{M}$ and $ab=\frac{\zeta}{M}$ yields the values of $a$ and $b$
offered by Eq. (10).

By virtue of Eqs. (17) and (18), the inverse Laplace transform of Eq. (16)
yields the result offered by Eq. (9). Consequently, the analysis presented
herein by implementing the viscous–viscoelastic correspondence principle [25,
26] confirms that the velocity autocorrelation function of Brownian particles
suspended in a viscous Newtonian fluid where effects of hydrodynamic memory
are accounted is given by:

$\left\langle v(\text{0})v(t)\right\rangle=\frac{NK_{B}T}{\text{6}\pi
R}\psi(t)$ (20)

where
$\psi(t)=\mathcal{L}^{-\text{1}}\left\\{\phi(s)\right\\}=\mathcal{L}^{-\text{1}}\left\\{s\mathcal{J}(s)\right\\}$
is the impulse strain-rate response function of the mechanical network shown
in Fig. 1 (the dashpot–inertpot–inerter parallel connection), defined as the
resulting strain-rate output, $\accentset{\mbox{\large.}}{\gamma}(t)$, due to
an impulsive stress input, $\tau(t)=\delta(t-\text{0})$.

In terms of the two-parameter Mittag–Leffler function given by the identity
(18), an alternative expression to Eq. (9) for the velocity autocorrelation
function of Brownian particles immersed in a viscous, Newtonian fluid with
hydrodynamic memory is

$\displaystyle\left\langle v(\text{0})v(t)\right\rangle$
$\displaystyle=\frac{\text{1}}{\text{2}}\frac{\mathrm{d}^{\text{2}}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t^{\text{2}}}=\frac{NK_{B}T}{M(b-a)}\frac{\text{1}}{\sqrt{t}}$
(21)
$\displaystyle\times\left[E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{1}}}{{\text{2}}}}\left(-a\sqrt{t}\right)-E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{1}}}{{\text{2}}}}\left(-b\sqrt{t}\right)\right]$

The constants $a$ and $b$ offered by Eq. (10) involve the radical

$\displaystyle\sqrt{z^{\text{2}}-\text{4}\zeta M}$
$\displaystyle=z\sqrt{\text{1}-\frac{\text{8}}{\text{9}}\left(\frac{1}{2}+\frac{\rho_{p}}{\rho_{f}}\right)}$
(22) $\displaystyle=\text{6}\pi
R^{\text{2}}\sqrt{\rho_{f}\eta}\sqrt{\text{1}-\frac{\text{8}}{\text{9}}\left(\frac{1}{2}+\frac{\rho_{p}}{\rho_{f}}\right)}$

which remains real as long as $\rho_{p}<\frac{5}{8}\rho_{f}$. For most
practical cases $\rho_{p}>\rho_{f}$. As an example, when melanine rasin
microspheres ($\rho_{p}=\text{1570}$ kg/m${}^{\text{3}}$) are immensed in
water ($\rho_{f}=\text{1000}$ kg/m${}^{\text{3}}$),
$\rho_{p}=\text{1.57}\rho_{f}$ [35, 36]; whereas, when they are immensed in
acetone ($\rho_{f}=\text{784}$ kg/m${}^{\text{3}}$),
$\rho_{p}=\text{2}\rho_{f}$. Accordingly, for most practical cases, the
radical $\sqrt{z^{2}-\text{4}\zeta M}$ is an imaginary quantity

$\sqrt{z^{\text{2}}-\text{4}\zeta
M}=\operatorname{i}z\sqrt{\frac{\text{8}}{\text{9}}\left(\frac{1}{2}+\frac{\rho_{p}}{\rho_{f}}\right)-\text{1}}=\operatorname{i}z\sqrt{\sigma-\text{1}}$
(23)

where $\operatorname{i}=\sqrt{-\text{1}}$ and

$\sigma=\frac{\text{8}}{\text{9}}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)$
(24)

is a parameter of the Brownian particle–Newtonian fluid system described by
the Langevin Eq. (4) of (8). With the introduction of the parameter $\sigma$
defined by Eq. (24), the constants $a$ and $b$ introduced by [13] and offered
by Eq. (10) are expressed as
$a=\frac{z}{\text{2}M}(\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}})$ and
$b=\frac{z}{\text{2}M}(\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}})$,
whereas, the factor

$\frac{z}{\text{2}M}=\frac{\text{9}}{\text{4}}\frac{\text{1}}{R}\frac{\sqrt{\eta}}{\sqrt{\rho_{f}}}\frac{\text{1}}{\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)}=\frac{\text{2}}{\sigma\sqrt{\tau_{f}}}$
(25)

where $\sigma$ is defined by Eq. (24) and
$\tau_{f}=\frac{R^{\text{2}}\rho_{f}}{\eta}$ is a time-scale which describes
the time needed by the perturbed fluid flow to diffuse over one Brownian
microsphere radius [12, 35, 37]. Accordingly, for
$\sigma>\text{1}\left(\frac{\rho_{p}}{\rho_{f}}>\frac{\text{5}}{\text{8}}\right)$
the constants $a$ and $b$ introduced by Eq. (10) assume the expressions

$a=\frac{\text{2}}{\sigma\sqrt{\tau_{f}}}\left(\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}}\right)\enskip\text{and}\enskip
b=\frac{\text{2}}{\sigma\sqrt{\tau_{f}}}\left(\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}}\right)$
(26)

and the expression for the velocity autocorrelation function offered by Eq.
(9) gives

Figure 2: Normalized velocity autocorrelation function of Brownian particles
immersed in a Newtonian, viscous fluid with viscosity $\eta$; when effects
from hydrodynamic memory are accounted. The curves for
$\sigma=\frac{\text{8}}{\text{9}}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)=$
0.9 and 0.99 are computed with Eq. (29), whereas the curves for $\sigma=$
1.01, 1.5, 2 and 2.5 are computed with Eq. (27) showing the continuity of the
solution as $\sigma$ crosses the value of $\sigma=$ 1.

$\displaystyle\left\langle v(\text{0})v(t)\right\rangle$
$\displaystyle=\frac{NK_{B}T}{M}\frac{e^{-\frac{\text{4}(\sigma-\text{2})}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}}{\text{2}\sqrt{\sigma-\text{1}}}\Bigg{[}\left(\sqrt{\sigma-\text{1}}+\operatorname{i}\right)e^{-\frac{\operatorname{i}\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)$
(27)
$\displaystyle+\left(\sqrt{\sigma-\text{1}}-\operatorname{i}\right)e^{+\frac{\operatorname{i}\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\Bigg{]}\text{,}\quad\sigma>\text{1}$

The velocity autocorrelation function offered by Eq. (27) for
$\sigma>\text{1}$ is a real-valued function given that the imaginary parts of
the right-hand side of Eq. (27) cancel. For the case when
$\sigma<\text{1}\left(\frac{\rho_{p}}{\rho_{f}}<\frac{\text{5}}{\text{8}}\right)$,
the parameters $a$ and $b$ introduced by [13] assume the expressions

$a=\frac{\text{2}}{\sigma\sqrt{\tau_{f}}}\left(\text{1}+\sqrt{\text{1}-\sigma}\right)\quad\text{and}\quad
b=\frac{\text{2}}{\sigma\sqrt{\tau_{f}}}\left(\text{1}-\sqrt{\text{1}-\sigma}\right)$
(28)

and the expession for the velocity autocorrelation function given by Eq. (9)
yields

$\displaystyle\left\langle v(\text{0})v(t)\right\rangle$
$\displaystyle=\frac{NK_{B}T}{M}\frac{e^{\frac{\text{4}(\text{2}-\sigma)}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}}{\text{2}\sqrt{\text{1}-\sigma}}\Bigg{[}-\left(\text{1}-\sqrt{\text{1}-\sigma}\right)e^{-\frac{\text{8}\sqrt{\text{1}-\sigma}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}-\sqrt{\text{1}-\sigma}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)$
(29)
$\displaystyle+\left(\text{1}+\sqrt{\text{1}-\sigma}\right)e^{+\frac{\text{8}\sqrt{\text{1}-\sigma}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}+\sqrt{\text{1}+\sigma}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\Bigg{]}\text{,}\quad\sigma<\text{1}$

Figure 2 plots the normalized velocity autocorrelation function
$\frac{M}{NK_{B}T}\left\langle v(\text{0})v(t)\right\rangle=m_{R}\psi(t)$ for
various values of the parameter $\sigma$. For $\sigma=\text{0.9}$ and 0.99 the
cuves are computed with Eq. (29); whereas, for values of $\sigma=\text{1.01}$,
1.5, 2 and 2.5 the curves are computed with Eq. (27).

##  III Time Derivative of the Mean-Square Displacement of Brownian Particles
in a Viscous Fluid with Hydrodynamic Memory

The time derivative of the mean square displacement of Brownian particles has
been occasionally interpreted as a time-dependent diffusion coefficient
$\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}=\text{2}ND(t)$ [12, 38, 39, 40, 41]
and has a central role in the interpretation of Brownian motion at various
time scales [42] and confined spacings [43].

With reference to Eq. (15) in association that the velocity autocorrelation
function $\left\langle
v(\text{0})v(t)\right\rangle=\frac{\text{1}}{\text{2}}\frac{\mathrm{d}^{\text{2}}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t^{\text{2}}}$ is the second time-
derivative of the mean-square displacement [44, 45]; in a recent paper Makris
[26] uncovered that the complex dynamic compliance, $\mathcal{J}(s)$, of the
rheological analogue for Brownian motion, such as the one shown in Fig. 1 is
proportional to the Laplace transform of the time-derivative of the mean-
square displacement

$\mathcal{J}(s)=\frac{\text{3}\pi
R}{NK_{B}T}\mathcal{L}\left\\{\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}\right\\}$ (30)

The inverse Laplace transform of the complex dynamic compliance,
$\mathcal{J}(s)$ given by Eq. (14) is the impulse fluidity,
$\phi(t)=\mathcal{L}^{-\text{1}}\left\\{\mathcal{J}(s)\right\\}$ [46, 47, 48],
defined as the resulting strain $\gamma(t)$, due to an impulse stress input,
$\tau(t)=\delta(t-\text{0})$. The impulse function $\delta(t-\text{0})$ is the
Dirac delta function [8] with the property
$\mathcal{L}\left\\{\delta(t-\xi)\right\\}=\int_{\text{0}^{-}}^{\infty}\delta(t-\xi)e^{-st}\,\mathrm{d}t=e^{-\xi
s}$. The equivalent of the impulse fluidity, $\phi(t)$, at the
displacement–force level is the impulse response function often expressed with
$h(t)$. Given that the term “impulse response function” (rather than the term
“impulse fluidity”) is widely known and used in dynamics [49], structural
mechanics [50], electrical signal processing [51, 52] and economics [53],
Makris [26] adopted the term “impulse response function $=h(t)$”, rather than
the term impulse fluidity” used narrowly in the viscoelasticity literature
alone. Accordingly, from Eq. (14) the impulse response function of the
dashpot–interpot–inerter parallel connection is

$h(t)=\mathcal{L}^{-\text{1}}\left\\{\mathcal{J}(s)\right\\}=\frac{\text{1}}{m_{R}}\mathcal{L}^{-\text{1}}\left\\{\frac{\text{1}}{s}\frac{\text{1}}{\left(s+\frac{\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}}{m_{R}}s^{\nicefrac{{\text{1}}}{{\text{2}}}}+\frac{\eta}{m_{R}}\right)}\right\\}$
(31)

where again $m_{R}=\frac{M}{\text{6}\pi
R}=\frac{\text{2}}{\text{9}}R^{\text{2}}\rho_{f}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)$
and
$\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}=R\sqrt{\rho_{f}\eta}=\frac{z}{\text{6}\pi
R}$ with $z$ being the Clercx and Schram [13] parameter appearing in Eq. (10).
By following the Clercx and Schram [13] notation in association with Eq. (30),
the impulse response function given by Eq. (31) is

$\displaystyle h(t)$ $\displaystyle=\frac{\text{3}\pi
R}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}$ (32)
$\displaystyle=\frac{\text{6}\pi
R}{M}\mathcal{L}^{-\text{1}}\left\\{\frac{\text{1}}{s}\frac{\text{1}}{\left(s+\frac{z}{M}s^{\nicefrac{{\text{1}}}{{\text{2}}}}+\frac{\zeta}{M}\right)}\right\\}$

The inverse Laplace transform of Eq. (31) is evaluated by employing the
fraction expansion

$\displaystyle\frac{\text{1}}{s\left[s+(a+b)s^{\nicefrac{{\text{1}}}{{\text{2}}}}+ab\right]}$
(33)
$\displaystyle=\frac{\text{1}}{b-a}\left[\frac{\text{1}}{s(s^{\nicefrac{{\text{1}}}{{\text{2}}}}+a)}-\frac{\text{1}}{s(s^{\nicefrac{{\text{1}}}{{\text{2}}}}+b)}\right]$

together with the result [32]

$\displaystyle\mathcal{L}^{-\text{1}}\left\\{\frac{\text{1}}{s(s^{\nicefrac{{\text{1}}}{{\text{2}}}}+a)}\right\\}$
$\displaystyle=\frac{\text{1}}{a}\left[\text{1}-e^{a^{\text{2}}t}\text{Erfc}(a\sqrt{t})\right]$
(34)
$\displaystyle=\sqrt{t}E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{3}}}{{\text{2}}}}(-a\sqrt{t})$

With reference to the denominators of the left-hand side of Eq. (33) and the
right-hand side of Eq. (32), the solution of the system of equations
$a+b=\frac{z}{M}$ and $ab=\frac{\zeta}{M}$ yields the values of $a$ and $b$
offered by Eq. (10). Accordingly, by virtue of Eqs. (33) and (34), the impulse
response function for Brownian motion in a Newtonian viscous fluid with
hydrodynamic memory is

$h(t)=\frac{\text{3}\pi R}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}=\frac{\text{6}\pi
R}{M}\frac{\text{1}}{b-a}\left\\{\frac{\text{1}}{a}\left[\text{1}-e^{a^{\text{2}}t}\text{Erfc}(a\sqrt{t})\right]-\frac{\text{1}}{b}\left[\text{1}-e^{b^{\text{2}}t}\text{Erfc}(b\sqrt{t})\right]\right\\}$
(35)

When $\frac{\rho_{p}}{\rho_{f}}>\frac{\text{5}}{\text{8}}$, parameter
$\sigma>\text{1}$ and the parameters $a$ and $b$ appearing in Eq. (35) are
complex-valued and are offered by Eq. (26). In this case ($\sigma>\text{1}$)
Eq. (35) becomes

$\displaystyle h(t)=\frac{\text{3}\pi
R}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}$
$\displaystyle=\frac{\text{1}}{\eta}\frac{\sigma}{\text{2}\sqrt{\sigma-\text{1}}}\Bigg{\\{}\frac{\operatorname{i}}{\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}}}\Bigg{[}\text{1}-e^{-\frac{\text{4}(\sigma-\text{2})}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}e^{\frac{\operatorname{i}\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\Bigg{]}$
(36)
$\displaystyle-\frac{\operatorname{i}}{\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}}}\Bigg{[}\text{1}-e^{-\frac{\text{4}(\sigma-\text{2})}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}e^{\frac{\operatorname{i}\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\Bigg{]}\Bigg{\\}}\text{,}\quad\sigma>\text{1}$

The impulse response function, $h(t)$, offered by Eq. (36) for
$\sigma>\text{1}$ is a real-valued function given that the imaginary parts of
the right-hand side of Eq. (36) cancel.

For the case when $\sigma<\text{1}$
$\left(\frac{\rho_{p}}{\rho_{f}}<\frac{\text{5}}{\text{8}}\right)$, the
parameters $a$ and $b$ appearing in Eq. (35) are real-valued and are offered
by Eq. (28). In this case ($\sigma<\text{1}$), Eq. (35) becomes

$\displaystyle h(t)=\frac{\text{3}\pi
R}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}$
$\displaystyle=\frac{\text{1}}{\eta}\frac{\sigma}{\text{2}\sqrt{\text{1}-\sigma}}\Bigg{\\{}-\frac{\text{1}}{\text{1}+\sqrt{\text{1}-\sigma}}\Bigg{[}\text{1}-e^{\frac{\text{4}}{\sigma^{\text{2}}}\left(\text{2}-\sigma+\text{2}\sqrt{\text{1}-\sigma}\right)\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}+\sqrt{\text{1}-\sigma}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\Bigg{]}$
(37)
$\displaystyle+\frac{\text{1}}{\text{1}-\sqrt{\text{1}-\sigma}}\Bigg{[}\text{1}-e^{\frac{\text{4}}{\sigma^{\text{2}}}\left(\text{2}-\sigma-\text{2}\sqrt{\text{1}-\sigma}\right)\frac{t}{\tau_{f}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}-\sqrt{\text{1}-\sigma}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\Bigg{]}\Bigg{\\}}\text{,}\quad\sigma<\text{1}$

Figure 3 plots the normalized time-derivative of the mean-square displacement
$\frac{\text{3}\pi R\eta}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}=\eta h(t)$ of Brownian particles
with mass $m$ and radius $R$ immersed in a Newtonian, viscous fluid with
viscosity $\eta$ for various values of the parameter $\sigma$ which accounts
for hydrodynamic memory. The curves for $\sigma=$ 0.9 and 0.99 are computed
with Eq. (37); whereas the curves for $\sigma=$ 1.01, 1.5, 2 and 2.5 are
computed with Eq. (36) and follow qualitatively the trend of the
experimentally measured time-derivative of the mean-square displacement
reported by [12] and [38].

Figure 3: Normalized time-derivative of the mean-square displacement of
Brownian particles immersed in a Newtonian, viscous fluid with viscosity
$\eta$; when effects from hydrodynamic memory are accounted. The curves for
$\sigma=\frac{\text{8}}{\text{9}}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)=$
0.9 and 0.99 are computed with Eq. (37); whereas, the curves for $\sigma=$
1.01, 1.5, 2 and 2.5 are computed with Eq. (36) showing the continuity of the
solution as parameter $\sigma$ crosses the value of $\sigma=\text{1}$.

In terms of the two-parameter Mittag–Leffler function given by the identity
(34), an alternative expression to Eq. (35) for the time-derivative of the
mean-square displacement of Brownian particles in a Newtonian fluid with
hydrodynamic memory is:

$\displaystyle\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}$
$\displaystyle=\frac{NK_{B}T}{\text{3}\pi
R}h(t)=\frac{\text{2}NK_{B}T}{M}\frac{\sqrt{t}}{b-a}$ (38)
$\displaystyle\times\left[E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{3}}}{{\text{2}}}}\left(-a\sqrt{t}\right)-E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{3}}}{{\text{2}}}}\left(-b\sqrt{t}\right)\right]$

where again $a$ and $b$ are given by Eq. (10).

##  IV Mean-Square Displacement of Brownian Particles in a Newtonian Viscous
Fluid with Hydrodynamic Memory

According to the properties of the Laplace transform of the derivative of a
function

$\mathcal{L}\left\\{\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}\right\\}=s\left\langle\Delta
r^{\text{2}}(s)\right\rangle-\left\langle\Delta
r^{\text{2}}(\text{0})\right\rangle$ (39)

Given that at the time origin $(t=\text{0})$, $\left\langle\Delta
r^{\text{2}}(\text{0})\right\rangle=\text{0}$, substitution of Eq. (39) into
Eq. (30) gives

$\left\langle\Delta r^{\text{2}}(s)\right\rangle=\frac{NK_{B}T}{\text{3}\pi
R}\frac{\mathcal{J}(s)}{s}=\frac{NK_{B}T}{\text{3}\pi
R}\mathcal{L}\left\\{J(t)\right\\}$ (40)

where $\mathcal{J}(s)$ is the complex dynamic compliance of the rheological
network shown in Fig. 1 and expressed by Eq. (14); whereas, $J(t)$ is the
creep compliance of the rheological network in the time-domain. Accordingly
from Eq. (40) in association with Eq. (14)

$\displaystyle\mathcal{L}\left\\{J(t)\right\\}$
$\displaystyle=\frac{\text{1}}{m_{R}}\frac{\text{1}}{s^{\text{2}}}\frac{\text{1}}{\left(s+\frac{\mu_{\text{3}/\text{2}}}{m_{R}}s^{\text{1}/\text{2}}+\frac{\eta}{m_{R}}\right)}$
(41) $\displaystyle=\frac{\text{6}\pi
R}{M}\frac{\text{1}}{s^{\text{2}}}\frac{\text{1}}{\left(s+\frac{z}{M}s^{\text{1}/\text{2}}+\frac{\zeta}{M}\right)}$

The inverse Laplace transform of Eq. (41) is evaluated by employing the
fraction expansion

$\displaystyle\frac{\text{1}}{s^{\text{2}}\left(s+(a+b)s^{\text{1}/\text{2}}+ab\right)}$
(42)
$\displaystyle=\frac{\text{1}}{b-a}\left[\frac{\text{1}}{s^{\text{2}}\left(s^{\text{1}/\text{2}}+a\right)}-\frac{\text{1}}{s^{\text{2}}\left(s^{\text{1}/\text{2}}+b\right)}\right]$

together with the known result [34]

$\mathcal{L}^{-\text{1}}\left\\{\frac{s^{p-q}}{s^{p}+a}\right\\}=t^{q-\text{1}}E_{p,q}(-at^{p})\text{,}\quad
p<q\in\mathbb{R}.$ (43)

Upon taking $p=\frac{\text{1}}{\text{2}}$ and $q=\frac{\text{5}}{\text{2}}$,
the Laplace transform of any of the two fractions appearing on the right-hand
side of Eq. (42) is

$\displaystyle\mathcal{L}^{-\text{1}}\left\\{\frac{s^{\text{1}/\text{2}-\text{5}/\text{2}}}{s^{\text{1}/\text{2}}+a}\right\\}$
$\displaystyle=\mathcal{L}^{-\text{1}}\left\\{\frac{\text{1}}{s^{\text{2}}\left(s^{\text{1}/\text{2}}+a\right)}\right\\}$
(44)
$\displaystyle=t^{\text{3}/\text{2}}E_{\text{1}/\text{2}\text{,}\,\text{5}/\text{2}}(-at^{\text{1}/\text{2}})$

With reference to the denominators of the left-hand side of Eq. (42) and the
right-hand side of Eq. (41), the solution of the system of equations
$a+b=\frac{z}{M}$ and $ab=\frac{\zeta}{M}$ yields the values of $a$ and $b$
offered by Eq. (10). Accordingly, by virtue of Eqs. (42) and (44), the mean-
square displacement of Brownian particles in a Newtonian viscous fluid with
hydrodynamic memory is

$\displaystyle\left\langle\Delta r^{\text{2}}(t)\right\rangle$
$\displaystyle=\frac{NK_{B}T}{\text{3}\pi
R}J(t)=\frac{\text{2}NK_{B}T}{M}\frac{t^{\text{3}/\text{2}}}{b-a}$ (45)
$\displaystyle\times\left[E_{\text{1}/\text{2}\text{,}\,\text{5}/\text{2}}(-a\sqrt{t})-E_{\text{1}/\text{2}\text{,}\,\text{5}/\text{2}}(-b\sqrt{t})\right]$

where the parameters $a$ and $b$ are offered by Eq. (10). Equation (45) is a
compact expression of the mean-square displacement in terms of the two-
parameter Mittag–Leffler function and resembles the structure of its time-
derivative given by Eq. (38) and that of the velocity autocorrelation function
given by Eq. (21). As we go up to higher-order time-derivatives of Eq. (45),
the power law of the time in front of the brackets reduces by one (from
$t^{\nicefrac{{\text{3}}}{{\text{2}}}}$ to
$t^{\nicefrac{{\text{1}}}{{\text{2}}}}$ to
$t^{-\nicefrac{{\text{1}}}{{\text{2}}}}$) together with the second index of
the Mittag–Leffler function (from $\nicefrac{{\text{5}}}{{\text{2}}}$ to
$\nicefrac{{\text{3}}}{{\text{2}}}$ to $\nicefrac{{\text{1}}}{{\text{2}}}$).

By employing the recurrence relation of the Mittag–Leffler function [33, 34]:

$E_{p,q}(x)=\frac{\text{1}}{x}E_{p,q-p}(x)-\frac{\text{1}}{x\Gamma(q-p)}$ (46)
$E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{5}}}{{\text{2}}}}\left(-at^{\nicefrac{{\text{1}}}{{\text{2}}}}\right)=-\frac{\text{1}}{a\sqrt{t}}\left[E_{\nicefrac{{\text{1}}}{{\text{2}}},\text{2}}\left(-a\sqrt{t}\right)-\text{1}\right]$
(47)

where $E_{\nicefrac{{\text{1}}}{{\text{2}}},\text{2}}\left(-a\sqrt{t}\right)$
in (47) is expressed as:

$E_{\nicefrac{{\text{1}}}{{\text{2}}},\text{2}}\left(-a\sqrt{t}\right)=-\frac{\text{1}}{a\sqrt{t}}\left[E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{3}}}{{\text{2}}}}\left(-a\sqrt{t}\right)-\frac{\text{1}}{\Gamma\left(\nicefrac{{\text{3}}}{{\text{2}}}\right)}\right]$
(48)

Substitution of Eq. (48) into (47) after using that
$\Gamma\left(\nicefrac{{\text{3}}}{{\text{2}}}\right)=\nicefrac{{\sqrt{\pi}}}{{\text{2}}}$,
an alternative expression of the mean-square displacement given by Eq. (45)
is:

$\left\langle\Delta r^{\text{2}}(t)\right\rangle=\frac{NK_{B}T}{\text{3}\pi
R}J(t)=\frac{\text{2}NK_{B}T}{M}\left\\{\frac{t}{ab}+\frac{\sqrt{t}}{b-a}\frac{\text{2}}{\sqrt{\pi}}\left(\frac{\text{1}}{b^{\text{2}}}-\frac{\text{1}}{a^{\text{2}}}\right)+\frac{\sqrt{t}}{b-a}\left[\frac{\text{1}}{a^{\text{2}}}E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{3}}}{{\text{2}}}}\left(-a\sqrt{t}\right)-\frac{\text{1}}{b^{\text{2}}}E_{\nicefrac{{\text{1}}}{{\text{2}}},\nicefrac{{\text{3}}}{{\text{2}}}}\left(-b\sqrt{t}\right)\right]\right\\}$
(49)

By employing once again the recurrence relation given by Eq. (46) in
association with the identity:

$E_{\nicefrac{{\text{1}}}{{\text{2}}},\text{1}}(x)=e^{x^{\text{2}}}\text{Erfc}(-x),$
(50)

the mean-square displacement of Brownian particles immersed in a Newtonian
viscous fluid with hydrodynamic memory can be expressed as:

$\left\langle\Delta r^{\text{2}}(t)\right\rangle=\frac{NK_{B}T}{\text{3}\pi
R}J(t)=\frac{NK_{B}T}{\text{3}\pi
R}\frac{\text{1}}{\eta}\left\\{t-\frac{\text{2}}{\sqrt{\pi}}\frac{z}{\zeta}\sqrt{t}+\frac{z^{\text{2}}-\zeta
M}{\zeta^{\text{2}}}+\frac{\zeta}{M}\frac{\text{1}}{b-a}\left[\frac{\text{1}}{b^{\text{3}}}e^{b^{\text{2}}t}\text{Erfc}\left(b\sqrt{t}\right)-\frac{\text{1}}{a^{\text{3}}}e^{a^{\text{2}}t}\text{Erfc}\left(a\sqrt{t}\right)\right]\right\\}$
(51)

which is the classical expression known in the literature [38, 13]. In
deriving Eq. (51) we used the identities:

$\frac{a+b}{ab}=\frac{z}{\zeta}=\sqrt{\tau_{f}}\enskip\text{and}\enskip\frac{a^{\text{2}}+ab+b^{\text{2}}}{a^{\text{2}}b^{\text{2}}}=\frac{z^{\text{2}}-\zeta
M}{\zeta^{\text{2}}}$ (52)

which derive from the definitions of $a$ and $b$ given by Eq. (10).

When $\frac{\rho_{p}}{\rho_{f}}>\frac{\text{5}}{\text{8}}$, parameter
$\sigma>\text{1}$ and the parameters $a$ and $b$ appearing in Eq. (51) are
complex valued and are offered by Eq. (26). In this case ($\sigma>\text{1}$)
Eq. (51) becomes:

$\displaystyle\left\langle\Delta r^{\text{2}}(t)\right\rangle$
$\displaystyle=\frac{NK_{B}T}{\text{3}\pi R}J(t)=\frac{NK_{B}T}{\text{3}\pi
R}\frac{\tau_{f}}{\eta}\left\\{\frac{t}{\tau_{f}}-\frac{\text{2}}{\sqrt{\pi}}\sqrt{\frac{t}{\tau_{f}}}+\left(\text{1}-\frac{\sigma}{\text{4}}\right)+\frac{\operatorname{i}}{\text{8}}\frac{\sigma^{\text{3}}}{\sqrt{\sigma-\text{1}}}e^{-\frac{\text{4}(\sigma-\text{2})}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\right.$
(53)
$\displaystyle\times\left[\frac{e^{-\frac{\operatorname{i}\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}}{\left(\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}}\right)^{\text{3}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}-\operatorname{i}\sqrt{\sigma-\text{1}}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\left.-\frac{e^{\frac{\operatorname{i}\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}}{\left(\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}}\right)^{\text{3}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}+\operatorname{i}\sqrt{\sigma-\text{1}}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\right]\right\\},\enskip\sigma>\text{1}$

The mean-square displacement offered by Eq. (53) for $\sigma>$ 1 is a real-
valued function given that the imaginary parts of the right-hand side of Eq.
(53) cancel. For the case when $\sigma<\text{1}$
$\left(\frac{\rho_{p}}{\rho_{f}}<\frac{\text{5}}{\text{8}}\right)$, the
parameters $a$ and $b$ appearing in Eq. (51) are real-valued and are offered
by Eq. (28). In this case ($\sigma<$ 1), Eq. (51) becomes

$\displaystyle\left\langle\Delta r^{\text{2}}(t)\right\rangle$
$\displaystyle=\frac{NK_{B}T}{\text{3}\pi R}J(t)=\frac{NK_{B}T}{\text{3}\pi
R}\frac{\tau_{f}}{\eta}\left\\{\frac{t}{\tau_{f}}-\frac{\text{2}}{\sqrt{\pi}}\sqrt{\frac{t}{\tau_{f}}}+\left(\text{1}-\frac{\sigma}{\text{4}}\right)-\frac{\sigma^{\text{3}}}{\text{8}\sqrt{\text{1}-\sigma}}e^{\frac{\text{4}(\sigma-\text{2})}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}\right.$
(54)
$\displaystyle\times\left[\frac{e^{-\frac{\text{8}\sqrt{\text{1}-\sigma}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}}{\left(\text{1}-\sqrt{\text{1}-\sigma}\right)^{\text{3}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}-\sqrt{\text{1}-\sigma}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\left.-\frac{e^{\frac{\text{8}\sqrt{\sigma-\text{1}}}{\sigma^{\text{2}}}\frac{t}{\tau_{f}}}}{\left(\text{1}+\sqrt{\text{1}-\sigma}\right)^{\text{3}}}\text{Erfc}\left(\frac{\text{2}\left(\text{1}+\sqrt{\text{1}-\sigma}\right)}{\sigma}\sqrt{\frac{t}{\tau_{f}}}\right)\right]\right\\},\quad\sigma<\text{1}$

Figure 4 plots the normalized mean-square displacement
$\frac{\text{3}\pi\eta^{\text{2}}}{NK_{B}TR\rho_{f}}\left\langle\Delta
r^{\text{2}}(t)\right\rangle=\frac{\eta^{\text{2}}}{R^{\text{2}}\rho_{f}}J(t)$
of Brownian particles with mass $m$ and radius $R$ immersed in a Newtonian,
viscous fluid with viscosity $\eta$ for various values of the parameter
$\sigma$ which accounts for hydrodynamic memory.

Figure 4: Normalized mean-square displacement of Brownian particles immersed
in a Newtonian, viscous fluid with viscosity $\eta$; when effects of
hydrodynamic memory are accounted. The curves for
$\sigma=\frac{\text{8}}{\text{9}}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)=$
0.9 and 0.99 are computed with Eq. (54); whereas, the curves for $\sigma=$
1.01, 1.5, 2 and 2.5 are computed with Eq. (53) showing the continuity of the
solution as parameter $\sigma$ crosses the value of $\sigma=\text{1}$.

##  V Brownian Motion in Viscoelastic Materials with Hydrodynamic Memory

The motion of a solid particle within a viscoelastic material is described by
nonlinear continuum mechanics equations [54, 55] together with rheological
models expressed in a corotating coordinate frame for constructing objective
rheological equations of state. Nevertheless, when the particle motion happens
over small distances in association with low Reynolds numbers, the
nonlinearities that manifest in the dynamic equilibrium and constitutive
equations can be neglected [3, 55, 9, 13, 56, 57]. Accordingly, the unsteady
motion of Brownian particles suspended in a linear, isotropic viscoelastic
material is described by a generalized Langevin equation

$m\frac{\mathrm{d}v(t)}{\mathrm{d}t}=-\int_{\text{0}^{-}}^{t}\zeta(t-\xi)v(\xi)\,\mathrm{d}\xi+f_{R}(t)$
(55)

where $m$ is the particle mass, $v(t)=\frac{\mathrm{d}r(t)}{\mathrm{d}t}$ is
the particle velocity and $f_{R}(t)$ are the random forces acting on the
Brownian particle. The convolution integral in Eq. (55) represents the drag
force on the particle as it moves through the viscoelastic fluid over small
distances and $\zeta(t-\xi)$ is the relaxation kernel of the Brownian
particle–viscoelastic material system.

In their seminal paper Zwanzig and Bixon [9] built upon Stoke’s and
Boussinesq’s solution for the unsteady motion of a sphere within a Newtonian
viscous fluid and upon decomposing the particle velocity into its Fourier
components, $v(\omega)=\int_{-\infty}^{\infty}v(t)e^{-i\omega
t}\,\mathrm{d}t$, reached a solution for the Fourier transform of the
relaxation kernel
$\mathcal{Z}(\omega)=\int_{-\infty}^{\infty}\zeta(t)e^{-i\omega
t}\,\mathrm{d}t$. The mathematical expression for the impedance
$\mathcal{Z}(\omega)$ derived by [9] involves longitudinal and transverse
Hankel functions after applying boundary conditions to the solution of a
vector Helmholtz equation expressed in spherical harmonics. While the solution
is valid for the general linear viscoelastic material, its complexity does not
allow for closed form results other than at selective limited cases [9].

Clercx and Schram [13] examined the motion of a Brownian particle suspended in
a Newtonian viscous fluid under the influence of a harmonic potential that
exerts a linear restoring force $kr(t)$, where $r(t)$ is the particle
displacement and $k$ is the constant of proportionality (stiffness). With this
configuration the drag force on the Brownian particle including the effects of
backflow are associated solely with the surrounding viscous fluid appearing in
the right-hand side of the generalized Langevin Eq. (4), or of its alternative
form given by Eq. (8).

With reference to Eq. (8) and after expressing the particle velocity
$v(t)=\frac{\mathrm{d}r(t)}{\mathrm{d}t}$, the generalized Langevin equation
for the Brownian particle in a harmonic trap (Kelvin-Voigt solid) when effects
from hydrodynamic memory are accounted simplifies to

$M\frac{\mathrm{d}^{\text{2}}r(t)}{\mathrm{d}t^{\text{2}}}+\text{6}\pi
R^{\text{2}}\sqrt{\rho_{f}\eta}\frac{\mathrm{d}^{\nicefrac{{\text{3}}}{{\text{2}}}}r(t)}{\mathrm{d}t^{\nicefrac{{\text{3}}}{{\text{2}}}}}+\text{6}\pi
R\eta\frac{\mathrm{d}r(t)}{\mathrm{d}t}+kr(t)=f_{R}(t)$ (56)

where again $M=\frac{\text{4}}{\text{3}}\pi
R^{\text{3}}(\rho_{p}+\frac{\text{1}}{\text{2}}\rho_{f})=\frac{\text{3}}{\text{2}}\pi
R^{\text{3}}\rho_{f}\sigma$ with $\sigma$ offered by Eq. (24). Clercx and
Schram [13], after deriving the appropriate correlation function $\left\langle
f_{R}(t_{\text{1}})f_{R}(t_{\text{2}})\right\rangle$ of the random forces,
derived closed form expressions of the velocity autocorrelation function and
of the mean-square displacements by computing ensemble averages of the random
Brownian process described by Eq. (56).

The concept of a rheological analogue advanced in this paper can be also
applied for the case of Brownian motion in a harmonic trap (Kelvin-Voigt
solid) when hydrodynamic memory is accounted. The synthesis of the proposed
macroscopic mechanical network is suggested from the left-hand side of the
Langevin equation (56), which consists of a inertia term; a
$\nicefrac{{\text{3}}}{{\text{2}}}$-order fractional derivative of the
displacement term, a viscous term and an elastic term acting in parallel to
balance the random force $f_{R}(t)$. Accordingly, by applying the
viscous–viscoelastic correspondence principle for Brownian motion [25, 26] the
mean-square displacement and its time derivatives of Brownian particles in a
harmonic trap when hydrodynamic memory is accounted are evaluated by merely
computing the creep compliance, $J(t)=\frac{\text{3}\pi
R}{NK_{B}T}\left\langle\Delta r^{\text{2}}(t)\right\rangle$, impulse response
function, $h(t)=\frac{\text{3}\pi
R}{NK_{B}T}\frac{\mathrm{d}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t}$ and impulse strain rate response
function $\psi(t)=\frac{\text{3}\pi
R}{NK_{B}T}\frac{\mathrm{d}^{\text{2}}\left\langle\Delta
r^{\text{2}}(t)\right\rangle}{\mathrm{d}t^{\text{2}}}=\frac{\text{6}\pi
R}{NK_{B}T}\left\langle v(\text{0})v(t)\right\rangle$ of the mechanical
network shown in Fig. 5. Given the parallel connection of the
spring–dashpot–interpot–inerter network shown in Fig. 5, its constitutive law
is

$\tau(t)=G\gamma(t)+\eta\frac{\mathrm{d}\gamma(t)}{\mathrm{d}t}+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}\frac{\mathrm{d}^{\nicefrac{{\text{3}}}{{\text{2}}}}\gamma(t)}{\mathrm{d}t^{\nicefrac{{\text{3}}}{{\text{2}}}}}+m_{R}\frac{\mathrm{d}^{\text{2}}\gamma(t)}{\mathrm{d}t^{\text{2}}}$
(57)

Figure 5: The rheological analogue for Brownian motion in a viscously damped
harmonic trap with hydrodynamic memory. A spring with shear modulus $G$, a
dashpot with viscosity $\eta$, a Scott-Blair element with material constant
$\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}=R\sqrt{\rho_{f}\eta}$ and an inerter
with distributed inertance $m_{R}=\frac{M}{\text{6}\pi
R}=\frac{\text{2}}{\text{9}}R^{\text{2}}\rho_{f}\left(\frac{\rho_{p}}{\rho_{f}}+\frac{\text{1}}{\text{2}}\right)=\frac{\text{1}}{\text{4}}\sigma
R^{\text{2}}\rho_{f}$ are connected in parallel; the
spring–dashpot–inertpot–inerter parallel connection.

The Laplace transform of Eq. (57) gives

$\tau(s)=\mathcal{G}(s)\gamma(s)=(G+\eta
s+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+m_{R}s^{\text{2}})\gamma(s)$
(58)

where $\mathcal{G}(s)=G+\eta
s+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+m_{R}s^{\text{2}}$
is the complex dynamic modulus. The complex dynamic compliance of the
rheological network shown in Fig. 5 is

$\displaystyle\mathcal{J}(s)$
$\displaystyle=\frac{\text{1}}{\mathcal{G}(s)}=\frac{\text{1}}{(G+\eta
s+\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+m_{R}s^{\text{2}})}$
(59)
$\displaystyle=\frac{\text{1}}{m_{R}}\frac{\text{1}}{s^{\text{2}}+\frac{\mu_{\nicefrac{{\text{3}}}{{\text{2}}}}}{m_{R}}s^{\nicefrac{{\text{3}}}{{\text{2}}}}+\frac{\eta}{m_{R}}s+\omega^{\text{2}}_{R}}$

where $\omega_{R}=\sqrt{\frac{G}{m_{R}}}=\sqrt{\frac{k}{M}}$ is the undamped
natural frequency of the trapped Brownian particle. The Laplace transform of
the velocity autocorrelation function of Brownian particles in a harmonic trap
when hydrodynamic memory is accounted is offered by Eq. (15) where now the
complex dynamic compliance $\mathcal{J}(s)$ is offered by Eq. (59). Inverse
Laplace transform of Eq. (15) yields the evolved expression for the velocity
autocorrelation function first presented by [13].

##  VI Summary

In this paper we first recognized that for Brownian particles immersed in a
Newtonian viscous fluid, the hydrodynamic memory term associated with the
backflow is essentially a $\nicefrac{{\text{1}}}{{\text{2}}}$-fractional
derivative of the velocity of the Brownian particle. Accordingly, the
hydrodynamic memory term can be expressed with a macroscopic Scott-Blair
fractional element. Subsequently, we built on a recently published
viscous–viscoelastic correspondence principle for Brownian motion [25, 26] and
presented macroscopic rheological networks for Brownian motion of particles
either immersed in a Newtonian viscous fluid or trapped in a harmonic
potential. The synthesis of the rheological network is suggested from the
structure of the generalized Langevin equation and simplifies appreciably the
calculations for evaluating the mean-square displacement and its time
derivatives.

## References

  * Einstein [1905] A. Einstein, Annalen der Physik 17, 549 (1905).
  * Langevin [1908] P. Langevin, Compt. Rendus 146, 530 (1908).
  * Landau and Lifshitz [1959] L. D. Landau and E. M. Lifshitz,  _Course of Theoretical Physics Vol. 6 Fluid Mechanies_ (Pergamon Press, 1959).
  * Attard [2012] P. Attard,  _Non-equilibrium thermodynamics and statistical mechanics: Foundations and applications_ (OUP Oxford, 2012).
  * Kalmykov and Coffey [2017] Y. P. Kalmykov and W. T. Coffey,  _The Langevin Equation_ (World Scientific Publishing Company, 2017).
  * Uhlenbeck and Ornstein [1930] G. E. Uhlenbeck and L. S. Ornstein, Physical Review 36, 823 (1930).
  * Wang and Uhlenbeck [1945] M. C. Wang and G. E. Uhlenbeck, Reviews of Modern Physics 17, 323 (1945).
  * Lighthill [1958] M. J. Lighthill,  _An introduction to Fourier analysis and generalised functions_ (Cambridge University Press, 1958).
  * Zwanzig and Bixon [1970] R. Zwanzig and M. Bixon, Physical Review A 2, 2005 (1970).
  * Widom [1971] A. Widom, Physical Review A 3, 1394 (1971).
  * Hinch [1975] E. J. Hinch, Journal of Fluid Mechanics 72, 499 (1975).
  * Paul and Pusey [1981] G. L. Paul and P. N. Pusey, Journal of Physics A: Mathematical and General 14, 3301 (1981).
  * Clercx and Schram [1992] H. J. H. Clercx and P. P. J. M. Schram, Physical Review A 46, 1942 (1992).
  * Oldham and Spanier [1974] K. Oldham and J. Spanier,  _The Fractional Calculus. Mathematics in science and engineering_ , Vol. III (Academic Press Inc., San Diego, CA, 1974).
  * Samko  _et al._ [1974] S. G. Samko, A. A. Kilbas, and O. I. Marichev,  _Fractional Integrals and Derivatives; Theory and Applications_ , Vol. 1 (Gordon and Breach Science Publishers, Amsterdam, 1974).
  * Miller and Ross [1993] K. S. Miller and B. Ross,  _An introduction to the fractional calculus and fractional differential equations_ (Wiley, New York, NY, 1993).
  * Podlubny [1998] I. Podlubny,  _Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications_ (Elsevier, 1998).
  * Mainardi [2010] F. Mainardi,  _Fractional calculus and waves in linear viscoelasticity: An introduction to mathematical models_ (Imperial College Press - World Scientific, London, UK, 2010).
  * Makris [2021a] N. Makris, Fractal and Fractional 5, 18 (2021a).
  * Nutting [1921] P. G. Nutting, Proceedings American Society for Testing Materials 21, 1162 (1921).
  * Gemant [1936] A. Gemant, Physics 7, 311 (1936).
  * Scott Blair [1944] G. W. Scott Blair,  _A survey of general and applied rheology_ (Isaac Pitman & Sons, 1944).
  * Scott Blair [1947] G. W. Scott Blair, Journal of Colloid Science 2, 21 (1947).
  * Li and Raizen [2013] T. Li and M. G. Raizen, Annalen der Physik 525, 281 (2013).
  * Makris [2020] N. Makris, Physical Review E 101, 052139 (2020).
  * Makris [2021b] N. Makris, Soft Matter 17, 5410 (2021b).
  * Gemant [1938] A. Gemant, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 25, 540 (1938).
  * Smith [2002] M. C. Smith, IEEE Transactions on Automatic Control 47, 1648 (2002).
  * Papageorgiou and Smith [2005] C. Papageorgiou and M. C. Smith, in  _Proceedings of the 44th IEEE Conference on Decision and Control_ (IEEE, 2005) pp. 3351–3356.
  * Makris [2017] N. Makris, Journal of Engineering Mechanics 143, 04017123 (2017).
  * Makris [2018] N. Makris, Meccanica 53, 2237 (2018).
  * Erdélyi [1954] A. Erdélyi, ed.,  _Bateman Manuscript Project, Tables of Integral Transforms Vol I_ (McGraw-Hill, New York, NY, 1954).
  * Erdélyi [1953] A. Erdélyi, ed.,  _Bateman Manuscript Project, Higher Transcendental Functions Vol III_ (McGraw-Hill, New York, NY, 1953).
  * Gorenflo  _et al._ [2014] R. Gorenflo, A. A. Kilbas, F. Mainardi, S. V. Rogosin, _et al._ ,  _Mittag-Leffler functions, related topics and applications Vol. 2_ (Springer, 2014).
  * Franosch  _et al._ [2011] T. Franosch, M. Grimm, M. Belushkin, F. M. Mor, G. Foffi, L. Forró, and S. Jeney, Nature 478, 85 (2011).
  * Domínguez-García  _et al._ [2014] P. Domínguez-García, F. Cardinaux, E. Bertseva, L. Forró, F. Scheffold, and S. Jeney, Physical Review E 90, 060301 (2014).
  * Jannasch  _et al._ [2011] A. Jannasch, M. Mahamdeh, and E. Schäffer, Physical Review Letters 107, 228301 (2011).
  * Weitz  _et al._ [1989] D. A. Weitz, D. J. Pine, P. N. Pusey, and R. Tough, Physical Review Letters 63, 1747 (1989).
  * Segre and Pusey [1996] P. Segre and P. Pusey, Physical Review Letters 77, 771 (1996).
  * Sperl [2005] M. Sperl, Physical Review E 71, 060401 (2005).
  * Safdari  _et al._ [2017] H. Safdari, A. G. Cherstvy, A. V. Chechkin, A. Bodrova, and R. Metzler, Physical Review E 95, 012120 (2017).
  * Khan and Mason [2014] M. Khan and T. G. Mason, Physical Review E 89, 042309 (2014).
  * Ghosh and Krishnamurthy [2018] K. Ghosh and C. V. Krishnamurthy, Physical Review E 98, 052115 (2018).
  * Kenkre  _et al._ [1981] V. Kenkre, R. Kühne, and P. Reineker, Zeitschrift für Physik B Condensed Matter 41, 177 (1981).
  * Bian  _et al._ [2016] X. Bian, C. Kim, and G. E. Karniadakis, Soft Matter 12, 6331 (2016).
  * Giesekus [1995] H. Giesekus, Rheologica Acta 34, 2 (1995).
  * Makris and Kampas [2009] N. Makris and G. Kampas, Rheologica Acta 48, 815 (2009).
  * Makris and Efthymiou [2020] N. Makris and E. Efthymiou, Rheologica Acta 59, 849 (2020).
  * Clough and Penzien [1970] R. W. Clough and J. Penzien,  _Dynamics of structures_ (McGraw-Hill, New York, NY, 1970).
  * Harris and Crede [1976] C. M. Harris and C. E. Crede,  _Shock and vibration handbook_ , 2nd ed. (McGraw-Hill, New York, NY, 1976).
  * Oppenheim and Schafer [1975] A. V. Oppenheim and R. W. Schafer,  _Digital signal processing_ (Prentice-Hall, Inc., Englewood Cliffs, NJ, 1975).
  * Reid [1983] G. J. Reid, McGraw-Hill Series in Electrical Engineering (1983).
  * Borovička and Hansen [2016] J. Borovička and L. P. Hansen, in  _Handbook of Macroeconomics_ , Vol. 2 (Elsevier, 2016) pp. 1641–1696.
  * Gurtin [1982] M. E. Gurtin,  _An introduction to continuum mechanics_ (Academic press, 1982).
  * Bird  _et al._ [1987] R. B. Bird, R. C. Armstrong, and O. Hassager,  _Dynamics of polymeric liquids. Vol. 1: Fluid mechanics_ , 2nd ed. (Wiley, New York, NY, 1987).
  * Liverpool and MacKintosh [2005] T. Liverpool and F. MacKintosh, Physical review letters 95, 208303 (2005).
  * Atakhorrami  _et al._ [2008] M. Atakhorrami, D. Mizuno, G. Koenderink, T. Liverpool, F. MacKintosh, and C. Schmidt, Physical Review E 77, 061508 (2008).