published in Physical Review Fluids 3, 103301 (2018)

# Time-dependent active microrheology in dilute colloidal suspensions

Sebastian Leitmann  Institut für Theoretische Physik, Universität Innsbruck,
Technikerstraße 21A, A-6020 Innsbruck, Austria  Suvendu Mandal  Institut für
Theoretische Physik, Universität Innsbruck, Technikerstraße 21A, A-6020
Innsbruck, Austria  Matthias Fuchs  Fachbereich Physik, Universität Konstanz,
78457 Konstanz, Germany  Antonio M. Puertas  Group of Complex Fluids Physics,
Department of Applied Physics, University of Almería, 04120 Almería, Spain
Thomas Franosch  Institut für Theoretische Physik, Universität Innsbruck,
Technikerstraße 21A, A-6020 Innsbruck, Austria<EMAIL_ADDRESS>
###### Abstract

In a microrheological set-up a single probe particle immersed in a complex
fluid is exposed to a strong external force driving the system out of
equilibrium. Here, we elaborate analytically the time-dependent response of a
probe particle in a dilute suspension of Brownian particles to a large step-
force, exact in first order of the density of the bath particles. The time-
dependent drift velocity approaches its stationary state value exponentially
fast for arbitrarily small driving in striking contrast to the power-law
prediction of linear response encoded in the long-time tails of the velocity
autocorrelation function. We show that the stationary-state behavior depends
nonanalytically on the driving force and connect this behavior to the
persistent correlations in the equilibrium state. We argue that this relation
holds generically. Furthermore, we elaborate that the fluctuations in the
direction of the force display transient superdiffusive behavior.

##  I Introduction

While the static and dynamic properties of interacting many-particle systems
in equilibrium encode the linear response via the fluctuation-dissipation
theorem, corresponding principles generally applicable for systems driven far
from equilibrium remain a grand challenge in statistical physics. Soft matter
systems are ideally suited to study such nonequilibrium phenomena, since their
defining characteristic is that they are strongly susceptible to forces.
Conceptually, the simplest experimentally realizable system then consists of
an interacting colloidal suspension that is driven out of equilibrium by a
strong external force acting on a single probe particle. This setup
constitutes the basic paradigm for active microrheology with the principle
goal to infer material properties beyond the linear regime Squires and Brady
(2005); Khair and Brady (2006); Wilson  _et al._ (2009); Zia and Brady (2010);
Sriram  _et al._ (2010); Wilson and Poon (2011); Zia and Brady (2012); Swan
and Zia (2013); Puertas and Voigtmann (2014); Zia and Brady (2015); Furst and
Squires (2017).

The nonlinear mobility in the stationary state $\mu(\text{Pe})$ for a
suspension of interacting Brownian particles has been derived to first order
of the density of bath particles $n$ in the seminal work by Squires and Brady
Squires and Brady (2005) in terms of an asymptotic expansion

$\displaystyle\mu(\text{Pe})/\mu=1-\frac{2\pi
n\sigma^{3}}{3}\frac{D_{a}}{D_{r}}\biggl{[}1-\frac{2}{15}\text{Pe}^{2}+\frac{1}{8}|\text{Pe}|^{3}-\frac{128}{1575}\text{Pe}^{4}+\mathcal{O}(|\text{Pe}|^{5})\biggr{]},\quad\text{Pe}\to
0.$ (1)

where $\mu$ is the mobility of the probe particle at infinite dilution with
diffusion coefficient $D_{a}$, and the Péclet number Pe is a suitable
dimensionless measure for the driving force (see below). The motion of the
bath particles with diffusion coefficient $D_{b}$ enters in terms of the
diffusion coefficient of the relative motion, $D_{r}=D_{a}+D_{b}$.

The leading term of the corresponding asymptotic expansion of the fluctuations
around the average drift has also been achieved Zia and Brady (2010),
revealing a long-time diffusion coefficient that can become arbitrarily large
for strong driving. Active microrheology has found fruitful applications in
dense colloidal systems in the vicinity of the glass transition, in
particular, computer simulations have revealed (transient) superdiffusion and
enhanced diffusivities Winter  _et al._ (2012); Winter and Horbach (2013).
Certain phenomena in the nonlinear regime have also been rationalized within a
mode-coupling approach Gazuz  _et al._ (2009); Gnann  _et al._ (2011); Gnann
and Voigtmann (2012); Harrer  _et al._ (2012); Gazuz and Fuchs (2013); Wang
_et al._ (2014); Gruber  _et al._ (2016), continuous-time random walks Jack
_et al._ (2008); Schroer and Heuer (2013a, b); Burioni  _et al._ (2014),
Langevin equations Démery  _et al._ (2014); Démery (2015), and kinetic theory
Wang and Sperl (2016).

Exact results beyond the stationary state have been accomplished for driven
transport in lattice models, e.g., for a biased intruder in a dense crowded
environment of mobile hard-core obstacles Bénichou  _et al._ (2013a); Illien
_et al._ (2013); Bénichou  _et al._ (2013b, 2014); Illien  _et al._ (2014, )
and the complementary limit of a tracer in a dilute quenched array of
obstacles Leitmann and Franosch (2013, 2017); Leitmann  _et al._ (2018). Due
to repeated encounters with the same obstacle one finds that the nonlinear
force-dependent mobility in the stationary state becomes a nonanalytic
function of the driving force.

For the case of colloids in continuum, one infers that the mobility also
becomes a nonanalytic function in the Péclet number, signaled by
$|\text{Pe}|^{3}$ in the asymptotic expansion [Eq. (1)], since the mobility
has to be an even function of Pe.

The frequency-dependent linear mobility $\hat{\mu}(\omega)$ can be obtained
via the fluctuation-dissipation theorem from the velocity-autocorrelation
function in equilibrium, which has been calculated earlier to the same order
in the packing fraction Hanna  _et al._ (1982, 1981); Ackerson and Fleishman
(1982); Felderhof and Jones (1983a); Khair and Brady (2005)

$\displaystyle-\mathsf{i}\omega\hat{\mu}(\omega)/\mu=1-\frac{2\pi
n\sigma^{3}}{3}\frac{D_{a}}{D_{r}}\frac{1+\sqrt{-\mathsf{i}\omega\tau}}{1+\sqrt{-\mathsf{i}\omega\tau}-\mathsf{i}\omega\tau/2},$
(2)

with diffusive time scale $\tau=\sigma^{2}/D_{r}$. Here, the nonanalytic
contribution in the frequency $\omega$ reflects the well-known long-time tails
in the velocity-autocorrelation function $Z(t)\simeq-At^{-5/2}$, $A>0$, also
familiar from the Lorentz model Ernst and Weyland (1971); van Beijeren (1982);
Nieuwenhuizen  _et al._ (1986); Höfling and Franosch (2007); Franosch  _et
al._ (2010); Bauer  _et al._ (2010).

These results immediately raise the question how the nonanalytic contributions
in the nonlinear mobility emerge from a perturbative scheme and whether they
are related to the persistent correlations in equilibrium? How fast is the
nonequilibrium steady state approached from an initial equilibrium state in
comparison to the predictions of linear response? What is the intermediate
time-dependent behavior of the fluctuations connecting the short-time motion
to the drastically enhanced long-time diffusion? These questions will be
answered by solving the two-particle Smoluchowski equation for the time-
dependent dynamics in the presence of a strong force for the first time.
Thereby, we unify the previous approaches for the driven stationary state
Squires and Brady (2005); Zia and Brady (2010) and the time-dependent
equilibrium dynamics Hanna  _et al._ (1982, 1981); Ackerson and Fleishman
(1982); Felderhof and Jones (1983a) and reveal the interplay between
persistent correlations and nonequilibrium driving. The solution enables us to
fully address the time-dependent approach to the stationary state, in
principle for all moments of the displacement along the force. The two lowest
moments, the mobility and the fluctuations, will be elaborated and compared to
computer simulations.

This work is organized as follows. In Sec. II the underlying model is
introduced followed by the complete solution strategy. The main result of this
section is the self-energy encoding the dynamics between probe and the bath
particles. Readers who are primarily concerned about the results may skip this
section and jump directly to Sec. III where we determine the time-dependent
behavior in terms of the mobility and the fluctuations along the applied
force. The obtained results from the analytic solution are compared to
computer simulations and the phenomena involved in the transition from the
initial equilibrium state to the new stationary state are discussed. In Sec.
IV, the key results of this work are summarized followed by general
conclusions.

##  II Model and solution strategy

We solve for the time-dependent dynamics of a probe particle pulled by a force
in the presence of other bath particles. In first order of the density of bath
particles, the dynamics is completely encoded by the interactions of the probe
particle with a single bath particle Felderhof and Jones (1983b); Squires and
Brady (2005). Discarding inertial effects it is sufficient to consider the
two-particle Smoluchowski equation for probe and bath particle. The problem
can be expressed via the independent motion of the center of diffusion and the
relative distance between both particles. The dynamics of the relative
distance is described in the frequency domain in terms of a self-energy which
encodes all corrections to the dynamics of the probe particle due to
interactions.

###  II.1 Two-particle Smoluchowski equation

We consider a single probe particle $a$ interacting with a bath particle $b$
with bare diffusion coefficients $D_{a}$ and $D_{b}$, respectively. The
particles interact by mutual hard-core exclusion with exclusion distance
$\sigma$. At time $t=0$ an external constant force $\vec{F}$ is switched on,
pulling the probe particle and driving the system from its initial equilibrium
state into a nonequilibrium stationary state. We describe the state of the
system by the conditional probability density
$\Psi(\vec{r}_{a},\vec{r}_{b},t|\vec{r}_{a}^{\,\prime},\vec{r}_{b}^{\,\prime})$
for probe particle $a$ and bath particle $b$ to be at positions $\vec{r}_{a}$,
$\vec{r}_{b}$ at time $t$ provided they start at initial positions
$\vec{r}_{a}^{\,\prime}$, $\vec{r}_{b}^{\,\prime}$ at time $t=0$. Since the
system is initially in equilibrium, the corresponding initial condition reads
$\Psi(\vec{r}_{a},\vec{r}_{b},t=0|\vec{r}_{a}^{\,\prime},\vec{r}_{b}^{\,\prime})=\vartheta(|\vec{r}_{a}^{\,\prime}-\vec{r}_{b}^{\,\prime}|-\sigma)\delta(\vec{r}_{a}-\vec{r}_{a}^{\,\prime})\delta(\vec{r}_{b}-\vec{r}_{b}^{\,\prime})/V$,
where the Heaviside-function $\vartheta(\cdot)$ accounts for the mutual
exclusion. The limit of large box sizes $V\to\infty$ is anticipated
throughout. At time $t=0$, the force $\vec{F}$ is switched on and the
probability density evolves according to the Smoluchowski equation

$\displaystyle\partial_{t}\Psi=(D_{a}\nabla_{a}^{2}+D_{b}\nabla_{b}^{2})\Psi-\mu\vec{F}\cdot\vec{\nabla}_{a}\Psi,$
(3)

with bare mobility $\mu=D_{a}/k_{\text{B}}T$ of the probe particle and the
thermal scale $k_{\text{B}}T$. The hard-core interaction between probe and
bath particle is encoded in the no-flux boundary condition

$\displaystyle(\vec{r}_{a}-\vec{r}_{b})\cdot\Bigl{[}\mu\vec{F}\Psi-(D_{a}\vec{\nabla}_{a}-D_{b}\vec{\nabla}_{b})\Psi\Bigr{]}=0,\quad\text{for
}|\vec{r}_{a}-\vec{r}_{b}|=\sigma.$ (4)

We introduce new coordinates for the center of diffusion
$\vec{R}=(D_{b}\vec{r}_{a}+D_{a}\vec{r}_{b})/(D_{a}+D_{b})$, and the relative
distance $\vec{r}=\vec{r}_{a}-\vec{r}_{b}$. After transformation, the
Smoluchowski equation in these adapted coordinates reads

$\displaystyle\begin{split}\partial_{t}\Psi&=\frac{D_{a}D_{b}}{D_{r}}\nabla^{2}_{R}\Psi-\frac{\mu
D_{b}}{D_{r}}\vec{F}\cdot\vec{\nabla}_{R}\Psi+D_{r}\nabla^{2}_{r}\Psi-\mu\vec{F}\cdot\vec{\nabla}_{r}\Psi,\end{split}$
(5)

with the diffusion coefficient for the relative motion, $D_{r}=D_{a}+D_{b}$.
Similarly, the no-flux boundary condition transforms to

$\displaystyle\begin{split}\vec{r}\cdot\left[\mu\vec{F}\Psi-
D_{r}\vec{\nabla}_{r}\Psi\right]=0,\quad\text{for
}|\vec{r}|=\sigma.\end{split}$ (6)

The Smoluchowski equation and the no-flux boundary condition in the new
coordinates [Eq. (5) and (6)] reveal that the motion of the center of
diffusion $\vec{R}$ and the dynamics of the relative motion $\vec{r}$, are
independent. Thus, the conditional probability factorizes into a simple
Gaussian with diffusion coefficient $D_{a}D_{b}/D_{r}$ and drift
$\mu\vec{F}D_{b}/D_{r}$ for the center of diffusion $\vec{R}$, and the
conditional probability $\psi(\vec{r},t|\vec{r}\,^{\prime})$ to find relative
distances $\vec{r}$ and $\vec{r}\,^{\prime}$ at times $t$ and $0$. It fulfills
the reduced Smoluchowski equation

$\displaystyle\partial_{t}\psi=D_{r}\nabla^{2}_{r}\psi-\mu\vec{F}\cdot\vec{\nabla}_{r}\psi,$
(7)

and the initial condition is provided by the equilibrium state
$\psi(\vec{r},t=0|\vec{r}\,^{\prime})=\vartheta(|\vec{r}\,^{\prime}|-\sigma)\delta(\vec{r}-\vec{r}\,^{\prime})$.
Similarly, the no-flux boundary condition for the conditional probability
$\psi$ reads

$\displaystyle\vec{r}\cdot\left[\mu\vec{F}\psi-
D_{r}\vec{\nabla}_{r}\psi\right]=0,\qquad\text{for }|\vec{r}|=\sigma.$ (8)

It is natural to measure the force in terms of the dimensionless Péclet number
$\text{Pe}=\mu F\sigma/D_{r}$, that already appeared in Eq. (1) Squires and
Brady (2005); Zia and Brady (2010), which can be written also as
$\text{Pe}=(F\sigma/k_{\text{B}}T)D_{a}/D_{r}$ by using the Stokes-Einstein
relation. This definition encompasses the case of equal sized colloids,
$\text{Pe}=F\sigma/2k_{\text{B}}T$, as well as the Lorentz model,
$\text{Pe}=F\sigma/k_{\text{B}}T$, where the bath particles are fixed in space
($D_{b}=0$).

Our main quantity of interest is the intermediate scattering function $\langle
e^{-\mathsf{i}\vec{q}\cdot\Delta\vec{r}_{a}(t)}\rangle$ for the displacement
$\Delta\vec{r}_{a}(t)=\vec{r}_{a}(t)-\vec{r}_{a}(0)$ of the probe particle,
from which in principle all moments of the displacement can be extracted by
derivatives with respect to the wave vector $\vec{q}$. In new coordinates,
$\Delta\vec{r}_{a}(t)=\Delta\vec{R}(t)+D_{a}\Delta\vec{r}(t)/D_{r}$, the
intermediate scattering function can be expressed as

$\displaystyle\langle
e^{-\mathsf{i}\vec{q}\cdot\Delta\vec{r}_{a}(t)}\rangle=\langle
e^{-\mathsf{i}\vec{q}\cdot\Delta\vec{R}(t)}\rangle\langle
e^{-\mathsf{i}(D_{a}/D_{r})\vec{q}\cdot\Delta\vec{r}(t)}\rangle,$ (9)

where we used the independence of $\vec{R}$ and $\vec{r}$ [Eq. (5)]. Since the
dynamics of the center of diffusion is Gaussian with diffusion coefficient
$D_{a}D_{b}/D_{r}$ and drift $\mu\vec{F}D_{b}/D_{r}$, one finds immediately

$\displaystyle\langle
e^{-\mathsf{i}\vec{q}\cdot\Delta\vec{R}(t)}\rangle=e^{-\mathsf{i}\vec{q}\cdot(\mu\vec{F}D_{b}/D_{r})t}e^{-q^{2}(D_{a}D_{b}/D_{r})t}.$
(10)

The remaining task is the calculation of the intermediate scattering function
$\langle e^{-\mathsf{i}(D_{a}/D_{r})\vec{q}\cdot\Delta\vec{r}(t)}\rangle$ for
the relative motion $\Delta\vec{r}(t)$ between probe and bath particle.

###  II.2 Dynamics of the relative motion

This part contains the detailed calculation of the relative motion and
contains the heart of our analytic approach. The solution strategy is adapted
from Felderhof’s approach for the equilibrium dynamics Felderhof and Jones
(1983a).

The intermediate scattering function, viz. the characteristic function of the
relative displacement $\Delta\vec{r}(t)=\vec{r}(t)-\vec{r}(0)$,

$\displaystyle\langle
e^{-\mathsf{i}\vec{q}\cdot\Delta\vec{r}(t)}\rangle=\int\\!\mathrm{d}^{3}r\int\\!\frac{\mathrm{d}^{3}r^{\prime}}{V}\
e^{-\mathsf{i}\vec{q}\cdot(\vec{r}-\vec{r}\,^{\prime})}\psi(\vec{r},t|\vec{r}\,^{\prime}),$
(11)

allows extracting the moments by series expansion in the wave vector
$\vec{q}$. To make analytic progress, we only consider the axially symmetric
case, where the wave vector and force are aligned $\vec{q}\parallel\vec{F}$.
This enables us to determine the motion in the direction of the force, the
time-dependent motion perpendicular to the force is left for future analysis.
We perform a temporal Fourier-Laplace transform and a spatial Fourier
transform,

$\displaystyle\hat{\psi}_{q}(\vec{r},\omega):=\int_{0}^{\infty}\\!\mathrm{d}t\,e^{\mathsf{i}\omega
t}\int\\!\frac{\mathrm{d}^{3}r^{\prime}}{\sqrt{V}}\,e^{\mathsf{i}\vec{q}\cdot\vec{r}\,^{\prime}}\psi(\vec{r},t|\vec{r}\,^{\prime}).$
(12)

The transformed quantity $\hat{\psi}_{q}(\vec{r},\omega)$ is connected to the
propagator $G(q,\omega)$ via a spatial Fourier transform
$G(q,\omega):=\int\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\hat{\psi}_{q}(\vec{r},\omega)/\sqrt{V}$.
From Eq. (7) and the initial condition one derives the equation of motion
Felderhof and Jones (1983a)

$\displaystyle(-\mathsf{i}\omega-
D_{r}\nabla^{2}+\mu\vec{F}\cdot\vec{\nabla})\hat{\psi}_{q}=\frac{1}{\sqrt{V}}e^{\mathsf{i}\vec{q}\cdot\vec{r}}\vartheta(|\vec{r}|-\sigma).$
(13)

In particular, $\hat{\psi}_{q}(\vec{r},\omega)=0$ for $|\vec{r}|<\sigma$,
reflects the hard-core exclusion. Without interaction, we denote the
respective conditional probability by $\hat{\psi}_{q}^{0}$ which evolves
according to

$\displaystyle(-\mathsf{i}\omega-
D_{r}\nabla^{2}+\mu\vec{F}\cdot\vec{\nabla})\hat{\psi}_{q}^{0}=\frac{1}{\sqrt{V}}e^{\mathsf{i}\vec{q}\cdot\vec{r}}.$
(14)

The free motion allows for a plane-wave solution of the form

$\displaystyle\hat{\psi}_{q}^{0}(\vec{r},\omega)=\frac{e^{\mathsf{i}\vec{q}\cdot\vec{r}}/\sqrt{V}}{-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu
Fq},$ (15)

which is valid for all $\vec{r}\in V$ and the respective free propagator can
be read off as

$\displaystyle
G_{0}(q,\omega)=\frac{1}{-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu Fq}.$ (16)

In particular, the plane-wave solution can be written as
$\hat{\psi}_{q}^{0}(\vec{r},\omega)=G_{0}(q,\omega)e^{\mathsf{i}\vec{q}\cdot\vec{r}}/\sqrt{V}$.
To make further progress, we observe that the difference of Eqs. (13) and
(14),

$\displaystyle(-\mathsf{i}\omega-
D_{r}\nabla^{2}+\mu\vec{F}\cdot\vec{\nabla})[\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)]=0,\quad\text{for
}|\vec{r}|>\sigma,$ (17)

vanishes for terminal distance $|\vec{r}|$ larger than the exclusion distance
$\sigma$. We are interested in the forward-scattering amplitude

$\displaystyle
G(q,\omega)-G_{0}(q,\omega)=\frac{1}{\sqrt{V}}\int\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}[\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)]$
(18)

which we separate into contributions outside and inside of the region of
overlap of probe and bath particle:

$\displaystyle\begin{split}G(q,\omega)-G_{0}(q,\omega)&=\frac{1}{\sqrt{V}}\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}[\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)]+\\\
&\phantom{=}\ +\frac{1}{\sqrt{V}}\int_{|\vec{r}|<\sigma}\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}[\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)].\end{split}$
(19)

In the second term we observe that the conditional probability
$\hat{\psi}_{\vec{q}}$ vanishes for terminal positions inside the obstacle.
Hence, there we can immediately compute the integral

$\displaystyle\frac{1}{\sqrt{V}}\int_{|\vec{r}|<\sigma}\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\hat{\psi}_{q}^{0}(\vec{r},\omega)$
$\displaystyle=\frac{4\pi\sigma^{3}/3V}{-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu
Fq}.$ (20)

To calculate the first term in Eq. (19), it is advantageous to convert the
volume integral into a surface integral relying on the equations of motion. We
use the auxiliary variable
$\chi(\vec{r})=\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)$
and investigate the identity

$\displaystyle\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}(-\mathsf{i}\omega-
D_{r}\nabla^{2}+\mu\vec{F}\cdot\vec{\nabla})\chi(\vec{r})=0,$ (21)

which follows from Eq. (17). Then, for the drift, we obtain

$\displaystyle\begin{split}\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}(-\mu\vec{F}\cdot\vec{\nabla})\chi(\vec{r})&=\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
(-\mu\vec{F}\cdot\vec{\nabla})\Bigl{[}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\chi(\vec{r})\Bigr{]}+\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
\chi(\vec{r})\mu\vec{F}\cdot\vec{\nabla}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\\\
&=\int\mathrm{d}S_{\vec{r}}\,\hat{\vec{r}}\cdot\mu\vec{F}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\chi(\vec{r})+\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
\chi(\vec{r})\mu\vec{F}\cdot\vec{\nabla}e^{-\mathsf{i}\vec{q}\cdot\vec{r}},\end{split}$
(22)

where we introduced the unit vector $\hat{\vec{r}}=\vec{r}/|\vec{r}|$ and the
surface element $\mathrm{d}S_{\vec{r}}$ of the sphere of radius $\sigma$.
Similarly, for the diffusive contribution, we derive

$\displaystyle\begin{split}\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}D_{r}\nabla^{2}\chi(\vec{r})&=\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
D_{r}\vec{\nabla}\cdot\bigl{[}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\vec{\nabla}\chi(\vec{r})\bigr{]}-\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
\bigl{[}D_{r}\vec{\nabla}\chi(\vec{r})\bigr{]}\cdot\vec{\nabla}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\\\
&=-\int\mathrm{d}S_{\vec{r}}\,D_{r}\hat{\vec{r}}\cdot
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\vec{\nabla}\chi(\vec{r})-\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
D_{r}\vec{\nabla}\cdot\bigl{[}\chi(\vec{r})\
\vec{\nabla}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\bigr{]}\ +\\\ &\phantom{=}\
+\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
\chi(\vec{r})D_{r}\nabla^{2}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\\\
&=-\int\mathrm{d}S_{\vec{r}}\,D_{r}\hat{\vec{r}}\cdot
e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\vec{\nabla}\chi(\vec{r})+\int\mathrm{d}S_{\vec{r}}\,D_{r}\hat{\vec{r}}\cdot\bigl{[}\chi(\vec{r})\
\vec{\nabla}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\bigr{]}\ +\\\ &\phantom{=}\
+\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}r\
\chi(\vec{r})D_{r}\nabla^{2}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}.\end{split}$
(23)

Together, this yields

$\displaystyle\begin{split}(-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu\vec{F}\cdot\vec{q})\int_{|\vec{r}|>\sigma}\mathrm{d}^{3}\vec{r}\
\chi(\vec{r})e^{-\mathsf{i}\vec{q}\cdot\vec{r}}&=\int\mathrm{d}S_{\vec{r}}\,\hat{\vec{r}}\cdot\bigl{[}\mu\vec{F}\chi(\vec{r})-D_{r}\vec{\nabla}\chi(\vec{r})\bigr{]}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\
+\\\ &\phantom{=}\
+\int\mathrm{d}S_{\vec{r}}\,D_{r}\hat{\vec{r}}\cdot\bigl{[}\chi(\vec{r})\
\vec{\nabla}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\bigr{]}.\end{split}$ (24)

The first term on the right-hand side simplifies by the no-flux boundary
condition [Eq. (8)] and the plane wave solution $\hat{\psi}_{q}^{0}$ [Eq.
(15)]:

$\displaystyle\begin{split}\int\mathrm{d}S_{\vec{r}}\,\hat{\vec{r}}\cdot\left[\mu\vec{F}\chi-
D_{r}\vec{\nabla}\chi\right]e^{-\mathsf{i}\vec{q}\cdot\vec{r}}&=-\int\mathrm{d}S_{\vec{r}}\,e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\hat{\vec{r}}\cdot\left[\mu\vec{F}-D_{r}\vec{\nabla}\right]\hat{\psi}_{\vec{q}}^{0}(\vec{r},\omega)\\\
&=-\int\mathrm{d}S_{\vec{r}}\,e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\hat{\vec{r}}\cdot\left[\mu\vec{F}-\mathsf{i}D_{r}\vec{q}\right]\hat{\psi}_{\vec{q}}^{0}(\vec{r},\omega)\\\
&=\frac{-1/\sqrt{V}}{-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu\vec{F}\cdot\vec{q}}\int\mathrm{d}S_{\vec{r}}\,\hat{\vec{r}}\cdot\left[\mu\vec{F}-\mathsf{i}D_{r}\vec{q}\right]=0,\end{split}$
(25)

where in the last line the flux integral of a constant vector vanishes.
Collecting results, we obtain an expression for the forward scattering
amplitude in terms of a surface integral of the difference of the conditional
probabilities:

$\displaystyle\begin{split}G(q,\omega)-G_{0}(q,\omega)&=G_{0}(q,\omega)\biggl{\\{}-\frac{4\pi\sigma^{3}}{3V}-\frac{1}{\sqrt{V}}\int\mathrm{d}S_{\vec{r}}\,\mathsf{i}D_{r}\vec{q}\cdot\hat{\vec{r}}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\bigl{[}\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)\bigr{]}\biggr{\\}}.\end{split}$
(26)

In order to determine the conditional probabilities, we return to Eq. (17).
For the solution of the homogeneous part of Eq. (17), we use the imaginary
“gauge transformation”
$\hat{X}(\vec{r},\omega)=[\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)]e^{-\mu\vec{F}\cdot\vec{r}/2D_{r}}$
such that the drift term is absorbed in the Laplacian by completing the
square, $\vec{\nabla}\to\vec{\nabla}-\mu\vec{F}/2D_{r}$ Nelson and Shnerb
(1998); Squires and Brady (2005). This step will have crucial implications for
the results at finite forces. Then, the new quantity $\hat{X}$ fulfills the 3d
source-free Helmholtz equation

$\displaystyle(\kappa^{2}-\nabla^{2})\hat{X}(\vec{r},\omega)=0,\quad\text{for
}|\vec{r}|>\sigma,$ (27)

with complex wavenumber
$\kappa^{2}\sigma^{2}=-\mathsf{i}\omega\tau+(\text{Pe}/2)^{2}$ and diffusive
time scale $\tau=\sigma^{2}/D_{r}$. We write the general axially symmetric
solution in the form

$\displaystyle\hat{\psi}_{q}(\vec{r},\omega)-\hat{\psi}^{0}_{q}(\vec{r},\omega)=\frac{e^{\mu\vec{F}\cdot\vec{r}/2D_{r}}/\sqrt{V}}{-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu
Fq}\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{k}_{\ell}(\kappa
r)}{\text{k}_{\ell}(\kappa\sigma)}\text{P}_{\ell}(\cos\vartheta),$ (28)

for $|\vec{r}|>\sigma$, with modified spherical Bessel functions of the second
kind $\text{k}_{\ell}(\cdot)$ and Legendre polynomials
$\text{P}_{\ell}(\cdot)$.

The expansion coefficients $a_{\ell}$ for the different angular channels
$\ell$ have to be determined via the no-flux boundary condition [Eq. (8)]:

$\displaystyle\begin{split}\left[\mu F\cos\vartheta-
D_{r}\frac{\partial}{\partial
r}\right](\hat{\psi}_{q}-\hat{\psi}_{q}^{0})\biggr{|}_{r=\sigma}&=-\left[\mu
F\cos\vartheta-D_{r}\frac{\partial}{\partial
r}\right]\hat{\psi}_{q}^{0}\biggr{|}_{r=\sigma}.\end{split}$ (29)

With the explicit solutions [Eqs. (15) and (28)], we find that the expansion
coefficients $a_{\ell}$ have to fulfill

$\displaystyle\left[\text{Pe}\,\eta-\sigma\frac{\partial}{\partial
r}\right]e^{\text{Pe}\,\eta
r/2\sigma}\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{k}_{\ell}(\kappa
r)}{\text{k}_{\ell}(\kappa\sigma)}\text{P}_{\ell}(\eta)\biggr{|}_{r=\sigma}=$
$\displaystyle-\left[\text{Pe}\,\eta-\sigma\frac{\partial}{\partial
r}\right]e^{\mathsf{i}q\eta r}\biggr{|}_{r=\sigma},$ (30)

where we abbreviated $\eta=\cos(\vartheta)$ and used the Péclet number
$\text{Pe}=\mu F\sigma/D_{r}$. Performing the derivatives the relation can be
expressed as

$\displaystyle\begin{split}e^{-\text{Pe}\,\eta/2}(\text{Pe}-\mathsf{i}q\sigma)\eta
e^{\mathsf{i}q\sigma\eta}&=\sum_{\ell=0}^{\infty}a_{\ell}\biggl{[}\frac{\kappa\sigma\text{k}_{\ell}^{\prime}(\kappa\sigma)}{\text{k}_{\ell}(\kappa\sigma)}-\frac{\text{Pe}}{2}\eta\biggr{]}\text{P}_{\ell}(\eta)\\\
&=\sum_{\ell=0}^{\infty}\biggl{[}\frac{\kappa\sigma\text{k}_{\ell}^{\prime}(\kappa\sigma)}{\text{k}_{\ell}(\kappa\sigma)}a_{\ell}\text{P}_{\ell}(\eta)-\frac{\text{Pe}}{2}a_{\ell}\biggl{(}\frac{\ell+1}{2\ell+1}\text{P}_{\ell+1}(\eta)+\frac{\ell}{2\ell+1}\text{P}_{\ell-1}(\eta)\biggr{)}\biggr{]}\\\
&=\sum_{\ell=0}^{\infty}\biggl{[}\frac{\kappa\sigma\text{k}_{\ell}^{\prime}(\kappa\sigma)}{\text{k}_{\ell}(\kappa\sigma)}a_{\ell}-\frac{\text{Pe}}{2}\frac{\ell}{2\ell-1}a_{\ell-1}-\frac{\text{Pe}}{2}\frac{\ell+1}{2\ell+3}a_{\ell+1}\biggr{]}\text{P}_{\ell}(\eta),\end{split}$
(31)

where we used the recursion formula
$(\ell+1)\text{P}_{\ell+1}(\eta)=(2\ell+1)\eta\text{P}_{\ell}(\eta)-\ell\text{P}_{\ell-1}(\eta)$.
Using the orthogonality relation of the Legendre polynomials,
$\int_{-1}^{1}\mathrm{d}\eta\
\text{P}_{\ell}(\eta)\text{P}_{\ell^{\prime}}(\eta)=2\delta_{\ell\ell^{\prime}}/(2\ell+1)$,
this becomes a tridiagonal matrix equation:

$\displaystyle\frac{\kappa\sigma\text{k}_{\ell}^{\prime}(\kappa\sigma)}{\text{k}_{\ell}(\kappa\sigma)}a_{\ell}-\frac{\text{Pe}}{2}\frac{\ell}{2\ell-1}a_{\ell-1}-\frac{\text{Pe}}{2}\frac{\ell+1}{2\ell+3}a_{\ell+1}=b_{\ell}.$
(32)

The inhomogeneity $b_{\ell}$ is defined as

$\displaystyle b_{\ell}$
$\displaystyle=\frac{2\ell+1}{2}(\text{Pe}-\mathsf{i}q\sigma)\int_{-1}^{1}\mathrm{d}\eta\
\eta\text{P}_{\ell}(\eta)e^{-\text{Pe}\,\eta/2}e^{\mathsf{i}q\sigma\eta}.$
(33)

The remaining integral can be be solved using the Rayleigh identity
$e^{-\mathsf{i}z\cos\vartheta}=\sum_{\ell=0}^{\infty}(-\mathsf{i})^{\ell}(2\ell+1)\text{j}_{\ell}(z)\text{P}_{\ell}(\cos\vartheta),$
with spherical Bessel function $\text{j}_{\ell}(\cdot)$ and relying again on
the orthogonality relation of the Legendre polynomials:

$\displaystyle
b_{\ell}=-(2\ell+1)(q\sigma+\mathsf{i}\text{Pe})\mathsf{i}^{\ell}\text{j}_{\ell}^{\prime}(q\sigma+\mathsf{i}\text{Pe}/2),$
(34)

where the prime indicates a derivative. For the special case of vanishing
wavenumber $q=0$ and the stationary state $\omega=0$, we recover the
tridiagonal matrix derived in Ref. Squires and Brady (2005).

For the forward scattering amplitude [Eq. (26)] we calculate the integral over
the spherical surface with radius $r=\sigma$:

$\displaystyle\begin{split}\frac{1}{\sqrt{V}}\int\mathrm{d}S_{\vec{r}}\,\mathsf{i}D_{r}\vec{q}\cdot\hat{\vec{r}}e^{-\mathsf{i}\vec{q}\cdot\vec{r}}\bigl{[}\hat{\psi_{q}}(\vec{r},\omega)-\hat{\psi}_{q}^{0}(\vec{r},\omega)\bigr{]}&=\frac{2\pi\mathsf{i}D_{r}\sigma^{2}q/V}{-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu
Fq}\times\\\ &\phantom{=}\
\times\sum_{\ell=0}^{\infty}a_{\ell}\int_{-1}^{1}\mathrm{d}\eta\,\eta
e^{-\mathsf{i}q\sigma\eta}e^{\mu
F\sigma\eta/2D_{r}}\text{P}_{\ell}(\eta),\end{split}$ (35)

The integral on the right-hand side can be evaluated again by inserting the
Rayleigh identity. As a result, we find

$\displaystyle\begin{split}\int_{-1}^{1}\mathrm{d}\eta\,\eta
e^{-\mathsf{i}q\sigma\eta}e^{\eta\text{Pe}/2}\text{P}_{\ell}(\eta)=\frac{\partial}{\partial(\text{Pe}/2)}\int_{-1}^{1}\mathrm{d}\eta\,e^{-\mathsf{i}q\sigma\eta}e^{\eta\text{Pe}/2}\text{P}_{\ell}(\eta)=2\mathsf{i}(-\mathsf{i})^{\ell}\,\text{j}_{\ell}^{\prime}(q\sigma+\mathsf{i}\text{Pe}/2),\end{split}$
(36)

Collecting results, one finds

$\displaystyle
G(q,\omega)-G_{0}(q,\omega)=-\frac{4\pi\sigma^{3}}{3V}G_{0}(q,\omega)+\frac{4\pi
D_{r}\sigma^{2}q}{V}G_{0}(q,\omega)^{2}\sum_{\ell=0}^{\infty}a_{\ell}(-\mathsf{i})^{\ell}\text{j}_{\ell}^{\prime}(q\sigma+\mathsf{i}\text{Pe}/2).$
(37)

The first term on the right-hand side merely reflects that the free propagator
$G_{0}(q,\omega)=(-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}\mu Fq)^{-1}$ allows
for particles starting also inside of the obstacle, thereby renormalizing the
residue of the perturbed propagator. We may safely drop this term.
Conventionally, the result is expressed in terms of a self-energy
$\Sigma(q,\omega)$ via the Dyson equation $G=G_{0}+G_{0}\Sigma G$. To first
order in the density, the self-energy is merely proportional to the number
density $n=N/V$ of the bath particles, and we obtain the self-energy as our
main result of the analytic calculations

$\displaystyle\Sigma(q,\omega)$ $\displaystyle=4\pi
nD_{r}\sigma^{2}q\sum_{\ell=0}^{\infty}a_{\ell}(-\mathsf{i})^{\ell}\text{j}_{\ell}^{\prime}(q\sigma+\mathsf{i}\text{Pe}/2).$
(38)

The self-energy encodes the density-induced corrections of all moments of the
relative motion along the force. The leading factor $q$ reflects the particle-
conservation law.

Without external driving, $\text{Pe}=0$, the tridiagonal matrix for the
coefficients $a_{\ell}$ [Eq. (32)] becomes diagonal yielding

$\displaystyle
a_{\ell}=-(2\ell+1)q\sigma\mathsf{i}^{\ell}\text{j}_{\ell}^{\prime}(q\sigma)\frac{\text{k}_{\ell}(\kappa\sigma)}{\kappa\sigma\text{k}_{\ell}^{\prime}(\kappa\sigma)}.$
(39)

Thus, we recover the known result for the self-energy for equilibrium Hanna
_et al._ (1982, 1981); Ackerson and Fleishman (1982); Felderhof and Jones
(1983a):

$\displaystyle\Sigma(q,\omega)$ $\displaystyle=-4\pi
nD_{r}\sigma^{3}q^{2}\sum_{\ell=0}^{\infty}(2\ell+1)\bigl{[}\text{j}_{\ell}^{\prime}(q\sigma)\bigr{]}^{2}\frac{\text{k}_{\ell}(\kappa\sigma)}{\kappa\sigma\text{k}_{\ell}^{\prime}(\kappa\sigma)}.$
(40)

In the case of driving, only terms of order $\ell=\mathcal{O}(\text{Pe})$
significantly contribute for small wave numbers and the matrix [Eq. (32)] may
be safely truncated for numerical evaluation. The matrix inversion does not
generate nonanalytic behavior and one infers that the coefficients
$a_{\ell}\equiv a_{\ell}(q\sigma,\kappa\sigma,\text{Pe})$ should be analytic
functions in the arguments. In particular, the matrix is suited for a
perturbative approach for small Péclet numbers, in particular, the linear
response results can be derived.

The matrix in Eq. (32) determines the coefficients $a_{\ell}$ and thereby the
complete solution of the self-energy [Eq. (38)]. Up to matrix inversion and a
temporal Fourier back-transform, which have to be implemented numerically, we
have elaborated a complete time-dependent analytic solution for the relative
motion of the probe particle in the presence of bath particles for the dilute
case.

##  III Cumulants of the displacement

The intermediate scattering function generates the moments of the displacement
along the field by taking derivatives with respect to the wavenumber. In this
section we elaborate explicitly how the mean displacement, respectively the
time-dependent nonlinear mobility, and the fluctuations along the field can be
calculated from the self-energy of the relative motion.

The starting point are the intermediate scattering functions for the
displacement [Eqs. (9) and (10)] where we put the wave vector along the field,
$\vec{q}\parallel\vec{F}$, which is chosen as the $z$ direction. The
corresponding cumulant generating function is obtained by taking the
logarithm:

$\displaystyle\begin{split}\ln\langle e^{-\mathsf{i}q\Delta
z_{a}(t)}\rangle&=\ln\langle e^{-\mathsf{i}q\Delta Z(t)}\rangle+\ln\langle
e^{-\mathsf{i}(D_{a}q/D_{r})\Delta z(t)}\rangle\\\ &=-\mathsf{i}q(\mu
FD_{b}/D_{r})t-q^{2}(D_{a}D_{b}/D_{r})t+\ln\langle
e^{-\mathsf{i}(D_{a}q/D_{r})\Delta z(t)}\rangle.\end{split}$ (41)

This formula establishes the connection between the relative motion $z(t)$ and
the displacement of the probe particle $z_{a}(t)$.

###  III.1 Average motion

First, we discuss the time-dependent nonlinear mobility defined by

$\displaystyle\mu(t,\text{Pe}):=\frac{1}{F}\frac{\mathrm{d}}{\mathrm{d}t}\langle\Delta
z_{a}(t)\rangle,$ (42)

characterizing the mean motion of the probe particle, upon switching on the
force in the $z$ direction. From the cumulant generating function of the
displacement $\Delta z_{a}(t)$ [Eq. (41)], the mean-displacement
$\langle\Delta z_{a}(t)\rangle$ is obtained as

$\displaystyle\begin{split}\langle\Delta
z_{a}(t)\rangle&=\mathsf{i}\frac{\partial}{\partial q}\ln\langle
e^{-\mathsf{i}q\Delta z_{a}(t)}\rangle\Bigr{|}_{q=0}=\mu
FtD_{b}/D_{r}+\frac{\mathsf{i}D_{a}}{D_{r}}\frac{\partial}{\partial(D_{a}q/D_{r})}\ln\langle
e^{-\mathsf{i}(D_{a}q/D_{r})\Delta z(t)}\rangle\Bigr{|}_{q=0}\\\ &=\mu
FtD_{b}/D_{r}+\frac{D_{a}}{D_{r}}\langle\Delta z(t)\rangle.\end{split}$ (43)

The contribution from the relative motion can be obtained from the
intermediate scattering function [Eq. (11)] as the first derivative
$\partial/\partial q|_{q=0}$. The dynamics of the relative distance $\Delta
z(t)$ between probe and bath particle is contained in the propagator

$\displaystyle
G(q,\omega)=G_{0}(q,\omega)+G_{0}(q,\omega)^{2}\Sigma(q,\omega)+\mathcal{O}(n)^{2},$
(44)

with free propagator
$G_{0}(q,\omega)=(-\mathsf{i}\omega+D_{r}q^{2}+\mathsf{i}q\mu F)^{-1}$ and
self-energy $\Sigma(q,\omega)$ [Eq. (38)]. Taking the $q$-derivative leads to
the appearance of the frequency dependent expression
$\mathcal{L}\\{\langle\Delta
z(t)\rangle\\}(\omega):=\int_{0}^{\infty}\mathrm{d}t\ e^{\mathsf{i}\omega
t}\langle\Delta z(t)\rangle$. It is obtained by considering the propagator
$G(q,\omega)$ [Eq. (44)]:

$\displaystyle\begin{split}\mathcal{L}\\{\langle\Delta
z(t)\rangle\\}(\omega)&=\mathsf{i}\frac{\partial}{\partial
q}G(q,\omega)\Big{|}_{q=0}=\biggl{[}\mathsf{i}\frac{\partial G_{0}}{\partial
q}+\mathsf{i}G_{0}^{2}\frac{\partial\Sigma}{\partial q}\biggr{]}_{q=0}\\\
&=\frac{\mu F}{(-\mathsf{i}\omega)^{2}}+2\pi n\sigma^{3}\frac{\mu
F}{(-\mathsf{i}\omega)^{2}}\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{i}_{\ell}^{\prime}(\text{Pe}/2)}{\text{Pe}/2},\end{split}$
(45)

where, in the second line, we used the relation $D_{r}/\sigma=\mu F/\text{Pe}$
and expressed the spherical Bessel functions $\text{j}_{\ell}(\cdot)$ in terms
of
$(-\mathsf{i})\text{i}^{\prime}_{\ell}(\text{Pe}/2)=(-\mathsf{i})^{\ell}\text{j}^{\prime}_{\ell}(\mathsf{i}\text{Pe}/2)$
with modified spherical Bessel function $\text{i}_{\ell}(\cdot)$. The
coefficients $a_{\ell}=a_{\ell}(q\sigma,\kappa\sigma,\text{Pe})$ are to be
evaluated for vanishing wavenumber $q=0$. In the above equation, the first
term corresponds to the free motion and the second encodes the density-induced
response. Then, the mean displacement of the probe particle along the field,
$\langle\Delta z_{a}(t)\rangle$, is determined via Eq. (43) leading to

$\displaystyle\begin{split}\mathcal{L}\\{\langle\Delta
z_{a}(t)\rangle\\}(\omega)&=\frac{D_{b}}{D_{r}}\frac{\mu
F}{(-\mathsf{i}\omega)^{2}}+\frac{D_{a}}{D_{r}}\mathcal{L}\\{\langle\Delta
z(t)\rangle\\}(\omega)\\\ &=\frac{\mu F}{(-\mathsf{i}\omega)^{2}}+2\pi
n\sigma^{3}\frac{\mu
F}{(-\mathsf{i}\omega)^{2}}\frac{D_{a}}{D_{r}}\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{i}_{\ell}^{\prime}(\text{Pe}/2)}{\text{Pe}/2}.\end{split}$
(46)

Figure 1:  Time-dependent approach of the mobility $\mu(t,\text{Pe})$ to its
stationary-state for different Péclet numbers, Pe, and the linear response,
$\text{Pe}=0$. Solid lines represent the analytic solution and symbols
correspond to Brownian-dynamics simulations for equal sized colloids
($D_{b}=D_{a}$). Dashed lines indicate negative values. The black dashed-
dotted lines show the long-time behavior
${\sim-t^{-3/2}\exp(-\text{Pe}^{2}D_{a}t/2\sigma^{2})}$ for the approach to
the stationary state.

By a one-sided Fourier transform of the time-dependent nonlinear mobility [Eq.
(42)], we obtain the nonlinear frequency-dependent mobility
$\hat{\mu}(\omega,\text{Pe}):=-\mathsf{i}\omega\mathcal{L}\\{\langle\Delta
z_{a}(t)\rangle\\}(\omega)/F$ with

$\displaystyle\begin{split}-\mathsf{i}\omega\hat{\mu}(\omega,\text{Pe})/\mu=1+2\pi
n\sigma^{3}\frac{D_{a}}{D_{r}}\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{i}_{\ell}^{\prime}(\text{Pe}/2)}{\text{Pe}/2}\end{split}$
(47)

For the time-dependent case, we perform a numerical Fourier inversion of Eq.
(47) and compare to event-driven Brownian dynamics simulations at low density
[Fig. 1]. Initially, the probe particle experiences only the drag of the pure
solvent, $\mu(t\to 0,\text{Pe})=\mu$. Then the dynamics slows down due to the
interaction with the bath particles and approaches its stationary mobility
$\mu(\text{Pe}):=\mu(t\to\infty,\text{Pe})$. Its dependence on the Péclet
number is discussed in Sec. III.3 below.

In linear response, $\text{Pe}=0$, we recover the equilibrium solution [Eq.
(49)] as anticipated by the fluctuation-dissipation theorem. In particular,
for the equilibrium case, the approach to the stationary state is completely
monotone and the time dependence can be worked out explicitly by Laplace
inversion of Eq. (49) Hanna  _et al._ (1981), leading to

$\displaystyle\mu(t)/\mu=1-\frac{2\pi
n\sigma^{3}}{3}\frac{D_{a}}{D_{r}}\bigl{\\{}1-\cos(2t/\tau)[1-2\text{C}(\sqrt{4t/\pi\tau})]-\sin(2t/\tau)[1-2\text{S}(\sqrt{4t/\pi\tau})]\bigr{\\}},$
(48)

with Fresnel integrals $\text{S}(x)=\int_{0}^{x}\mathrm{d}u\sin(\pi u^{2}/2)$
and $\text{C}(x)=\int_{0}^{x}\mathrm{d}u\cos(\pi u^{2}/2)$ Olver  _et al._
(2010); DLMF . Inserting the asymptotic expansions for the Fresnel integrals
for large arguments, we recover the algebraic approach of $\sim t^{-3/2}$ to
the stationary mobility which is due to repeated collisions of the probe
particle with the same bath particle:

$\displaystyle\mu(t)/\mu=1-\frac{2\pi
n\sigma^{3}}{3}\frac{D_{a}}{D_{r}}\Bigl{\\{}1-\frac{1}{4\sqrt{\pi}}(t/\tau)^{-3/2}+\mathcal{O}(t^{-7/2})\Bigr{\\}},\quad
t\to\infty.$ (49)

The persistent memory in the system also emerges as a long-time tail of the
form $\sim t^{-5/2}$ in the velocity-autocorrelation function,
$k_{\text{B}}T\mathrm{d}\mu(t)/\mathrm{d}t$.

For finite driving, expansion of Eq. (47) yields

$\displaystyle-\mathsf{i}\omega\hat{\mu}(\omega,\text{Pe})/\mu=$
$\displaystyle 1-2\pi
n\sigma^{3}\frac{D_{a}}{D_{r}}\Big{[}1-\frac{\text{Pe}^{2}}{120}+\frac{\text{Pe}^{2}}{8}\kappa\sigma-\frac{1}{2}(\kappa\sigma)^{2}+\frac{1}{2}(\kappa\sigma)^{3}+{\cal
O}(\text{Pe}^{4},\kappa^{2}\text{Pe}^{2},\kappa^{4})\Big{]}$ (50)

which is analytic in the complex wavenumber $\kappa$ and the Péclet number.
Moreover only even terms in Pe appear. The nonanalytic dependence on the
frequency and the Péclet number arises via the square root
$\kappa\sigma=\sqrt{-\mathsf{i}\omega\tau+(\text{Pe}/2)^{2}}$. For vanishing
forces $\kappa\sigma=\sqrt{-\mathsf{i}\omega\tau}$ and the term ${\cal
O}(\kappa)^{3}$ generates the long-time tail in Eq. (49). For small but finite
forces, the singularity is shifted in the complex plane and this term yields
the long-time behavior in the time domain

$\displaystyle\frac{\mu(t,\text{Pe})-\mu(\text{Pe})}{2\pi
n\sigma^{3}D_{a}/3D_{r}}=\frac{1}{4\sqrt{\pi}}(t/\tau)^{-3/2}\exp(-\text{Pe}^{2}t/4\tau),\quad\text{Pe}\to
0\quad t\to\infty,$ (51)

Therefore the initial decay becomes more rapid and the long-time tail is
followed only for small forces up to some driving-dependent crossover time
$\tau_{F}:=\tau/\text{Pe}^{2}$, where the tail is decorated by a decaying
exponential. Yet, there is a second non-analytic contribution at finite
driving ${\cal O}(\text{Pe}^{2}\kappa)$ (see also Ref. Leitmann  _et al._
(2018) for the corresponding lattice case) which yields in the temporal domain

$\displaystyle\frac{\mu(t,\text{Pe})-\mu(\text{Pe})}{2\pi
n\sigma^{3}D_{a}/3D_{r}}=-\frac{1}{8\sqrt{\pi}}\,\text{Pe}^{2}(t/\tau)^{-1/2}\exp(-\text{Pe}^{2}t/4\tau),\quad
t\to\infty,\quad\text{Pe}\to 0$ (52)

Comparing expressions in Eqs. (51),(52) reveals that for times
$t\gtrsim\tau_{F}=\tau/\text{Pe}^{2}$ the latter decays more slowly and is the
relevant one. Therefore, the divergent time scale $\tau_{F}$ in the problem
separates two different regimes. For $t\lesssim\tau_{F}$ the linear response
prediction remains qualitatively correct, while for $t\gtrsim\tau_{F}$ the
non-equilibirium driving dominates. The second regime also explains the
counterintuitive sign change in Fig. 1 such that in the terminal regime the
velocity of the probe particle speeds up again for small and intermediate
Péclet numbers. For large Péclet number the approach to the stationary state
becomes monotonic again.

The simulations quantitatively confirm the theory for a low density of bath
particles (see Appendix A for simulation details). Note that the simulations
were performed at a finite density of $n\sigma^{3}=0.01$ [Fig. 1].
Nevertheless, the agreement for the approach to the stationary behavior of the
mobility extends to relative order $10^{-3}$. In particular, for equilibrium,
we observe the persistent power-law correlations due to repeated interactions
over roughly one decade in time.

###  III.2 Fluctuations along the force

The next interesting quantity for the motion of the probe particle are the
time-dependent fluctuations around the drift motion,

$\displaystyle\text{Var}_{z}(t):=\langle\Delta
z_{a}(t)^{2}\rangle-\langle\Delta z_{a}(t)\rangle^{2},$ (53)

which is the second cumulant of the fluctuating probe displacement $\Delta
z_{a}(t)$. From the cumulant generating function [Eq. (41)], it is obtained
via

$\displaystyle\begin{split}\text{Var}_{z}(t)&=-\frac{\partial^{2}}{\partial
q^{2}}\ln\langle e^{-\mathsf{i}q\Delta
z_{a}(t)}\rangle=2(D_{a}D_{b}/D_{r})t-\frac{D_{a}^{2}}{D_{r}^{2}}\frac{\partial^{2}}{\partial(D_{a}q/D_{r})^{2}}\ln\langle
e^{-\mathsf{i}(D_{a}q/D_{r})\Delta z(t)}\rangle\\\
&=2(D_{a}D_{b}/D_{r})t+\frac{D_{a}^{2}}{D_{r}^{2}}\bigl{[}\langle\Delta
z(t)^{2}\rangle-\langle\Delta z(t)\rangle^{2}\bigr{]},\end{split}$ (54)

where, in the second line, we already inserted the explicit expression for the
second cumulant of the relative displacement $\Delta z(t)$. First, we
calculate the second moment along the field in the frequency-domain:

$\displaystyle\begin{split}\mathcal{L}\\{\langle\Delta
z(t)^{2}\rangle\\}(\omega)&=-\frac{\partial^{2}}{\partial
q^{2}}G(q,\omega)\biggr{|}_{q=0}=\biggl{[}-\frac{\partial^{2}G_{0}}{\partial
q^{2}}-4G_{0}\frac{\partial G_{0}}{\partial q}\frac{\partial\Sigma}{\partial
q}-G_{0}^{2}\frac{\partial^{2}\Sigma}{\partial q^{2}}\biggr{]}_{q=0}\\\
&=\frac{2D_{r}}{(-\mathsf{i}\omega)^{2}}+\frac{2(\mu
F)^{2}}{(-\mathsf{i}\omega)^{3}}+2\pi n\sigma^{3}\frac{4(\mu
F)^{2}}{(-\mathsf{i}\omega)^{3}}\Delta\hat{\mu}(\omega,\text{Pe})+2\pi
n\sigma^{3}\frac{4D_{r}}{(-\mathsf{i}\omega)^{2}}\Delta\hat{R}(\omega,\text{Pe}).\end{split}$
(55)

where we introduce auxiliary functions

$\displaystyle\Delta\hat{\mu}(\omega,\text{Pe})$
$\displaystyle:=\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{i}_{\ell}^{\prime}(\text{Pe}/2)}{\text{Pe}/2},$
(56) $\displaystyle\Delta\hat{R}(\omega,\text{Pe})$
$\displaystyle:=\sum_{\ell=0}^{\infty}\Bigl{[}\frac{\mathsf{i}}{\sigma}\frac{\partial
a_{\ell}}{\partial
q}\text{i}_{\ell}^{\prime}(\text{Pe}/2)+a_{\ell}\text{i}_{\ell}^{\prime\prime}(\text{Pe}/2)\Bigr{]}.$
(57)

In particular, with the auxiliary function for the mobility,
$\Delta\hat{\mu}(\omega,\text{Pe})$, the frequency-dependent nonlinear
mobility [Eq. (47)] can be written as

$\displaystyle-\mathsf{i}\omega\hat{\mu}(\omega,\text{Pe})/\mu=1+2\pi
n\sigma^{3}\frac{D_{a}}{D_{r}}\Delta\hat{\mu}(\omega,\text{Pe}).$ (58)

In the time-domain, the mean-square displacement then follows as

$\displaystyle\begin{split}\langle\Delta z(t)^{2}\rangle&=2D_{r}t+(\mu
F)^{2}t^{2}+8\pi n\sigma^{3}D_{r}\int_{0}^{t}\mathrm{d}t^{\prime}\
\mathcal{L}^{-1}\biggl{\\{}\frac{\Delta\hat{R}(\omega,\text{Pe})}{(-\mathsf{i}\omega)}+\frac{(\mu
F)^{2}}{D_{r}}\frac{\Delta\hat{\mu}(\omega,\text{Pe})}{(-\mathsf{i}\omega)^{2}}\biggr{\\}}(t^{\prime}),\end{split}$
(59)

where we introduced the inverse Laplace transform
$\mathcal{L}^{-1}\\{\cdot\\}(t)$. To calculate the square of the mean
displacement $\langle\Delta z(t)\rangle$, we first express the frequency-
dependent mean displacement given in Eq. (45) in the time-domain leading to

$\displaystyle\langle\Delta z(t)\rangle=\mu Ft+2\pi n\sigma^{3}\mu
F\int_{0}^{t}\mathrm{d}t^{\prime}\
\mathcal{L}^{-1}\biggl{\\{}\frac{\Delta\hat{\mu}(\omega,\text{Pe})}{(-\mathsf{i}\omega)}\biggr{\\}}(t^{\prime}).$
(60)

Then, the square of the mean-displacement to first order in the density is
given by

$\displaystyle\langle\Delta z(t)\rangle^{2}=(\mu F)^{2}t^{2}+4\pi
n\sigma^{3}(\mu F)^{2}t\int_{0}^{t}\mathrm{d}t^{\prime}\
\mathcal{L}^{-1}\biggl{\\{}\frac{\Delta\hat{\mu}(\omega,\text{Pe})}{(-\mathsf{i}\omega)}\biggr{\\}}(t^{\prime})+\mathcal{O}(n^{2}).$
(61)

Collecting results, the time-dependent diffusion coefficient along the field,
$D_{z}(t,\text{Pe}):=(1/2)\mathrm{d}\text{Var}_{z}(t)/\mathrm{d}t$, can be
written as

$\displaystyle\begin{split}D_{z}(t,\text{Pe})&=D_{a}+2\pi
n\sigma^{3}\frac{D_{a}^{2}}{D_{r}}\mathcal{L}^{-1}\biggl{\\{}\frac{2\Delta\hat{R}(\omega,\text{Pe})}{(-\mathsf{i}\omega)}+\frac{\text{Pe}^{2}}{(-\mathsf{i}\omega)}\frac{\partial\Delta\hat{\mu}(\omega,\text{Pe})}{\partial(-\mathsf{i}\omega\tau)}\biggr{\\}}(t),\end{split}$
(62)

where we used the relation
$t\mathcal{L}^{-1}\\{(\cdot)\\}(t)=\mathcal{L}^{-1}\\{\partial(\cdot)/\partial(-\mathsf{i}\omega)\\}(t)$
of the Laplace transform.

Figure 2:  Time-dependent diffusion coefficient $D_{z}(t)$ for different
strength of the driving, $\text{Pe}=\mu F\sigma/D_{r}$. Lines correspond to
the analytic solution and symbols represent Brownian-dynamics simulation for
equal sized colloids ($D_{b}=D_{a}$). Inset: Zoom of the same quantity for the
smaller Péclet numbers.

For small times the probe particle does not interact with the bath particles
and the diffusion coefficient is given by the bare one $D_{z}(t\to
0,\text{Pe})=D_{a}$. Only in the equilibrium case, the diffusion coefficient
decreases monotonically [see Fig. 2] to its stationary-state value
$D_{z}^{\text{eq}}/D=1-(D_{a}/D_{r})2\pi n\sigma^{3}/3$.

For any finite driving, a minimum of least diffusivity emerges at intermediate
times such that the growth of the fluctuations speeds up again until the
stationary diffusion coefficient is reached [Fig. 2 (inset)]. With increasing
driving the time of least diffusivity becomes smaller and smaller. The time-
dependent growth becomes arbitrarily large as the Péclet number is increased,
even at small densities. A growing time-dependent diffusion coefficient is a
fingerprint of transient superdiffusive behavior.

The superdiffusion becomes more and more pronounced upon increasing the force
beyond values that are feasible for numerical implementation of our analytical
solution ($\text{Pe}\gtrsim 10^{3}$). Nevertheless, in simulations the regime
of these large Péclet numbers can still be accessed [Fig. 3]. The window of
superdiffusion opens for stronger driving and the local exponent for the
variance

$\displaystyle\alpha(t):=\frac{\mathrm{d}\\!\ln[\text{Var}_{z}(t)]}{\mathrm{d}\\!\ln(t)}=\frac{2D_{z}(t,\text{Pe})t}{\text{Var}_{z}(t)}$
(63)

approaches a value of $3$ at intermediate times. The long-time behavior is for
all forces diffusive with a strongly enhanced diffusion coefficient
$D_{z}(\text{Pe})/D\simeq\pi n\sigma^{3}\text{Pe}[\ln(2)-1/4]/6$, for
$\text{Pe}\to\infty$ Zia and Brady (2010).

Let us rationalize the superdiffusive behavior in terms of an asymptotic model
for the limit of strong driving. Therefore we adapt an earlier asymptotic
model valid for the driven lattice Lorentz model Leitmann and Franosch (2017).
Here, the probe particle’s motion is dominated by the drift with constant
velocity $\mu F$ along the force until it hits a bath particle for the first
time and then slowly slides along its surface. In this time regime, the motion
becomes essentially one-dimensional such that the probability distribution for
the relative motion, $\mathbb{P}(\Delta z,t)$, consists of freely moving
particles with fixed velocity or particles that are transiently blocked by
bath particles. The free path lengths are exponentially distributed, since at
low density the positions of the bath particles are independent, and the
probability distribution can be estimated directly to

$\displaystyle\mathbb{P}(\Delta z,t)=\delta(\Delta z-vt)e^{-\Delta
z/l_{*}}+(1/l_{*})e^{-\Delta z/l_{*}}\vartheta(vt-\Delta z),$ (64)

where $l_{*}=1/n\Sigma_{*}$ denotes a mean-free path length. The first term
corresponds to the freely moving probe (with $v=\mu F$), while the second term
accounts for the blocked probe particles.

From the distribution of the displacements [Eq. (64)], it is straight forward
to calculate the mean and the mean-square displacement. Then, the growth of
the fluctuations for the probe particle are obtained from Eq. (54), where we
discarded the diffusive contribution since we are only interested in the drift
motion. As a result, the distribution of the displacements implies a strong
initial growth of the variance with

$\displaystyle\text{Var}_{z}(t)=\Bigl{(}\frac{D_{a}}{D_{r}}\Bigr{)}^{2}\frac{(\mu
Ft)^{3}}{3l_{*}},$ (65)

which is also confirmed by Brownian dynamcis simulation [Fig. 3]. Empirically
we find, that the relevant scattering cross section
$\Sigma_{*}\approx\sigma^{2}$ is smaller than the geometric cross section
$\pi\sigma^{2}$, leading to a mean-free path length of $l_{*}\approx
1/n\sigma^{2}$. The smaller scattering cross section can be readily
interpreted since only head-on collisions effectively stop the directed motion
[in the lattice variant no empirical correction was necessary, the scattering
cross section coincides with the geometrical one Leitmann and Franosch
(2017)]. Matching the superdiffusion with the short-time asymptote $2D_{a}t$
yields as crossover time $\sim\text{Pe}^{-3/2}$. Similarly the terminal time
of superdiffusion is set by the time the probe particle needs to pass a bath
particle $\sigma/\mu F\sim 1/\text{Pe}$ leading by crossover matching to the
scaling prediction $D_{z}(\text{Pe})~{}\sim n\,\text{Pe}$ for the diffusion
coefficient consistent with Ref. Zia and Brady (2010). Thus, the window of
superdiffusion expands as $\text{Pe}^{1/2}$ as the force is increased.

Figure 3:  Time-dependent fluctuations $\text{Var}_{z}(t)=\langle\Delta
z_{a}(t)^{2}\rangle-\langle\Delta z_{a}(t)\rangle^{2}$ and local exponent
$\alpha(t)=\mathrm{d}\\!\ln[\text{Var}_{z}(t)]/\mathrm{d}\\!\ln(t)$ (inset) of
the probe particle along the applied force obtained from Brownian dynamcis
simulation of equal sized colloids ($D_{b}=D_{a}$) for different strength of
the driving. The solid lines correspond to the diffusive asymptote
$\text{Var}_{z}(t)=2D_{z}(\text{Pe})t$ and the dashed lines are the asymptotic
model $\text{Var}_{z}(t)=(D_{a}/D_{r})^{2}(\mu Ft)^{3}/3l_{*}$ with the
emperical mean-free path length $l_{*}\approx 1/n\sigma^{2}$.

###  III.3 Long-time behavior

Let us specialize our time-dependent solution to the stationary state. The
stationary mobility $\mu(\text{Pe})$ follows as special case from Eq. (47) via

$\displaystyle\mu(\text{Pe})/\mu=\lim_{\omega\to
0}(-\mathsf{i}\omega)\hat{\mu}(\omega,\text{Pe})/\mu=1+2\pi
n\sigma^{3}\frac{D_{a}}{D_{r}}\Delta\mu_{0}(\text{Pe}),$ (66)

where the coefficient $\Delta\mu_{0}(\text{Pe})$ is defined via a low-
frequency expansions of the auxiliary function

$\displaystyle\Delta\hat{\mu}(\omega,\text{Pe})=\Delta\mu_{0}(\text{Pe})+(-\mathsf{i}\omega\tau)\Delta\mu_{1}(\text{Pe})+\dotsb,\quad\omega\to
0.$ (67)

For small driving and small frequency, we can determine the coefficients
$a_{\ell}$ of the auxiliary function $\Delta\hat{\mu}(\omega,\text{Pe})$ [Eq.
(56)] by an inversion of the tridiagonal matrix equation [Eq. (32)] and by
considering the frequency-independent contribution:

$\displaystyle\Delta\mu_{0}(\text{Pe})=-\frac{1}{3}+\frac{2}{45}\text{Pe}^{2}-\frac{1}{24}|\text{Pe}|^{3}+\frac{128}{4725}\text{Pe}^{4}+\mathcal{O}(|\text{Pe}|^{5}),\quad\text{Pe}\to
0.$ (68)

Inserting this result into Eq. (66), we recover the previously derived
asymptotic expansion [Eq. (1)] Squires and Brady (2005). In particular, in
equilibirium, $\Delta\mu_{0}(\text{Pe}\to 0)=-1/3$ and for equal-sized
colloids with relative diffusion coefficient $D_{r}=2D_{a}$ and packing
fraction $\varphi=\pi n\sigma^{3}/6$, we obtain the known result from
equilibrium, $\mu(\text{Pe}=0)/\mu=1-2\varphi$ Hanna  _et al._ (1982). In the
limit of large forces, $\text{Pe}\to\infty$, the analytic solution for the
mobility [Eq. (66)] approaches the limit
$\mu(\text{Pe}\to\infty)/\mu=1-(D_{a}/D_{r})\pi n\sigma^{3}/3$ [Fig. 4],
derived earlier in terms of a boundary layer analysis Squires and Brady
(2005).

Figure 4:  Density-induced suppression of the stationary mobility
$\mu(\text{Pe})$. Dashed lines indicate the asymptotic expansion for
increasing order ($2$, $4$, $6$, $10$, and $20$). Lines correspond to the
theory and symbols represent Brownian-dynamics simulations of equal sized
colloids, $D_{b}=D_{a}$, and Lorentz systems with $D_{b}=0$.

Similarly, for the stationary diffusion coefficient $D_{z}(\text{Pe})$, we
perform the limit of small frequencies in Eq. (62) and obtain

$\displaystyle D_{z}(\text{Pe})/D_{a}=\lim_{\omega\to
0}(-\mathsf{i}\omega)\hat{D}_{z}(\omega,\text{Pe})/D_{a}=1+2\pi
n\sigma^{3}\frac{D_{a}}{D_{r}}[2\Delta
R_{0}(\text{Pe})+\text{Pe}^{2}\Delta\mu_{1}(\text{Pe})],$ (69)

where the coefficient $\Delta\mu_{1}(\text{Pe})$ is defined in Eq. (67) and
the coefficient $\Delta R_{0}(\text{Pe}):=\Delta\hat{R}(\omega\to
0,\text{Pe})$ is the small-frequeny limit of the auxiliary function
$\Delta\hat{R}(\omega,\text{Pe})$ [Eq. (57)]. Similar to the coefficient
$\Delta\mu_{0}(\text{Pe})$, we determine the coefficients $a_{\ell}$ and
$\partial a_{\ell}/\partial q$ by matrix inversion of the tridiagonal matrix
equation [Eq. (32)] in powers of the the Péclet number and the frequency. The
coefficients are then obtained as

$\displaystyle\Delta R_{0}(\text{Pe})$
$\displaystyle=-\frac{1}{6}+\frac{17}{270}\text{Pe}^{2}-\frac{11}{144}|\text{Pe}|^{3}+\frac{1679}{28350}\text{Pe}^{4}+\mathcal{O}(|\text{Pe}|)^{5},\quad\text{Pe}\to
0,$ (70) $\displaystyle\Delta\mu_{1}(\text{Pe})$
$\displaystyle=\frac{1}{6}-\frac{1}{6}|\text{Pe}|+\frac{101}{810}\text{Pe}^{2}-\frac{17}{216}|\text{Pe}|^{3}+\frac{310571}{6804000}\text{Pe}^{4}+\mathcal{O}(|\text{Pe}|)^{5},\quad\text{Pe}\to
0,$ (71)

and the series expansion of the diffusion coefficient parallel to the applied
field is calculated to

$\displaystyle\begin{split}D_{z}(\text{Pe})/D_{a}=1-\frac{2\pi
n\sigma^{3}}{3}\frac{D_{a}}{D_{r}}\biggl{[}1&-\frac{79}{90}\text{Pe}^{2}+\frac{23}{24}|\text{Pe}|^{3}-\frac{6893}{9450}\text{Pe}^{4}+\mathcal{O}(|\text{Pe}|)^{5}\biggr{]}.\end{split}$
(72)

With our solution, we recover the leading correction
$\mathcal{O}(\text{Pe}^{2})$ to the equilibrium case, which has been
calculated earlier Zia and Brady (2010). Furthermore, our calculation extends
this result to arbitrary order and reveals the emergence of nonanalytic
contributions similar to the mobility in the stationary state.

As can be inferred from the stationary mobility [Fig. 4] and the stationary
diffusion coefficient [Fig. 5], both series expansions [Eqs. (1) and (72)]
break down already at moderate driving. However, the numerical evaluation via
the matrix inversion [Eq. (32)] is valid for the full range of Péclet numbers
which has been shown first by Khair and Brady Khair and Brady (2006) for the
mobility. The limiting value of the stationary mobility as well as the
asymptotic behavior $\mathcal{O}(\text{Pe})$ of the stationary diffusion
coefficient for large Péclet numbers [Fig. 5] Zia and Brady (2010) is nicely
corroborated in the Brownian dynamics simulations.

Figure 5:  Stationary-state force-induced diffusion coefficient
$D_{z}^{\text{ind}}:=D_{z}(\text{Pe})-D_{z}(\text{Pe}\to 0)$. Dashed lines
represent the asymptotic behavior for small (increasing order $2$, $3$, $5$,
$9$, $31$) and large driving. Lines correspond to the theory and symbols
represent Brownian-dynamics simulations of equal sized colloids,
$D_{b}=D_{a}$, and Lorentz systems $D_{b}=0$.

###  III.4 Time-dependent pair-distribution function

Our solution provides also the time-dependent pair-distribution function
$g(\vec{r},t)$ for the relative distance $\vec{r}$ of the tracer and the bath
particles by integrating the conditional probability distribution
$\psi(\vec{r},t|\vec{r}^{\prime})$ over all $\vec{r}^{\prime}$. In the
frequency domain,
$\hat{g}(\vec{r},\omega)=\hat{\psi}_{q=0}(\vec{r},\omega)\sqrt{V}$ and with
Eq. (28), we find for $|\vec{r}|>\sigma$

$-\mathsf{i}\omega\hat{g}(\vec{r},\omega)=1+e^{\mu\vec{F}\cdot\vec{r}/2D_{r}}\sum_{\ell=0}^{\infty}a_{\ell}\frac{\text{k}_{\ell}(\kappa
r)}{\text{k}_{\ell}(\kappa\sigma)}\text{P}_{\ell}(\cos\vartheta),$ (73)

where the coefficients $a_{\ell}\equiv
a_{\ell}(q\sigma,\kappa\sigma,\text{Pe})$ are to be evaluated at zero
wavenumber $q=0$. In particular, in equilibrium all coefficients vanish
$a_{\ell}=0$ and the pair-distribution function reduces to the stationary step
function.

Simulation results for the time-evolution of the pair-distribution function
are displayed in Fig. 6. Shortly after switching on the force, the pair-
distribution function is still almost spherically symmetric, there has been no
time to propagate the information of the strong perturbing force to the
surroundings. Quickly the tracer’s motion is obstructed by the bath particles
such that probability accumulates in front of the tracer particle. For longer
times, a trailing wake evolves where probability is depleted behind the tracer
particle. The numerical results for the pair-distribution function in the
frequency domain shown in Fig. 6 corroborate this picture. For small
frequencies we recover the stationary distribution calculated in Refs. Squires
and Brady (2005); Khair and Brady (2006). For high Péclet numbers the
stationary two-particle distribution function becomes strongly asymmetric as
has been shown first Ref. Squires and Brady (2005) and more systematically in
Ref. Khair and Brady (2006). Probability piles up in a narrow boundary layer
in front of the pulled particle particle and leaves behind a wake extending to
larger and larger distances ${\cal O}(\text{Pe})$ as the driving force is
increased.

Figure 6:  Top: Simulation results for the time-dependent pair-distribution
function $g(\vec{r},t)$ for Péclet number $\text{Pe}=32$ and density
$n\sigma^{3}=0.01$. As time progresses $t\cdot D_{a}/\sigma^{2}=4\cdot
10^{-4},4\cdot 10^{-2},4$ (left to right) the tracer piles up probability in
front and leaves a trail of depleted probability. Bottom: Analytic result for
the real part of the pair-distribution function in the frequency domain
$\operatorname{Re}[-\mathsf{i}\omega\hat{g}(\vec{r},\omega)]$ at the same
driving $\text{Pe}=32$ for frequencies $\omega\cdot\sigma^{2}/2\pi
D_{a}=1/(4\cdot 10^{-4}),1/(4\cdot 10^{-2}),1/4$ (left to right).

##  IV Summary and conclusion

We have derived an analytic solution for the full time-dependent response of a
probe particle driven out of equilibrium by a step force in first order of the
density of bath particles (discarding inertial and hydrodynamic effects). The
response is completely encoded in the self-energy from which, in principle,
all moments along the force can be generated. In comparison to the known self-
energy in equilibrium, the force results in a complex shift of the frequency
$-\mathsf{i}\omega\tau\mapsto-\mathsf{i}\omega\tau+(\text{Pe}/2)^{2}=\kappa^{2}\sigma^{2}$
encoded in the complex wavenumber $\kappa\sigma$, an explicit regular
variation of the expansion coefficients $a_{\ell}$ via the boundary condition,
and a shift in the wavenumber $q\sigma\mapsto
q\sigma+\mathsf{i}\,\text{Pe}/2$. The shift in the frequency reveals that the
nonanalytic frequency behavior in equilibrium, as manifested in the long-time
tails, and the nonanalytic dependence on the driving, are merely two sides of
the same coin. In particular, this explains why at any finite driving the
approach of the nonlinear mobility to its stationary value becomes
exponentially fast. Furthermore, this reveals that for finite times, all
response functions are analytic functions in the driving, however, this does
not hold for infinite times since the limits do not commute.

The emergence of non-analytic behavior and a divergent time scale for
$\text{Pe}\to 0$ calls for an explanation in terms of physics. The stationary
Smoluchowski equation is non-uniform at small Péclet number [the same physics
has been discussed for the advection-diffusion equation in the seminal
contribution by Acrivos and Taylor Acrivos and Taylor (1962)]: it displays an
outer region at distances $r\gg\sigma/\text{Pe}$ where advection dominates,
and an inner one where diffusion is the dominant contribution. For the
transport coefficient already to order ${\cal O}(\text{Pe}^{2})$ both regions
need to be calculated and matching the solutions makes the non-analytic
contributions evident Khair and Brady (2006). For non-zero frequencies or
finite times the advected Smoluchowski equation is regular and solutions decay
exponentially fast on the ’Skin penetration’ depth
$\sqrt{D_{r}/\omega}\sim\sqrt{D_{r}t}$. Correspondingly, if the Skin
penetration depth is smaller than the inner region, i.e. for times
$\sqrt{D_{r}t},\lesssim\sigma/\text{Pe}$, i.e.
$t\lesssim\tau_{F}=\tau/\text{Pe}^{2}$, one can safely ignore the presence of
the outer region and the system behaves as in linear response. In contrast,
for large enough times, the Skin penetration depth covers the outer region and
the nonlinearities becomes important. The interplay of the divergent boundary
layer and the Skin penetration depth is the origin of non-commuting limits.
Mathematically related, but not quite identical, is the problem of the time-
dependent motion at small but finite Reynolds number Lovalenti and Brady
(1993).

The full solution provides the first direct access to the time-dependent
fluctuations along the force. Here we explicitly characterized the transient
superdiffusion which connects the short-time bare diffusion and the long-time
enhanced diffusion Zia and Brady (2010). The emergence of superdiffusion has
been rationalized by considering the distribution of the free path lengths.

Let us comment also on hydrodynamic interactions. Progress has been made for
certain limiting cases: In equilibrium, the stationary diffusion coefficients
have been estimated by including instantaneous hydrodynamic interactions at
the level of the Oseen tensor, as well as including near-field corrections.
Results for the density-induced suppression depend somewhat on the
approximation of the hydrodynamic interactions but are close to the diffusion
coefficient neglecting hydrodynamics Hanna  _et al._ (1982); Batchelor (1976);
Rallison and Hinch (1986); Hoh and Zia (2016). For the driven case,
hydrodynamic interactions have been accounted for in the stationary state both
for the structure deformation as well as the mobility for all Péclet numbers
Khair and Brady (2005). The force-induced corrections to the stationary
diffusion coefficient have been elaborated only recently, and it has been
shown that hydrodynamic interactions do not lead to qualitatively new
behavior, although the numerical values change. In particular, the low-force
corrections are still $\mathcal{O}(\text{Pe}^{2})$, and in the regime of
strong forces they still scale as $\mathcal{O}(\text{Pe})$, however with a
slow convergence to the asymptotic result Hoh and Zia (2016).

A second effect due to hydrodynamics arises due to frequency-dependent
hydrodynamic interactions implying hydrodynamic memory. Here, the slow vortex
diffusion of transverse momentum in the fluid leads to a characteristic
algebraic decay of the form $\simeq Bt^{-3/2}$, $B>0$ for the velocity-
autocorrelation function of the particleAlder and Wainwright (1970); Zwanzig
and Bixon (1970). Hence, one may ask the question if the long-time tail due to
hydrodynamics dominates the long-time tail $\simeq-At^{-5/2}$, $A>0$ due to
repeated collisions of the probe particle with the bath particles. Taking
physical values from experiments Franosch  _et al._ (2011); Jeney  _et al._
(2008) and comparing both tails shows that, in principle, a window of time
opens where the algebraic decay due to the collision of the probe particle
with bath particles dominates before hydrodynamics becomes relevant at larger
times.

Our predictions for the time-dependent response of a driven colloid can be
tested, in principle, in laboratory experiments on a colloidal suspension via
particle tracking. The general scenario persists also for soft spheres and is
not restricted to dilute systems. To first order in the density, the different
diffusivities of probe and bath particles can be trivially accounted for,
e.g., by a rescaling of time. In particular, the case of a dilute and quenched
array of obstacles (Lorentz model) is also included. Simulation of the Lorentz
model, where the bath particles are pinned, were added to the comparison of
the mobility and diffusivity with theory in Figs. 4 and 5.

Our analysis of the driven colloid shows that the nonequilibrium stationary
state is inherently a nonanalytic function of the driving such that the
transport coefficients can not be expanded in a Taylor series beyond linear
response. Although this has been derived to first order in the densities only,
this behavior is anticipated to be generic and valid for arbitrary densities.
Arguably, such relations should hold universally in a general nonlinear
response framework. This view is supported by recent predictions for a two-
dimensional or three-dimensional driven lattice Lorentz gas Leitmann and
Franosch (2013, 2017); Leitmann  _et al._ (2018), where qualitatively the same
scenario applies.

###### Acknowledgements.

We gratefully acknowledge support by the DFG research unit FOR1394 “Nonlinear
response to probe vitrification” and MINECO and ERDF under project
FIS2015-69022P (AMP). The computational results presented have been achieved
(in part) using the HPC infrastructure LEO of the University of Innsbruck.

##  Appendix A Computer Simulations

We have simulated the motion of the pulled probe particle in the presence of
bath particles interacting via a hard-core potential using event-driven pseudo
Brownian dynamics simulation Scala  _et al._ (2007) ignoring inertial effects
or hydrodynamic interactions. The starting point for the stochastic simulation
of the suspension is the Langevin equation

$\displaystyle\mathrm{d}\vec{r}=\sqrt{2D_{a}}\vec{\eta}(t)\mathrm{d}t+\mu
F\vec{e}_{z}\mathrm{d}t,$ (74)

which describes the change of position $\mathrm{d}\vec{r}$ of the probe
particle in terms of the Gaussian white noise process $\vec{\eta}$ with zero
mean and covariance
$\langle\eta_{i}(t)\eta_{j}(t)\rangle=\delta_{ij}\delta(t-t^{\prime})$. For
the bath particles, the same equation with $F=0$ holds. The Langevin equation
[Eq. (74)] is implemented by introducing a fixed Brownian time step
$\tau_{B}$, such that for every step, the pseudo-velocity

$\displaystyle\vec{v}=\sqrt{\frac{2D_{a}}{\tau_{B}}}\mathcal{N}_{\eta}+\mu
F\vec{e}_{z}$ (75)

is assigned to the probe and the bath particles ($F=0$). Between these
Brownian interrupts the particles move with constant velocities and collide
elastically Scala  _et al._ (2007). The normal distributed random variable
$\mathcal{N}_{\eta}$ arises from discretization of the white noise and has
zero mean and unit variance. The Brownian time step $\tau_{B}$ should be much
smaller than the diffusion time $\sigma^{2}/D_{a}$ and the drift time
$\sigma/\mu F$.

For the colloidal case, we have used equilibrated configurations consisting of
$1000$ particles at a fixed number density $n\sigma^{3}=0.01$ and a typical
value of the Brownian time step is $\tau_{B}\simeq
10^{-3}(\sigma^{2}/D_{a})/\text{max}(1,\text{Pe})$. For each data set, we have
simulated at least $10^{7}$ independent trajectories. For the case of the
Lorentz system, we freeze the dynamics of the bath particles and apply the
stochastic dynamics only to the tracer particle.

## References

  * Squires and Brady (2005) Todd M. Squires and John F. Brady, “A simple paradigm for active and nonlinear microrheology,” Phys. Fluids 17, 073101 (2005).
  * Khair and Brady (2006) Aditya S. Khair and John F. Brady, “Single particle motion in colloidal dispersions: a simple model for active and nonlinear microrheology,” J. Fluid Mech. 557, 73–117 (2006).
  * Wilson  _et al._ (2009) L. G. Wilson, A. W. Harrison, A. B. Schofield, J. Arlt, and W. C. K. Poon, “Passive and active microrheology of hard-sphere colloids,” J. Phys. Chem. B 113, 3806–3812 (2009).
  * Zia and Brady (2010) Roseanna N. Zia and John F. Brady, “Single-particle motion in colloids: force-induced diffusion,” J. Fluid Mech. 658, 188–210 (2010).
  * Sriram  _et al._ (2010) Indira Sriram, Alexander Meyer, and Eric M. Furst, “Active microrheology of a colloidal suspension in the direct collision limit,” Phys. Fluids 22, 062003 (2010).
  * Wilson and Poon (2011) Laurence G. Wilson and Wilson C. K. Poon, “Small-world rheology: an introduction to probe-based active microrheology,” Phys. Chem. Chem. Phys. 13, 10617–10630 (2011).
  * Zia and Brady (2012) Roseanna N. Zia and John. F. Brady, “Microviscosity, microdiffusivity, and normal stresses in colloidal dispersions,” Journal of Rheology 56, 1175–1208 (2012).
  * Swan and Zia (2013) James W. Swan and Roseanna N. Zia, “Active microrheology: Fixed-velocity versus fixed-force,” Phys. Fluids 25, 083303 (2013).
  * Puertas and Voigtmann (2014) Antonio M. Puertas and Thomas Voigtmann, “Microrheology of colloidal systems,” J. Phys. Condens. Matter 26, 243101 (2014).
  * Zia and Brady (2015) Roseanna N. Zia and John F. Brady, “Theoretical microrheology,” in  _Complex Fluids in Biological Systems: Experiment, Theory, and Computation_, edited by E. Saverio Spagnolie (Springer New York, New York, NY, 2015) pp. 113–157.
  * Furst and Squires (2017) Eric M. Furst and Todd M. Squires,  _Microrheology_ (Oxford University Press, Oxford, 2017).
  * Winter  _et al._ (2012) D. Winter, J. Horbach, P. Virnau, and K. Binder, “Active nonlinear microrheology in a glass-forming Yukawa fluid,” Phys. Rev. Lett. 108, 028303 (2012).
  * Winter and Horbach (2013) D. Winter and J. Horbach, “Nonlinear active micro-rheology in a glass-forming soft-sphere mixture,” J. Chem. Phys. 138, 12A512 (2013).
  * Gazuz  _et al._ (2009) I. Gazuz, A. M. Puertas, Th. Voigtmann, and M. Fuchs, “Active and nonlinear microrheology in dense colloidal suspensions,” Phys. Rev. Lett. 102, 248302 (2009).
  * Gnann  _et al._ (2011) M. V. Gnann, I. Gazuz, A. M. Puertas, M. Fuchs, and Th. Voigtmann, “Schematic models for active nonlinear microrheology,” Soft Matter 7, 1390–1396 (2011).
  * Gnann and Voigtmann (2012) M. V. Gnann and Th. Voigtmann, “Asymptotic analysis of mode-coupling theory of active nonlinear microrheology,” Phys. Rev. E 86, 011406 (2012).
  * Harrer  _et al._ (2012) Ch J Harrer, D Winter, J Horbach, M Fuchs, and Th Voigtmann, “Force-induced diffusion in microrheology,” J. Phys. Condens. Matter 24, 464105 (2012).
  * Gazuz and Fuchs (2013) I. Gazuz and M. Fuchs, “Nonlinear microrheology of dense colloidal suspensions: A mode-coupling theory,” Phys. Rev. E 87, 032304 (2013).
  * Wang  _et al._ (2014) Ting Wang, Matthias Grob, Annette Zippelius, and Matthias Sperl, “Active microrheology of driven granular particles,” Phys. Rev. E 89, 042209 (2014).
  * Gruber  _et al._ (2016) M. Gruber, G. C. Abade, A. M. Puertas, and M. Fuchs, “Active microrheology in a colloidal glass,” Phys. Rev. E 94, 042602 (2016).
  * Jack  _et al._ (2008) Robert L. Jack, David Kelsey, Juan P. Garrahan, and David Chandler, “Negative differential mobility of weakly driven particles in models of glass formers,” Phys. Rev. E 78, 011506 (2008).
  * Schroer and Heuer (2013a) Carsten F. E. Schroer and Andreas Heuer, “Anomalous diffusion of driven particles in supercooled liquids,” Phys. Rev. Lett. 110, 067801 (2013a).
  * Schroer and Heuer (2013b) Carsten F. E. Schroer and Andreas Heuer, “Microrheology of supercooled liquids in terms of a continuous time random walk,” J. Chem. Phys. 138, 12A518 (2013b).
  * Burioni  _et al._ (2014) R. Burioni, G. Gradenigo, A. Sarracino, A. Vezzani, and A. Vulpiani, “Scaling properties of field-induced superdiffusion in continuous time random walks,” Commun. Theor. Phys. 62, 514 (2014).
  * Démery  _et al._ (2014) Vincent Démery, Olivier Bénichou, and Hugo Jacquin, “Generalized Langevin equations for a driven tracer in dense soft colloids: construction and applications,” New J. Phys. 16, 053032 (2014).
  * Démery (2015) Vincent Démery, “Mean-field microrheology of a very soft colloidal suspension: Inertia induces shear thickening,” Phys. Rev. E 91, 062301 (2015).
  * Wang and Sperl (2016) Ting Wang and Matthias Sperl, “Thinning and thickening in active microrheology,” Phys. Rev. E 93, 022606 (2016).
  * Bénichou  _et al._ (2013a) Olivier Bénichou, Anna Bodrova, Dipanjan Chakraborty, Pierre Illien, Adam Law, Carlos Mejía-Monasterio, Gleb Oshanin, and Raphaël Voituriez, “Geometry-induced superdiffusion in driven crowded systems,” Phys. Rev. Lett. 111, 260601 (2013a).
  * Illien  _et al._ (2013) P. Illien, O. Bénichou, C. Mejía-Monasterio, G. Oshanin, and R. Voituriez, “Active transport in dense diffusive single-file systems,” Phys. Rev. Lett. 111, 038102 (2013).
  * Bénichou  _et al._ (2013b) Olivier Bénichou, Pierre Illien, Carlos Mejía-Monasterio, and Gleb Oshanin, “A biased intruder in a dense quiescent medium: looking beyond the force–velocity relation,” Journal of Statistical Mechanics: Theory and Experiment 2013, P05008 (2013b).
  * Bénichou  _et al._ (2014) O. Bénichou, P. Illien, G. Oshanin, A. Sarracino, and R. Voituriez, “Microscopic theory for negative differential mobility in crowded environments,” Phys. Rev. Lett. 113, 268002 (2014).
  * Illien  _et al._ (2014) Pierre Illien, Olivier Bénichou, Gleb Oshanin, and Raphaël Voituriez, “Velocity anomaly of a driven tracer in a confined crowded environment,” Phys. Rev. Lett. 113, 030603 (2014).
  * (33) Pierre Illien, Olivier Bénichou, Gleb Oshanin, and Raphaël Voituriez, “Distribution of the position of a driven tracer in a hardcore lattice gas,” J. Stat. Mech. (2015), P11016.
  * Leitmann and Franosch (2013) Sebastian Leitmann and Thomas Franosch, “Nonlinear response in the driven lattice Lorentz gas,” Phys. Rev. Lett. 111, 190603 (2013).
  * Leitmann and Franosch (2017) Sebastian Leitmann and Thomas Franosch, “Time-dependent fluctuations and superdiffusivity in the driven lattice Lorentz gas,” Phys. Rev. Lett. 118, 018001 (2017).
  * Leitmann  _et al._ (2018) Sebastian Leitmann, Olivier Bénichou, and Thomas Franosch, “Time-dependent dynamics of the three-dimensional driven lattice Lorentz gas,” J. Phys A: Math. and Theoretical (2018), in print.
  * Hanna  _et al._ (1982) Sarwat Hanna, Walter Hess, and Rudolf Klein, “Self-diffusion of spherical Brownian particles with hard-core interaction,” Physica A 111, 181 (1982).
  * Hanna  _et al._ (1981) S Hanna, W Hess, and R Klein, “The velocity autocorrelation function of an overdamped Brownian system with hard-core intraction,” J. Phys. A. 14, L493 (1981).
  * Ackerson and Fleishman (1982) Bruce J. Ackerson and L. Fleishman, “Correlations for dilute hard core suspensions,” J. Chem. Phys. 76, 2675–2679 (1982).
  * Felderhof and Jones (1983a) B.U. Felderhof and R.B. Jones, “Diffusion in hard sphere suspensions,” Physica A 122, 89 – 104 (1983a).
  * Khair and Brady (2005) Aditya S. Khair and John F. Brady, “Microviscoelasticity of colloidal dispersions,” Journal of Rheology 49, 1449 (2005).
  * Ernst and Weyland (1971) M. H. Ernst and A. Weyland, “Long time behaviour of the velocity auto-correlation function in a Lorentz gas,” Phys. Lett. A 34, 39 – 40 (1971).
  * van Beijeren (1982) Henk van Beijeren, “Transport properties of stochastic Lorentz models,” Rev. Mod. Phys. 54, 195–234 (1982).
  * Nieuwenhuizen  _et al._ (1986) Th. M. Nieuwenhuizen, P. F. J. van Velthoven, and M. H. Ernst, “Diffusion and long-time tails in a two-dimensional site-percolation model,” Phys. Rev. Lett. 57, 2477–2480 (1986).
  * Höfling and Franosch (2007) Felix Höfling and Thomas Franosch, “Crossover in the slow decay of dynamic correlations in the Lorentz model,” Phys. Rev. Lett. 98, 140601 (2007).
  * Franosch  _et al._ (2010) Thomas Franosch, Felix Höfling, Teresa Bauer, and Erwin Frey, “Persistent memory for a Brownian walker in a random array of obstacles,” Chem. Phys. 375, 540 – 547 (2010).
  * Bauer  _et al._ (2010) T. Bauer, F. Höfling, T. Munk, E. Frey, and T. Franosch, “The localization transition of the two-dimensional Lorentz model,” The European Physical Journal Special Topics 189, 103–118 (2010).
  * Felderhof and Jones (1983b) B.U. Felderhof and R.B. Jones, “Cluster expansion of the diffusion kernel of a suspension of interacting Brownian particles,” Physica A 121, 329 – 344 (1983b).
  * Nelson and Shnerb (1998) David R. Nelson and Nadav M. Shnerb, “Non-hermitian localization and population biology,” Phys. Rev. E 58, 1383–1403 (1998).
  * Olver  _et al._ (2010) F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, eds.,  _NIST Handbook of Mathematical Functions_ (Cambridge University Press, New York, NY, 2010) print companion to DLMF .
  * (51) DLMF, “NIST digital library of mathematical functions,” http://dlmf.nist.gov/, Release 1.0.15 of 2017-06-01, online companion to Olver  _et al._ (2010).
  * Acrivos and Taylor (1962) Andreas Acrivos and Thomas D. Taylor, “Heat and mass transfer from single spheres in Stokes flow,” The Physics of Fluids 5, 387 (1962).
  * Lovalenti and Brady (1993) Phillip M. Lovalenti and John F. Brady, “The hydrodynamic force on a rigid particle undergoing arbitrary time-dependent motion at small Reynolds number,” J. Fluid Mech. 256, 561 (1993).
  * Batchelor (1976) G. K. Batchelor, “Brownian diffusion of particles with hydrodynamic interaction,” Journal of Fluid Mechanics 74, 1–29 (1976).
  * Rallison and Hinch (1986) J. M. Rallison and E. J. Hinch, “The effect of particle interactions on dynamic light scattering from a dilute suspension,” Journal of Fluid Mechanics 167, 131–168 (1986).
  * Hoh and Zia (2016) N. J. Hoh and R. N. Zia, “Force-induced diffusion in suspensions of hydrodynamically interacting colloids,” J. Fluid Mech. 795, 739–783 (2016).
  * Alder and Wainwright (1970) B. J. Alder and T. E. Wainwright, “Decay of the velocity autocorrelation function,” Phys. Rev. A 1, 18–21 (1970).
  * Zwanzig and Bixon (1970) Robert Zwanzig and Mordechai Bixon, “Hydrodynamic theory of the velocity correlation function,” Phys. Rev. A 2, 2005–2012 (1970).
  * Franosch  _et al._ (2011) Thomas Franosch, Matthias Grimm, Maxim Belushkin, Flavio M. Mor, Giuseppe Foffi, László Forró, and Sylvia Jeney, “Resonances arising from hydrodynamic memory in Brownian motion,” Nature (London) 478, 85 (2011).
  * Jeney  _et al._ (2008) Sylvia Jeney, Branimir Lukić, Jonas A. Kraus, Thomas Franosch, and László Forró, “Anisotropic memory effects in confined colloidal diffusion,” Phys. Rev. Lett. 100, 240604 (2008).
  * Scala  _et al._ (2007) A. Scala, Th. Voigtmann, and C. De Michele, “Event-driven brownian dynamics for hard spheres,” J. Chem. Phys. 126, 134109 (2007).