# Dynamic Density Functional Theory with Inertia and Background Flow

R. D. Mills-Williams  Edinburgh Designs Ltd, 27 Ratcliffe Terrace, Edinburgh,
United Kingdom, EH9 1SX, United Kingdom  B. D. Goddard  School of Mathematics
and Maxwell Institute for Mathematical Sciences, University of Edinburgh, EH9
3FD, United Kingdom  A. J. Archer  Department of Mathematical Sciences,
Loughborough University, Loughborough, Leicestershire LE11 3TU, United Kingdom

###### Abstract

We present dynamic density functional theory (DDFT) incorporating general
inhomogeneous, incompressible, time dependent background flows and inertia,
describing externally driven passive colloidal systems out of equilibrium. We
start by considering the underlying nonequilibrium Langevin dynamics,
including the effect of the local velocity of the surrounding liquid bath, to
obtain the nonlinear, nonlocal partial differential equations governing the
evolution of the (coarse–grained) density and velocity fields describing the
dynamics of colloids. Additionally, we show both with heuristic arguments, and
by numerical solution, that our equations and solutions agree with existing
DDFTs in the overdamped (high friction) limit. We provide numerical solutions
that model the flow of hard spheres, in both unbounded and confined domains,
and compare to previously–derived DDFTs with and without the background flow.

††preprint: AIP/123-QED

##  I Introduction

Statistical mechanical methods such as dynamic density functional theories
(DDFTs) are increasingly popular approaches to model out of equilibrium
colloidal systems. Having been applied to a multitude of systems, DDFTs have
been shown to accurately describe many nonequilibrium fluid properties, such
as: phase separation in binary fluidsArcher (2005), inertia, withGoddard  _et
al._ (2012) and withoutArcher (2009) hydrodynamic interactions (HI),
orientationWittkowski and Löwen (2011); Durán-Olivencia, Goddard, and
Kalliadasis (2016), and thermal fluctuations on the hydrodynamic scaleDurán-
Olivencia  _et al._ (2017). Such methods predict the evolution of the one-body
density distribution out of equilibrium with remarkable accuracy when compared
to the underlying stochastic (Langevin) equationsHansen and McDonald (2013);
te Vrugt, Löwen, and Wittkowski (2020). Most existing DDFTs, however, are
derived from thermostated equilibrium molecular dynamics equations, such as
Langevin’s equations, in a _static_ heat bath, and therefore do not allow for
the bath in which the colloidal particles are immersed to be flowing. This
limitation for existing models of nonequilibrium colloidal systems is due to a
central assumption in deriving the Langevin equations: that the atoms or
molecules which make up the heat bath undergo Hamiltonian ballistic dynamics,
exchanging momentum with the larger colloidal particles via elastic
collisions. This is correct in so far as the resulting dynamics satisfies the
fluctuation–dissipation relation Callen and Welton (1951) for static baths.
But for a general nonequilibrium system with a given nonzero background flow,
a Hamiltonian description of the bath will not produce the correct dynamics or
fluctuation–dissipation relation. This has been shown rigorously by Dobson
_et al._ (2013) and is due to that fact the finite local strain rate on the
colloidal particles between bath atom/molecule collisions is neglected under a
Hamiltonian description of the bath dynamics.

Elsewhere, classical Nonequilibrium Molecular Dynamics (NEMD) has been applied
to colloidal dispersions driven via inhomogeneous boundary-driven flows, for
example Newton’s equations combined with Lees–Edwards boundary conditions
(BCs) in the case of shear flow Evans and Morriss (2007) or Kraynik–Reinelt
BCs in the case of elongation flow Todd and Daivis (1999). These methods
however are limited to periodic problems, making more general NEMD methods
increasingly desirable. Several strategies have been proposed to sample
molecular systems with general, inhomogeneous background flows, including the
SLLOD and g-SLLOD equations of motion for cold systemsTuckerman  _et al._
(1997); Todd and Daivis (2007) with proposals for thermostats (Nosé–Hoover,
Gaussian isokinetic) for temperature control at high strain rate. However,
such thermostats are problematic for, e.g., high strain rate planar-shear-flow
with S-shaped streaming velocity profiles. Being kinetic thermostats, they
interpret deviations away from linearity as excess thermal energy leading to
an erroneous ordering effect, whereby for a simple atomic fluid, the system
displays an accentuated alignment to the velocity streamlines Todd and Daivis
(2007). Additionally, Nosé–Hoover dynamics can be non-ergodicLegoll, Luskin,
and Moeckel (2007).

In the deterministic DDFT continuum theory setting, there have been several
previous approaches to study driven particle systems. Some of the earliest of
these considered the overdamped dynamics of the solvent particles in a polymer
solution, and their hydrodynamic interactions. Penna, Dzubiella, and Tarazona
(2003) determined a DDFT accounting for the drifting effect of a single
colloid particle in the case of both a non-interacting and an interacting
polymer. Dzubiella, Löwen, and Likos (2003), and separately Krüger and
Rauscher (2007), considered the case of two colloids, in order to calculate
the effective depletion force of the polymer solute, the latter considering
the deterministic Smoluchowski equation for the Brownian solute particles and
solving it in bispherical coordinates. As has been previously observed, such
models only modify the density of the solute and do not account for
perturbations of the flow itself due to the solvent interacting with the
colloid, essentially meaning the velocity field flows through the particle
boundary. To go beyond this, starting from the underlying Brownian dynamics,
Rauscher  _et al._ (2007) derived a DDFT for advected colloidal particles
externally driven by an inhomogeneous field field, based on a classical
solution for incompressible low Reynolds number flow past a sphere. Extensions
to this DDFT including HI with walls and between particles was again treated
by Rauscher (2010) (including the formal details of the moment closure), and
there have been many practical applications including: dynamics in confined
geometries Almenar and Rauscher (2011), advected colloids through DNA Gutsche
_et al._ (2008), phase transitions in microrheology Reinhardt, Scacchi, and
Brader (2014), as well as the investigation of the hard sphere interactions in
the presence of flow induced hydrodynamic lift in confined geometries Yu  _et
al._ (2017). The influence of solvent flow is also thought to be important in
the dynamics of the density in active systems such as micro swimmers Menzel
_et al._ (2016). However, care must be taken in the case of shear flow, since
for this specific choice of flow, the advection term in the DDFT necessarily
vanishes Brader and Krüger (2011). This breakdown of the DDFT is traced back
to the adiabatic assumption on the higher order ($n$-body) densities (assumed
to be in equilibrium), and is corrected by considering exact solutions to the
pair Smoluchowski equation for the correlation function in the presence of an
external flow, for small flow rate Scacchi, Krüger, and Brader (2016).
Additionally, there have been many extensions to DDFT specific to shear flow,
including: showing the importance of HI in computing the microscopic stress
tensor Aerov and Krüger (2015), the investigation of laning instabilities
Scacchi, Archer, and Brader (2017), and the emergence of nonequilibrium phase
transitions depending on shear rate Stopper and Roth (2018).

All of the above DDFTs treat the dynamics of the colloidal particles as being
overdamped (non-inertial), and there is a distinct shortcoming when it comes
to the study of inertial effects in the presence of external flows. In this
paper we present a generalisation of Refs. Rauscher  _et al._ , 2007;
Rauscher, 2010; Almenar and Rauscher, 2011 to include the effect of inertia
and HI. We show that the theory reduces to previous DDFTs in the overdamped
limit. Our DDFT is derived _ab initio_ , starting from the appropriately
stated Langevin equations rigorously derived by Dobson  _et al._ (2013). We
also present some numerical solutions demonstrating the applicability of the
derived DDFT to physically motivated problems.

##  II Nonequilibrium Langevin Dynamics

Our starting point is a thermostated nonequilibrium dynamics for colloidal
particles suspended in the flowing solvent (see, e.g., Dobson  _et al._ (2013,
Equation (3.1), Section 3)). In particular, we consider the system of
interacting stochastic differential equations (SDEs) on $\mathbb{R}^{3}$,
which govern the positions $\mathbf{r}_{i}$ and momenta $\mathbf{p}_{i}$ of
$i=1,\ldots,N$ hard, spherically symmetric colloidal particles in the presence
of a flowing (incompressible) background field of many more and much smaller
particles given by

$\displaystyle\frac{\mathrm{d}\mathbf{r}_{i}}{\mathrm{d}t}$
$\displaystyle=\frac{1}{m}\mathbf{p}_{i},$ (1a)
$\displaystyle\frac{\mathrm{d}\mathbf{p}_{i}}{\mathrm{d}t}$
$\displaystyle=-\nabla_{\mathbf{r}_{i}}V(\mathbf{r}^{N},t)-\gamma(\mathbf{p}_{i}-m\mathbf{u}(\mathbf{r}_{i},t))+\sigma\,\mathbf{f}_{i}(t),$
(1b)

where $m$ is the particle mass, $V$ is a potential energy (including one body
and multi-body interactions), $\mathbf{u}$ is the inhomogeneous solvent flow
(which may depend on both space and time), $\gamma$ is the friction
coefficient, $k_{\text{B}}$ is Boltzmann’s constant, $T$ is temperature, and
the constant $\sigma$ satisfies the fluctuation-dissipation relation
$\sigma=(mk_{\text{B}}T\gamma)^{1/2}$. Additionally,
$\mathbf{f}_{i}(t)=(\zeta^{x}_{i}(t),\zeta^{y}_{i}(t),\zeta^{z}_{i}(t))^{\top}$
is a Gaussian white noise term with mean $\langle\zeta_{i}^{a}(t)\rangle=0$
and autocorrelation
$\langle\zeta_{i}^{a}(t)\zeta_{j}^{b}(t)\rangle=2\delta_{ij}\delta^{ab}\delta(t-t^{\prime})$.

The two conditions on the well-posedness of equations (1a)–(1b) in $d$
dimensions are that (i) $\mathbf{p}\cdot\nabla\mathbf{u}$ is uniformly
Lipschitz continuous on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ and (ii)
$\nabla_{\mathbf{p}}\cdot(\mathbf{p}\cdot\nabla\mathbf{u})=\operatorname{Tr}\nabla\mathbf{u}=\nabla\cdot\mathbf{u}=0$.
The first we regard as being mathematically reasonable and the second we note
is true for incompressible background flows. Finally, by including a potential
$V$, the stochastic equations (1a)–(1b) are equivalent to the g-SLLOD
equations in the inertial reference frame at zero temperature $T=0$, and
Dobson  _et al._ (2013) provide a derivation of the thermostated g-SLLOD
equations from a heat bath model for a single large particle.

We note that the stochastic equations (1a)–(1b) have subtle implications on
the existence of steady states. Firstly, as discussed in Dobson  _et al._
(2013, Sec 2.4.3), the stochastic equations (1a)–(1b) possess the invariant
measure
$f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1})\propto\exp(-(2mk_{B}T)^{-1}|\mathbf{p}_{1}-m\mathbf{u}|^{2})$
_only_ with the inclusion of an acceleration term
$\mathbf{p}\cdot\nabla\mathbf{u}$ (arising from the flowing solvent’s local
strain rate) in the momentum equation. In this case however, the resulting
dynamics can be seen as an inertial frame transformation of the classical
static bath case (see, e.g., Ref. Archer, 2009). In contrast, an analytical
expression for stationary states of (1a)–(1b), are not known for general
choice of $\mathbf{u}$ (other than the trivial choice $\mathbf{u}=0$), and
hence we consider these Langevin dynamics in order to obtain non-trivial
solutions not obtainable by frame transformation. Secondly, no assumptions on
$\mathbf{u}$ (other than Lipschitz continuity and incompressibility) are
required. In particular, as we will see, solvent flows with finite vorticity
are permissible; this motivates us to believe that the present analysis is of
particular interest since the final equations of motion are not restricted to,
e.g., potential flows which may be trivially absorbed into an external one-
body potential. In general we assume that $\mathbf{u}$ is not the gradient of
some potential, and can have non-zero curl. Additionally, the angular momenta
of the individual colloids are assumed negligible. To go beyond the overdamped
assumption on the angular degrees of freedom of the colloids, one must couple
to (1a)–(1b) a dynamical equation for the individual colloid angular
acceleration, which has been presented before for orientable particles, albeit
without external solvent flows; see e.g., Ref. Durán-Olivencia, Goddard, and
Kalliadasis, 2016.

### Potential Energy

We use the standard many-body expansion for the potential

$\displaystyle
V(\mathbf{r}^{N},t)=\sum_{i=1}^{N}V_{1}(\mathbf{r}_{i},t)+\sum_{n=2}^{N}\frac{1}{n!}\sum_{i_{1}\neq\cdots\neq
i_{n}=1}^{N}V_{n}(\mathbf{r}_{i_{1}},\cdots\mathbf{r}_{i_{n}}),$ (2)

where $V_{1}$ is a one-body external potential and $V_{n}$ are the $n$-body
interparticle potentials. We assume that these many-body interactions have no
explicit time dependence.

### Overdamped Dynamics

It can be observed, by the usual heuristic arguments, that overdamped dynamics
may be reproduced by assuming that the acceleration of a colloid in the flow
is zero, $\mathrm{d}\mathbf{p}_{i}/\mathrm{d}t=0$, such that

$\displaystyle 0=$
$\displaystyle-\nabla_{\mathbf{r}_{i}}V(\mathbf{r}^{N},t)-\gamma(\mathbf{p}_{i}-m\mathbf{u}(\mathbf{r}_{i},t))+\sigma\,\mathbf{f}_{i}(t),$
(3)

which, after substituting the corresponding value of $\mathbf{p}_{i}$ into
(1a), gives the overdamped Langevin dynamics

$\displaystyle\frac{\mathrm{d}\mathbf{r}_{i}}{\mathrm{d}t}$
$\displaystyle=\mathbf{u}(\mathbf{r}_{i},t)-\frac{1}{m\gamma}\nabla_{\mathbf{r}^{N}}V(\mathbf{r}^{N},t)+\frac{\sigma}{m\gamma}\,\mathbf{f}_{i}.$
(4)

The dynamics (4) are recovered for DDFTs modelling colloids in a flowing
solvent driven externally by a velocity potential in the overdamped regime,
e.g., Rauscher  _et al._ (2007), Rauscher (2010) but here without the
assumption that $\mathbf{u}$ is the gradient of a potential. The overdamped
limit may be performed more rigorously through an adiabatic elimination
process that is based on a multiscale analysis applied to the Fokker-Planck
equation associated to Eqs. (1a)–(1b) (see Ref. Goddard, Pavliotis, and
Kalliadasis, 2012 for a version without a background flow but including
hydrodynamic interactions).

##  III Towards A DDFT With External Flow

The Fokker-Planck equation associated to the Langevin dynamics (1a)–(1b) is
the partial differential equation (PDE)

$\displaystyle\partial_{t}f^{(N)}(\mathbf{r}^{N},\mathbf{p}^{N},t)+\frac{1}{m}\sum_{i=1}^{N}\mathbf{p}\cdot\nabla_{\mathbf{r}_{i}}f^{(N)}(\mathbf{r}^{N},\mathbf{p}^{N},t)$
$\displaystyle\quad-\sum_{i=1}^{N}\nabla_{\mathbf{r}_{i}}V(\mathbf{r}^{N},t)\cdot\nabla_{\mathbf{p}_{i}}f^{(N)}(\mathbf{r}^{N},\mathbf{p}^{N},t)$
$\displaystyle=\sum_{i=1}^{N}\nabla_{\mathbf{p}_{i}}\cdot\Big{[}\gamma\Big{(}\mathbf{p}_{i}-m\mathbf{u}(\mathbf{r}_{i},t)+mk_{\text{B}}T\nabla_{\mathbf{p}_{j}}\Big{)}f^{(N)}(\mathbf{r}^{N},\mathbf{p}^{N},t)\Big{]}$
(5)

where $f^{(N)}(\mathbf{r}^{N},\mathbf{p}^{N},t)$ is the probability density of
finding each particle $i$ at position $\mathbf{r}_{i}$ with momentum
$\mathbf{p}_{i}$ at time $t$, referred to as the $N$-body density. Equation
(5) is a high dimensional PDE; indeed the number of discretisation points in a
standard numerical scheme, such as finite difference or pseudospectral
methods, increases exponentially in $N$. Hence it is numerically intractable
for all but a small number of particles.

To continue, we integrate (5) and derive equations for its moments, giving an
infinite hierarchy to be truncated with moment closure. This passage is well
discussed; see e.g. Refs. Archer, 2009; Goddard  _et al._ , 2012 for an in
depth derivation of the moment closure. Here we provide only the essential
details in order to connect the main quantities and equations; see Ref.
Goddard  _et al._ , 2012 for relevant definitions. Multiplying (5) by $N$ and
integrating over all $\mathbf{p}_{i}$ and all but one particle position,
$\mathbf{r}_{1}$, yields the continuity equation (6) for the one body density
$\varrho\equiv\varrho^{(1)}$ given by

$\displaystyle\partial_{t}\varrho(\mathbf{r}_{1},t)+\nabla_{\mathbf{r}_{1}}\cdot\mathbf{j}(\mathbf{r}_{1},t)=0,$
(6)

where $\varrho(\mathbf{r}_{1},t)$ is the one body density (Eq. (11) below,
with $n=1$) and the current $\mathbf{j}$ is defined by

$\displaystyle\mathbf{j}(\mathbf{r}_{1},t):=\int\mathrm{d}\mathbf{p}_{1}\,\frac{\mathbf{p}_{1}}{m}f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1},t).$
(7)

Multiplying (5) by $N\mathbf{p}_{1}/m$ and once again integrating over all
variables except for $\mathbf{r}_{1}$ yields the current evolution equation

$\displaystyle\partial_{t}\mathbf{j}(\mathbf{r}_{1},t)+\nabla_{\mathbf{r}_{1}}\cdot\int\mathrm{d}\mathbf{p}_{1}\,\Big{(}\frac{\mathbf{p}_{1}}{m}\otimes\frac{\mathbf{p}_{1}}{m}\Big{)}f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1},t)$
$\displaystyle\quad+\frac{1}{m}\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}V_{1}(\mathbf{r}_{1},t)$
$\displaystyle\quad+\frac{1}{m}\sum_{n=2}^{N}\int\mathrm{d}\mathbf{r}^{n-1}\nabla_{\mathbf{r}_{1}}V_{n}(\mathbf{r}^{n})\varrho^{(n)}(\mathbf{r}^{n},t)$
$\displaystyle\quad+\gamma(\mathbf{j}\big{(}\mathbf{r}_{1},t)-\varrho(\mathbf{r}_{1},t)\mathbf{u}(\mathbf{r}_{1},t)\big{)}=0,$
(8)

where $\mathrm{d}\mathbf{r}^{n-1}$ denotes
$\mathrm{d}\mathbf{r}^{2}\ldots\mathrm{d}\mathbf{r}^{n}$. To close the system
(6)–(8), equations are needed for the quantities involving
$\varrho^{(n)}(\mathbf{r}^{n},t)$ and
$f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1},t)$. For the former we employ the
standard adiabatic approximationMarconi and Tarazona (1999); Archer (2009);
Goddard  _et al._ (2012); Hansen and McDonald (2013); te Vrugt, Löwen, and
Wittkowski (2020); Schmidt (2022), while for the latter we use a modified
local equilibrium (more properly, local steady state) approximation, which
takes into account the background flow.

###  III.1 Adiabatic Approximation

We consider the free energy of the non-equilibrium system at constant
temperature and introduce the Helmholtz free energy functionalHansen and
McDonald (2013); te Vrugt, Löwen, and Wittkowski (2020); Evans (1979)

$\displaystyle\mathcal{F}[\varrho]$
$\displaystyle=k_{B}T\int\mathrm{d}\mathbf{r}_{1}\,\varrho(\log(\Lambda^{3}\varrho)-1)+\int\mathrm{d}\mathbf{r}_{1}\,\varrho
V_{1}(\mathbf{r}_{1},t)$ $\displaystyle\quad+\mathcal{F}_{ex}[\varrho].$ (9)

In the above functional $\Lambda$ is the de Broglie wavelength (which is
superfluous) and $\mathcal{F}_{ex}[\varrho]$ is the excess over ideal gas term
which is, in general, unknown; we will later make use of an equilibrium
identity (sum rule) connecting it to the many-body parts of the potential.
First we introduce the reduced phase space densities, defined by

$\displaystyle f^{(n)}(\mathbf{r}^{n},\mathbf{p}^{n},t)$
$\displaystyle:=\frac{N!}{(N-n)!}\int\mathrm{d}\mathbf{r}^{N-n}\,\mathrm{d}\mathbf{p}^{N-n}\,f^{(N)}(\mathbf{r}^{N},\mathbf{p}^{N},t),$
(10)

and the reduced configuration densities $\varrho^{(n)}(\mathbf{r}^{n},t)$,
given by

$\displaystyle\varrho^{(n)}(\mathbf{r}^{n},t)$
$\displaystyle:=\int\mathrm{d}\mathbf{p}^{n}\,f^{(n)}(\mathbf{r}^{n},\mathbf{p}^{n},t),$
(11)

where, again, $\varrho^{(1)}\equiv\varrho$ and
$\mathrm{d}\mathbf{x}^{N-n}=\mathrm{d}\mathbf{x}^{n+1}\ldots\mathrm{d}\mathbf{x}^{N}$.
It is assumed that the $\varrho^{(n)}(\mathbf{r}^{n},t)$ are well approximated
by their counterpart $n$–body densities for an equilibrium fluid and we will
now determine the required steady sum rule. To establish a well defined
equilibrium density distribution we must assume $\mathbf{u}\equiv 0$, which,
e.g., can be viewed as a valid initial condition for the Fokker–Planck
equation (5). In particular, at equilibrium, the one body phase space function
is assumed to be a Maxwell–Boltzmann distribution

$\displaystyle
f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1}):=\tfrac{\varrho(\mathbf{r}_{1})}{(2\pi
mk_{B}T)^{3/2}}\exp\left(-\tfrac{|\mathbf{p}_{1}|^{2}}{2mk_{B}T}\right).$ (12)

Under this assumption, we see that, as the momentum distribution equilibrates,
the kinetic pressure term may be calculated by

$\displaystyle\nabla_{\mathbf{r}_{1}}\cdot\int\mathrm{d}\mathbf{p}_{1}\,\Big{(}\frac{\mathbf{p}_{1}}{m}\otimes\frac{\mathbf{p}_{1}}{m}\Big{)}f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1})=\tfrac{k_{B}T}{m}\nabla_{\mathbf{r}_{1}}\varrho(\mathbf{r}_{1},t).$
(13)

Hence, at steady flow, using equations (12) and (13), equation (8) becomes

$\displaystyle\tfrac{k_{B}T}{m}\nabla_{\mathbf{r}_{1}}\varrho(\mathbf{r}_{1},t)+\frac{1}{m}\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}V_{1}(\mathbf{r}_{1},t)$
$\displaystyle\quad+\frac{1}{m}\sum_{n=2}^{N}\int\mathrm{d}\mathbf{r}^{n-1}\nabla_{\mathbf{r}_{1}}V_{n}(\mathbf{r}^{n})\varrho^{(n)}(\mathbf{r}^{n},t)=0.$
(14)

Now, by computing the gradient of the variational derivative of $\mathcal{F}$,
the Euler–Lagrange equation for the steady state density is

$\displaystyle 0$
$\displaystyle=\frac{1}{m}\varrho\nabla_{\mathbf{r}_{1}}\frac{\delta\mathcal{F}[\varrho]}{\delta\varrho}$
$\displaystyle=\frac{k_{B}T}{m}\nabla_{\mathbf{r}_{1}}\varrho+\frac{1}{m}\varrho\nabla_{\mathbf{r}_{1}}V_{1}+\frac{1}{m}\varrho\nabla_{\mathbf{r}_{1}}\frac{\delta\mathcal{F}_{\text{ex}}[\varrho]}{\delta\varrho}.$
(15)

By subtracting (14) from (15) one obtains the steady sum ruleEvans (1979);
Archer and Evans (2004)

$\displaystyle\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}_{\text{ex}}}{\delta\varrho}[\varrho]$
$\displaystyle=\sum_{n=2}^{N}\int\mathrm{d}\mathbf{r}^{n-1}\nabla_{\mathbf{r}_{1}}V_{n}(\mathbf{r}^{n})\varrho^{n}(\mathbf{r}^{n},t).$
(16)

We assume that (16) holds out of equilibrium to complete the adiabatic
approximation on the higher $n$-body densities.

###  III.2 Moment Closure and Local Equilibrium Approximation

Using the definition of $\mathcal{F}$ along with the (out of equilibrium)
approximation (16), equation (8) becomes

$\displaystyle\partial_{t}\mathbf{j}(\mathbf{r}_{1},t)+\bm{A}(\mathbf{r}_{1},t)+\frac{1}{m}\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}\frac{\delta\mathcal{F}[\varrho]}{\delta\varrho}$
$\displaystyle\quad+\gamma\big{(}\mathbf{j}(\mathbf{r}_{1},t)-\varrho(\mathbf{r}_{1},t)\mathbf{u}(\mathbf{r}_{1},t)\big{)}=0,$
(17)

where

$\displaystyle\bm{A}(\mathbf{r}_{1},t)$
$\displaystyle=\nabla_{\mathbf{r}_{1}}\cdot\int\mathrm{d}\mathbf{p}_{1}\,\Big{(}\frac{\mathbf{p}_{1}}{m}\otimes\frac{\mathbf{p}_{1}}{m}-\frac{k_{B}T}{m}\bm{1}\Big{)}f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1},t).$
(18)

Note that at steady state, $\bm{A}$ converges to zero. Away from equilibrium,
we assume that the momentum of each colloid particle is distributed according
to the Maxwell–Boltzmann distribution centred at the mean local velocity
$\mathbf{v}(\mathbf{r}_{1},t)$. We write $f^{(1)}$ as the sum of local steady
and non-steady parts

$\displaystyle
f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1},t)=f^{(1)}_{\text{ls}}(\mathbf{r}_{1},\mathbf{p}_{1},t)+f^{(1)}_{\text{ns}}(\mathbf{r}_{1},\mathbf{p}_{1},t),$
(19)

where

$\displaystyle
f^{(1)}_{\text{ls}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=\tfrac{\varrho(\mathbf{r}_{1},t)}{(2\pi
mk_{B}T)^{3/2}}\exp\left(-\tfrac{|\mathbf{p}_{1}-m\mathbf{v}(\mathbf{r}_{1},t)|^{2}}{2mk_{B}T}\right).$
(20)

The first few moments of $f^{(1)}_{\text{ls}}$ are

$\displaystyle\int\mathrm{d}\mathbf{p}_{1}\,f^{(1)}_{\text{ls}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=\varrho(\mathbf{r}_{1},t),$
(21)
$\displaystyle\int\mathrm{d}\mathbf{p}_{1}\,\mathbf{p}_{1}f^{(1)}_{\text{ls}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=m\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t),$
(22)
$\displaystyle\int\mathrm{d}\mathbf{p}_{1}\,|\mathbf{p}_{1}-m\mathbf{v}(\mathbf{r}_{1},t)|^{2}f^{(1)}_{\text{ls}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=mk_{B}T\varrho(\mathbf{r}_{1},t)$
(23)

and we impose the integral restrictions

$\displaystyle\int\mathrm{d}\mathbf{p}_{1}\,f^{(1)}_{\text{ns}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=0,$
(24)
$\displaystyle\int\mathrm{d}\mathbf{p}_{1}\,\mathbf{p}_{1}f^{(1)}_{\text{ns}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=0,$
(25)
$\displaystyle\int\mathrm{d}\mathbf{p}_{1}\,|\mathbf{p}_{1}-m\mathbf{v}(\mathbf{r}_{1},t)|^{2}f^{(1)}_{\text{ns}}(\mathbf{r}_{1},\mathbf{p}_{1},t)=0.$
(26)

By using the integral conditions (22) and (25), combined with the expansion
(19), the continuity equation (6) becomes

$\displaystyle\partial_{t}\varrho(\mathbf{r}_{1},t)+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t)\right)=0.$
(27)

We may also compute $\mathbf{A}(\mathbf{r}_{1},t)$ explicitly as

$\displaystyle\bm{A}(\mathbf{r}_{1},t)$
$\displaystyle=\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t)\otimes\mathbf{v}(\mathbf{r}_{1},t))\right)$
$\displaystyle\quad+\nabla_{\mathbf{r}_{1}}\cdot\int\mathrm{d}\mathbf{p}_{1}\,\frac{\mathbf{p}_{1}\otimes\mathbf{p}_{1}}{m^{2}}f^{(1)}_{\text{ns}}(\mathbf{r}_{1},\mathbf{p}_{1},t).$
(28)

We insert this expression for $\mathbf{A}(\mathbf{r}_{1},t)$ into (17) to
obtain

$\displaystyle\partial_{t}\left(\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t)\right)+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t)\otimes\mathbf{v}(\mathbf{r}_{1},t)\right)$
$\displaystyle\,+\nabla_{\mathbf{r}_{1}}\cdot\int\mathrm{d}\mathbf{p}_{1}\,\frac{\mathbf{p}_{1}\otimes\mathbf{p}_{1}}{m^{2}}f^{(1)}_{\text{ns}}(\mathbf{r}_{1},\mathbf{p}_{1},t)$
$\displaystyle\,+\frac{1}{m}\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}[\varrho]}{\delta\varrho}+\gamma\varrho(\mathbf{r}_{1},t)(\mathbf{v}(\mathbf{r}_{1},t)-\mathbf{u}(\mathbf{r}_{1},t))=0.$
(29)

Now consider the following identity, for an arbitrary, suitably well behaved
vector field $\mathbf{z}(\mathbf{r}_{1},t)$

$\displaystyle\partial_{t}\left(\varrho\mathbf{z}\right)+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho\mathbf{z}\otimes\mathbf{z}\right)$
$\displaystyle=\mathbf{z}\partial_{t}\varrho+\varrho\partial_{t}\mathbf{z}+\varrho\left(\mathbf{z}(\nabla\cdot\mathbf{z})+(\mathbf{z}\cdot\nabla)\mathbf{z}\right)+\mathbf{z}\left(\mathbf{z}\cdot\nabla\varrho\right).$
(30)

This result follows from the product rule of differentiation and, e.g.,
equation (27) of Ref. Archer, 2009. By setting
$\mathbf{z}=\mathbf{v}(\mathbf{r}_{1},t)$, using the continuity equation (27),
and neglecting the $f_{\text{ns}}$ term, (since it can be shown to be small,
at least in the overdamped limit Goddard, Pavliotis, and Kalliadasis (2012)
$\gamma\to\infty$), we obtain our equations of motion.

$\displaystyle\begin{cases}&\qquad\quad\partial_{t}\varrho(\mathbf{r}_{1},t)+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t)\right)=0,\\\\[5.69054pt]
&\partial_{t}\mathbf{v}(\mathbf{r}_{1},t)+(\mathbf{v}(\mathbf{r}_{1},t)\cdot\nabla_{\mathbf{r}_{1}})\mathbf{v}(\mathbf{r}_{1},t)\\\
&\qquad\quad+\tfrac{1}{m}\nabla_{\mathbf{r}_{1}}\frac{\delta\mathcal{F}[\varrho]}{\delta\varrho}(\mathbf{r}_{1},t)+\gamma\mathbf{w}(\mathbf{r}_{1},t)=0,\end{cases}$
(31)

where
$\mathbf{w}(\mathbf{r}_{1},t):=\mathbf{v}(\mathbf{r}_{1},t)-\mathbf{u}(\mathbf{r}_{1},t)$.
The momentum equation in (31) should be compared with the inertial DDFT given
in Eq. (30) of Ref. Archer, 2009. We note that the correct inertial dynamics
(31) cannot be obtained from the static bath case therein simply by a
transformation between the fixed and moving frame of reference (moving with
the streaming velocity $\mathbf{u}$). This comes as a consequence of the
nonlinear advection term in (31). One can see this by transforming the static
bath case in Ref. Archer, 2009 to the moving reference frame, which is
equivalent to shifting the transport of the velocity distribution
$\mathbf{v}(\mathbf{r}_{1},t)$ by the external flow
$\mathbf{u}(\mathbf{r}_{1},t)$. We see that the velocity advection transforms
to

$\displaystyle
D_{t}\mathbf{v}(\mathbf{r}_{1},t)\to\partial_{t}\mathbf{w}(\mathbf{r}_{1},t)+(\mathbf{w}(\mathbf{r}_{1},t)\cdot\nabla_{\mathbf{r}_{1}})\mathbf{w}(\mathbf{r}_{1},t),$
(32)

which differs from the advection term in (31). We note that it is possible,
via a careful choice of $\mathbf{u}(\mathbf{r}_{1},t)$, to make the steady
state velocity and density distribution resulting from (31) and a
frame–transformed Eq. (30) in Ref. Archer, 2009 to be made the same.
Nevertheless, the passage to equilibrium will clearly not be the same, due to
the difference in the nonlinear advection terms.

This concludes an _ab initio_ derivation of the inertial DDFT in the presence
of an inhomogeneous, incompressible background flow. We reemphasise that (31)
is not restricted to curl-free $\mathbf{u}$, and may therefore not be
trivially derived by appropriately perturbing the one-body field $V_{1}$ to
include a velocity potential accounting for the flowing solvent.

###  III.3 Overdamped Limit

We now consider the overdamped limit by attempting to close (31), into a
single equation for the dynamics of the density alone. Returning to Eqs. (6)
and (17), we have under the local-steady approximation, the conservation of
momentum

$\displaystyle\partial_{t}\mathbf{j}+\bm{A}_{ls}+\frac{1}{m}\varrho\nabla_{\mathbf{r}_{1}}\frac{\delta\mathcal{F}[\varrho]}{\delta\varrho}+\gamma(\mathbf{j}-\varrho\mathbf{u})=0,$
(33)

where $\bm{A}_{ls}$ is given by (28) with $f_{ns}=0$. Differentiating (6) we
obtain

$\displaystyle\partial^{2}_{t}\varrho+\partial_{t}\left(\nabla_{\mathbf{r}_{1}}\cdot\mathbf{j}\right)=0.$
(34)

Additionally, by taking the divergence of (33), we have

$\displaystyle\partial_{t}\left(\nabla_{\mathbf{r}_{1}}\cdot\mathbf{j}\right)+\nabla_{\mathbf{r}_{1}}\cdot\bm{A}_{ls}+\tfrac{1}{m}\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}[\varrho]}{\delta\varrho}\right)$
$\displaystyle\quad+\gamma\nabla_{\mathbf{r}_{1}}\cdot\mathbf{j}=0.$ (35)

Combining Eqs. (34) and (35), we obtain

$\displaystyle\partial^{2}_{t}\varrho+\gamma\left(\partial_{t}\varrho+\nabla_{\mathbf{r}_{1}}\cdot(\varrho\,\mathbf{u})\right)-\tfrac{1}{m}\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}[\varrho]}{\delta\varrho}\right)$
$\displaystyle\quad=\nabla_{\mathbf{r}_{1}}\cdot\bm{A}_{ls}.$ (36)

It remains to treat the term on the right hand side. We note that in Ref.
Goddard  _et al._ , 2012 we have $\bm{A}_{ls}=\bm{A}_{le}$, and, by the
results of Ref. Goddard, Pavliotis, and Kalliadasis, 2012, at least for a two-
body potential truncation of $\mathcal{F}_{ex}[\varrho]$, terms containing
$\bm{A}_{le}$ are negligible in the high-friction limit, $\gamma\to\infty$.
Hence, under a similar adiabatic elimination process, based on a multiscale
analysis of the Fokker-Planck equation (5), we expect that $\bm{A}_{ls}$ will
also be negligible in the presence of the solvent flow $\mathbf{u}$.
Therefore, (36) maybe approximated as

$\displaystyle\partial_{t}^{2}\varrho+\gamma\left(\partial_{t}\varrho+\nabla_{\mathbf{r}_{1}}\cdot(\varrho\,\mathbf{u})\right)-\tfrac{1}{m}\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}[\varrho]}{\delta\varrho}\right)\approx
0.$ (37)

The second-order time derivative $\partial_{t}^{2}\varrho$ can be
heuristically neglected in the high-friction limit, though its effect could be
important on short timescales. Hence, previous overdamped DDFTs with external
flow, e.g., Rauscher  _et al._ (2007); Krüger and Rauscher (2007); Rauscher
(2010); Almenar and Rauscher (2011), are recovered. Assuming $\varrho$ is
strictly positive, we may also establish that the free energy is a
monotonically decreasing function of time along the solvent flow

$\displaystyle D_{t}\mathcal{F}$ $\displaystyle=$
$\displaystyle\int\left(\partial_{t}\varrho+\mathbf{u}\cdot\nabla\varrho\right)\tfrac{\delta\mathcal{F}}{\delta\varrho}\mathrm{d}\mathbf{r}$
(38) $\displaystyle=$
$\displaystyle-\tfrac{1}{m\gamma}\int\varrho|\nabla\tfrac{\delta\mathcal{F}}{\delta\varrho}|^{2}\mathrm{d}\mathbf{r}\leq
0,$

where we have assumed that the boundary terms vanish; this is the case, e.g.,
when $\rho=0$ on the boundary, or when the system is equipped with for no-flux
or periodic boundary conditions. Thus, this is a generalisation of the
standard monotonicity of the free energy over time in overdamped DDFT.

##  IV Including Hydrodynamic Interactions

Hydrodynamic interactions (HIs) between the colloids may be introduced by
modifying the friction term in Eq. (31) Rex and Löwen (2008, 2009); Goddard
_et al._ (2012). By introducing the friction tensor
$\bm{\Gamma}\in\mathbb{R}^{3N\times 3N}$, the non-equilibrium Langevin
dynamics (1a)–(1b) becomes

$\displaystyle\frac{\mathrm{d}\mathbf{r}_{i}}{\mathrm{d}t}$
$\displaystyle=\frac{1}{m}\mathbf{p}_{i}$ (39a)
$\displaystyle\frac{\mathrm{d}\mathbf{p}_{i}}{\mathrm{d}t}$
$\displaystyle=-\nabla_{\mathbf{r}_{i}}V(\mathbf{r}^{N},t)-\sum_{j=1}^{N}\bm{\Gamma}_{ij}(\mathbf{r}^{N})\big{(}\mathbf{p}_{j}-m\mathbf{u}(\mathbf{r}_{j},t)\big{)}$
(39b)
$\displaystyle\quad+\sum_{j=1}^{N}\bm{B}_{ij}(\mathbf{r}^{N})\mathbf{f}_{j}(t)$
(39c)

where $\bm{B}=\left(mk_{\text{B}}T\bm{\Gamma}\right)^{1/2}$ and $\bm{\Gamma}$
comprises $N^{2}$ positive definite resistance matrices $\bm{\Gamma}_{ij}$
given by $\bm{\Gamma}_{ij}=\gamma\bm{1}+\gamma\tilde{\bm{\Gamma}}_{ij}$, where
$\bm{1}$ is the $3\times 3$ identity matrix, and
$\tilde{\bm{\Gamma}}_{ij}\in\mathbb{R}^{3\times 3}$. Physically we typically
regard the forces introduced by $\gamma\tilde{\bm{\Gamma}}_{ij}$ as short
range HIs mediated by the bath, (also known as lubrication forces) and
$\gamma$ corresponds to Stokes drag at infinite particle separation. We remark
that (39a)–(39c) are the natural generalisation to include HIs but we stress
that these equations have not been rigorously derived. However, as we will
see, they do lead to the expected set of continuum equations, consistent with
(31) and (37) upon taking appropriate limits, $\bm{Z}_{1}=\bm{Z}_{2}=\bm{0}$
(defined below) and $\gamma\to\infty$ respectively.

After carefully repeating the derivation from Section III onwards, we use the
Enskog approximation to resolve $f^{(2)}$ [see Eq. (10)]

$\displaystyle
f^{(2)}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{p}_{1},\mathbf{p}_{2},t)=f^{(1)}(\mathbf{r}_{1},\mathbf{p}_{1},t)f^{(1)}(\mathbf{r}_{2},\mathbf{p}_{2},t)g(\mathbf{r}_{1},\mathbf{r}_{2},[\varrho]),$
(40)

where $g$ is the pair correlation function, which must be supplied by
auxiliary means. Additionally we restrict the HIs to be two-body, thus

$\displaystyle\tilde{\bm{\Gamma}}_{ij}(\mathbf{r}^{N})=\delta_{ij}\sum_{l\neq
i}\bm{Z}_{1}(\mathbf{r}_{i},\mathbf{r}_{l})+(1-\delta_{ij})\bm{Z}_{2}(\mathbf{r}_{i},\mathbf{r}_{j}),$
(41)

where $\bm{Z}_{1}$ and $\bm{Z}_{2}$ are the diagonal and off diagonal blocks
respectively of the translational component of the grand resistance matrix
originating in the classical theory of low Reynolds’ number hydrodynamics
between suspended particles (see, e.g., Jeffrey and Onishi (1984), Happel and
Brenner (2012)). The corresponding DDFT is then

$\displaystyle\begin{cases}&\qquad\quad\partial_{t}\varrho(\mathbf{r}_{1},t)+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\mathbf{v}(\mathbf{r}_{1},t)\right)=0,\\\\[5.69054pt]
&\partial_{t}\mathbf{v}(\mathbf{r}_{1},t)+(\mathbf{v}(\mathbf{r}_{1},t)\cdot\nabla_{\mathbf{r}_{1}})\mathbf{v}(\mathbf{r}_{1},t)\\\
&\quad+\tfrac{1}{m}\nabla_{\mathbf{r}_{1}}\frac{\delta\mathcal{F}[\varrho]}{\delta\varrho}(\mathbf{r}_{1},t)+\gamma\mathbf{w}(\mathbf{r}_{1},t)\\\
&\quad+\gamma\int\mathrm{d}\mathbf{r}_{2}\,\Big{[}\bm{Z}_{1}(\mathbf{r}_{1},\mathbf{r}_{2})\mathbf{w}(\mathbf{r}_{1},t)+\bm{Z}_{2}(\mathbf{r}_{1},\mathbf{r}_{2})\mathbf{w}(\mathbf{r}_{2},t))\Big{]}\\\
&\qquad\qquad\times\varrho(\mathbf{r}_{2},t)g(\mathbf{r}_{1},\mathbf{r}_{2},[\varrho])=0.\end{cases}$
(42)

We note that by setting the solvent velocity to zero we obtain the
previously–derived DDFT including HIs in the case of a static background
bathGoddard  _et al._ (2012). In addition, setting the HI tensors to zero
($\bm{Z}_{1}=\bm{Z}_{2}=\bm{0}$) recovers the results of Ref. Archer, 2009.

###  IV.1 Overdamped Limit

We now consider the overdamped limit with the inclusion of HI. Assuming that
inter-particle HI are weak, that is $\bm{Z}_{2}\approx 0$, we obtain

$\displaystyle\tfrac{1}{m}\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}[\varrho]}{\delta\varrho}(\mathbf{r}_{1},t)+\gamma\bm{D}^{-1}(\mathbf{r}_{1},[\varrho],t)(\mathbf{v}(\mathbf{r}_{1},t)-\mathbf{u}(\mathbf{r}_{1},t))\approx
0.$ (43)

where $\bm{D}$ is the diffusion tensor, known rigorously to be positive
definite (see, e.g., Ref. Goddard, Pavliotis, and Kalliadasis, 2012) given by

$\displaystyle\bm{D}(\mathbf{r}_{1},[\varrho],t):=\left[{\mathbf{1}}+\int\mathrm{d}\mathbf{r}_{2}\,\bm{Z}_{1}(\mathbf{r}_{1},\mathbf{r}_{2})\varrho(\mathbf{r}_{2},t)g(\mathbf{r}_{1},\mathbf{r}_{2},[\varrho]\right]^{-1}.$
(44)

Now by adding and subtracting the same term to the continuity equation

$\displaystyle\partial_{t}\varrho+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{v}-\mathbf{u})+\varrho\mathbf{u}\right)=0,$
(45)

and combining (43) and (45) we obtain

$\displaystyle\partial_{t}\varrho(\mathbf{r}_{1},t)+\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\mathbf{u}\right)$
$\displaystyle\quad\approx\tfrac{1}{m\gamma}\nabla_{\mathbf{r}_{1}}\cdot\left(\varrho(\mathbf{r}_{1},t)\bm{D}(\mathbf{r}_{1},[\varrho],t)\nabla_{\mathbf{r}_{1}}\tfrac{\delta\mathcal{F}[\varrho]}{\delta\varrho}\right),$
(46)

which may be compared with the DDFT in Eq. (27) of Ref. Rauscher  _et al._ ,
2007. Equation (46) is the Smoluchowski equation, differing from the one found
in Rauscher  _et al._ (2007) and other approaches that start from overdamped
Langevin dynamics. We note that previous work has demonstrated that such
differences arise even in the case when there is no background flow Goddard,
Pavliotis, and Kalliadasis (2012).

##  V Results From Numerical Solution

Figure 1: Time snapshots of the density profile predicted by the DDFT (31),
with time increasing from top to bottom, as indicated on each subplot, and
external potential given in Eq. (50), with $(x_{a},\,y_{a})=(-3,\,0)$ and
$(x_{b},\,y_{b})=(3,\,0)$. In each case, we display both a surface plot and a
contour plot, with the local flux vector $\mathbf{j}$ shown by black arrows
displayed on top of the contour plots. The four plots on the left are for the
case of zero background flow $\mathbf{u}=0$, and the four on the right are for
a line source, i.e., uniform flow with $\mathbf{u}$ given by Eqs. (47a)–(47b).

The aim of this section is to illustrate that the inclusion of general,
incompressible background flows in our theory is a non-trivial addition to
DDFT, and can, in fact, lead to interesting differences in the dynamics of the
density and velocity of colloidal dispersions when compared to non-driven
systems. We also wish to demonstrate that the DDFTs described above are
numerically tractable, at least in two spatial dimensions. We do this by
presenting a selection of numerical computations of the derived DDFT (31) for
a system of hard spheres with excess Helmholtz free energy functional
$\mathcal{F}_{ex}[\varrho]$ given by fundamental measure theory (FMT), and
compare the resulting dynamics to previously derived DDFTs, specifically the
inertial DDFT of Ref. Archer, 2009 and the overdamped with external flow
theory of Ref. Rauscher  _et al._ , 2007. The choice of excess free energy
functional to model hard spheres is not restrictive, and, in particular, our
DDFT can also be used to describe the flow of other systems through the choice
of a different approximation for $\mathcal{F}_{ex}[\varrho]$. For example, we
have also tested our DDFT in conjunction with a simple mean-field
approximation for $\mathcal{F}_{ex}[\varrho]$ relevant to soft-core particles,
but do not present the results here.

All of the numerical simulations presented here are performed using
pseudospectral methods, available in the MATLAB 2DChebClass package Goddard,
Nold, and Kalliadasis (2017).

The first example we present concerns a solvent flow taking the form of a line
source in a uniform field that is inspired by two-dimensional irrotational
flow from classical fluid mechanics. The second example considers an external
flow field with finite vorticity and a local sink, thereby demonstrating the
general validity of the presented DDFT when subject to flows with a rotational
component as well as the effects of inertia. Such rotational flows can be
considered as a model for forced-flow layer liquid chromatography Sherma and
Fried (2003), whereas the addition of a local sink models the flow of colloids
around a plug-hole-like exit. The third example concerns flow in a confined
geometry (a periodic two-dimensional capillary) for a simple model system in
haemodynamics Yu  _et al._ (2017). Here we demonstrate the effect of inertia
in response to an instantaneous change in the direction of the applied
background flow, showing agreement with Ref. Rauscher  _et al._ , 2007 both in
the high friction limit, and in equilibrium.

For all our numerical experiments we non-dimensionalize the equations with the
units of length, mass and energy being $\varsigma$ (colloid diameter), $m$
(colloid mass) and $k_{B}T$ (Boltzmann’s constant times the temperature),
respectively. For examples in unbounded domains, we also artificially modify
the solvent velocity in the far field, in order to obtain finite solutions to
our computations. This is necessary when using a pseudo-spectral collocation
scheme, where on unconfined domains there exist computational points at
infinity. This is only done in regions of the computational domain where the
density $\rho(\mathbf{r},t)$ is infinitesimal or zero (in particular, at
infinity) and is achieved by smoothly reducing the solvent velocity to zero
away from regions with finite density (e.g., through multiplication by an
appropriately–chosen error function). The results of the simulations are
insensitive to the exact choice of this methodology.

We choose the colloidal particles to be hard spheres, with pair potential
$V_{2}(|\mathbf{r}-\mathbf{r}^{\prime}|)=\infty$ for
$|\mathbf{r}-\mathbf{r}^{\prime}|<\varsigma$ and zero otherwise, where
$\mathbf{r}$ and $\mathbf{r}^{\prime}$ denote the positions of two different
colloidal particles. This hard sphere exclusion is described by the excess
free energy $\mathcal{F}_{ex}[\varrho]$; here we use the three-dimensional FMT
functional of Rosenfeld (1989), averaged into two dimensions. This means that
though our numerical solutions are restricted to two dimensions, the physical
modelling of the hard spheres is maintained in three dimensions. Additionally,
we do not consider any Faxén (1922) type corrections to the solvent fluid
(c.f., Rauscher  _et al._ (2007); Rauscher (2010)) whereby we mean that, in
reality, the colloidal particle geometry is known to perturb the flowing
solvent field local to each particle. In effect, we assume that the solvent
flows through the particles as it carries them in the flow. We also assume
only moderately low density particle systems, due to the known deficiencies in
the local steady approximation to accurately predict the dynamics far from
steady state Almenar and Rauscher (2011).

We remark now on initial conditions: In each of the numerical examples
presented we assume $\mathbf{u}=0$ at $t=0$, but thereafter ($t>0$)
$\mathbf{u}$ instantaneously changes to a non-zero flow. The reason for this
is twofold: (i) it permits the existence of an initial density distribution
compatible with the present system of PDEs (31); (ii) it allows us to examine
the transient behaviour between the different DDFTs without including dynamic
effects arising due to different initial data. Obtaining the initial
equilibrium density distribution involves solving the DFTEvans (1979); Wu and
Li (2007) problem $(\delta\mathcal{F}[\varrho]/\delta\varrho)=0$. We obtain
the initial condition via Picard iteration, a standard technique within the
DFT community, see e.g., Roth (2010), and which can also be viewed as a damped
version of standard fixed point algorithms. The numerical examples presented
here were determined to be converged in the number of collocation points $M$.
In particular, we repeated the computations under increasing $M$, until the
time dependent dynamics appeared constant in $M$.

###  V.1 Line Source in a Uniform Flow

Figure 2: Time snapshots of the density profile predicted by the DDFT (31),
with time increasing from top to bottom, as indicated on each subplot, and
external potential given in Eq. (52), with $(x_{a},\,y_{a})=(0,\,4)$. In each
case, we display both a surface plot and a contour plot, with the local flux
vector $\mathbf{j}$ shown by black arrows displayed on top of the contour
plots. The four plots on the left are for the case of zero background flow
$\mathbf{u}=0$, and the four on the right are for the case where $\mathbf{u}$
given by (51a)–(51b).

This example is a model for colloidal particles subject to an optical trap in
the presence of a background solvent flow emanating from a point source, and
moving uniformly throughout the $x$-$y$ plane. The background flow is defined
as follows: Let $\mathbf{e}_{x}$ and $\mathbf{e}_{y}$ denote the canonical
basis vectors of the two-dimensional Euclidean plane in the $x$ and $y$
directions, respectively. We consider an external flow of the form
$\mathbf{u}(x,y)=\left(u(x,y),\,v(x,y)\right)^{\top}$ where $u(x,y)$ and
$v(x,y)$ are the scalar velocity fields of the solvent in the $\mathbf{e}_{x}$
and $\mathbf{e}_{y}$ directions respectively, given by

$\displaystyle u$ $\displaystyle=U+\frac{Qx}{2\pi(x^{2}+y^{2})},$ (47a)
$\displaystyle v$ $\displaystyle=\frac{Qy}{2\pi(x^{2}+y^{2})},$ (47b)

where $U$ and $Q$ are the magnitudes of the uniform and line flow sources,
respectively. We note that $\mathbf{u}(x,y)$ is derived from the gradient of
the potential $\phi(x,y)=Ux+Q/(2\pi)\log\sqrt{x^{2}+y^{2}}$. There is a
stagnation point (i.e. point where the velocity is zero) in the flow, at
$(x,\,y)=(-Q/(2\pi U),\,0)$, and a separatrix (crossing the $x$-axis through
this point) separates the fluid emanating from the source (at
$(x,\,y)=(0,\,0)$) from the fluid carried by the uniform flow. To illustrate
the influence of this flow field, we choose $U=0.1$ and $Q=0.35$ and use a
confining quadratic external potential $V_{1}(x,y)=\lambda(x^{2}+y^{2})$ to
force the particles to remain in the vicinity of the origin, a choice which
has been shown to not induce layering (i.e. oscillations in the particle
density profile due to packing against the boundary). For $t>0$, the external
potential $V_{1}$ is modified to

$\displaystyle V_{1}(x,y,t)=$ $\displaystyle\lambda(x^{2}+y^{2})$ (48)
$\displaystyle-\lambda
e^{-\left(\frac{t}{\tau}\right)^{2}}e^{-\left(\frac{(x-x_{a})^{2}}{\sigma^{2}}+\frac{(y-y_{a})^{2}}{\sigma^{2}}\right)}$
(49)
$\displaystyle-\lambda\left(1-e^{-\left(\frac{t}{\tau}\right)^{2}}\right)e^{-\left(\frac{(x-x_{b})^{2}}{\sigma^{2}}+\frac{(y-y_{b})^{2}}{\sigma^{2}}\right)}.$
(50)

We choose $(x_{a},\,y_{a})=(-3,\,0)$ and $(x_{b},\,y_{b})=(3,\,0)$, thereby
moving the minimum in $V_{1}$ from near $(x_{a},\,y_{a})$ to near
$(x_{b},\,y_{b})$, which drives the density across the stagnation point. We
also choose $\lambda=0.01$, $\sigma=50$, $\tau=0.01$ and $\gamma=2$ for $N=50$
hard spheres. In Figure 1 we show the transient densities
$\varrho(\mathbf{r},t)$ and velocities $\mathbf{v}(\mathbf{r},t)$ obtained
from (31), alongside the solutions obtained with external flow turned off
($\mathbf{u}=0$), for comparison. We see substantial differences between
density solutions with and without external flow. We find that the chosen
termination time for the dynamics, $t=12$, is sufficient to observe both
systems reach steady state (not shown). In the presence of external flow
(yellow-pink density plots) we observe that the density peak is advected
across the stagnation point at $(-Q/(2\pi U),\,0)\approx(-0.56,\,0)$, which
generates a dip in the density profile in the vicinity of that point.

###  V.2 Flow in a Localised Vortex and Sink

Figure 3: DDFT results for the density and current for flow in a channel. In
each row the first and third plots show the densities resulting from
overdamped and underdamped DDFT, respectively; the second plot shows the
difference between the two densities ($\rho_{o}-\rho_{u}$); the fourth plot
shows the velocity of the underdamped dynamics (blue arrows), and the
background flow (grey arrows). Each row corresponds to a different time:
$t=0,0.4,0.6,1.2,5$.

To more fully understand the effect of including an external flow under
inertial dynamics, we now study a simple model including finite vorticity. We
consider an external flow field of the form

$\displaystyle u$ $\displaystyle=R(y-x_{r_{0}})+\frac{Qx}{2\pi(x^{2}+y^{2})},$
(51a) $\displaystyle v$
$\displaystyle=-R(x-y_{r_{0}})+\frac{Qy}{2\pi(x^{2}+y^{2})},$ (51b)

where $R=0.9$ is a scalar parameter controlling the strength and direction of
the vortex, $Q=-1.57$ is the strength of the local sink, and
$(x_{r_{0}},\,y_{r_{0}})=(0,\,0)$ is the centre of rotation. This flow
diverges as $x,y\to\pm\infty$, so we include an additional multiplicative term
to decay the solvent velocity to zero such that
$\tilde{\mathbf{u}}(x,y)=\mathbf{u}(x,y)e^{-(x^{2}+y^{2})/\alpha}$, with the
choice $\alpha=8.5$. We have maintained the same external potential width
$\sigma$, time scale $\tau$, friction coefficient $\gamma$, as well as the
number of particles $N=50$ from the previous example, and once again we have
used Rosenfeld (1989) FMT as the excess over ideal gas free energy functional
$\mathcal{F}_{ex}$. For this example we use a slowly decaying confining
potential

$\displaystyle V_{1}(x,y,t)=\lambda(x^{2}+y^{2})-\lambda
e^{-\left(\frac{t}{\tau}\right)^{2}}e^{-\left(\frac{(x-x_{a})^{2}}{\sigma^{2}}+\frac{(y-y_{a})^{2}}{\sigma^{2}}\right)},$
(52)

with $\lambda=0.001$ and $(x_{a},\,y_{a})=(0,\,4)$. Here,
$V_{1}(x,y,t)\to\lambda(x^{2}+y^{2})$ as $t\to\infty$, which is a weak
background potential necessary to maintain a nonzero density over the infinite
2D plane. In Figure 2 we plot time snap-shots of the numerical solution. A
non-trivial density evolution due to the combined effect of inertia and the
external flow may be observed. We observe in the $\mathbf{u}\neq 0$ case on
the right that the flow field (51a)–(51b) drags the colloids around in a
swirling circular motion. Eventually, in the steady state, the effects of hard
sphere exclusion become apparent, with the colloids forming a stable ring
shaped density distribution around the sink location.

###  V.3 Flow in a Periodic Strip

To illustrate the effects of externally flowing solvents on the dynamics of
the density combined with geometrical confinement we consider a system of hard
spheres, of radius 0.5, confined in a channel of width 4 ($x\in[-2,2]$) and
length $2\pi$ (between $y\in[0,2\pi]$). The domain is periodic in the
$y$-direction, and we impose no-flux boundary conditions on the density and
current in the $x$-direction. In other words,
$\mathbf{j}\cdot\mathbf{e}_{x}=0$ on $x=-2$, $x=2$ and $t\geq 0$, and
$\varrho(\mathbf{r},t)=\varrho(\mathbf{r}+2\pi\mathbf{e}_{y},t)$ for
$(x,y)\in[-2,2]\times[0,2\pi]$ and $t\geq 0$. We impose a parabolic external
flow of the form

$\displaystyle u$ $\displaystyle=0,$ (53a) $\displaystyle v$
$\displaystyle=A(t)(y-W/2)(y+W/2),$ (53b)

where $W=4$ is the width of the channel, and $A(t)$ controls the strength and
direction of the flow via:

$A(t)=\begin{cases}0,\quad t=0\\\ \alpha,\quad t\in(0,t_{s}]\\\ -\alpha,\quad
t>t_{s}\end{cases}$

where $\alpha=1$ and $t_{s}=0.5$. This represents a parabolic flow which
(instantaneously) switches direction at a given time ($t_{s}$).

We choose an initial condition of the form

$\rho=Z(\cos(y+y_{0})+\rho_{0}),$

where $y_{0}=\pi$ and $\rho_{0}=2$. We choose $Z$ such that $\int\rho dxdy=4$.
Note that this is not the total number of particles in the system, but rather
a measure of the (average) area of the 2D slice containing particles; the true
number of particles in the 3D system can be arbitrarily large.

In Figure 3 we clearly see significant effects of both the channel walls
(which causes non-trivial packing of particles at the walls; note that we do
not show the excluded-volume region where the density is zero), and inertia.
The latter is demonstrated both at times before the switch of the external
flow (at $t=0.4$, it can be seen that the peak of the underdamped distribution
lags behind that of the overdamped dynamics), and shortly after the switch of
the external flow (at $t=0.6$), where the velocity predicted by the DDFT
incorporating inertial effects remains in the opposite direction to the
external flow. In contrast, the flux of the overdamped dynamics
instantaneously switches direction at $t=0.5$. For longer times, the
velocities of the inertial system align with the external flow, until, for
very long times, the two densities become indistinguishable, where the
underdamped velocity aligns with the external flow (not shown).

##  VI Summary and Outlook

We have derived a general DDFT for systems of colloidal particles, including
the effects of inertia, general time-dependent, nonhomogeneous externally
flowing solvents, and hydrodynamic interactions. Additionally, in the
overdamped limit, the total free energy in the system may be shown to be
decreasing along the external flow. For the trivial flow field $\mathbf{u}=0$,
the DDFT relaxes towards the equilibrium density given by the free energy
functional used for the dynamics. The _ab initio_ derivation of the presented
DDFT, relies on three assumptions for the novel inclusion of general
nonhomogeneous external flows: (i) the externally flowing solvent is
incompressible; (ii) in the case of including HI, the $n$–body distributions
governing the HI terms remain well approximated by functionals of
$\varrho(\mathbf{r},t)$ and $\mathbf{v}(\mathbf{r},t)$; (iii) the density of
the particle systems is relatively low. For the latter, in the overdamped
limit, Almenar and Rauscher (2011) have shown that the steady state sum rule
(16) is violated in dense systems, by computing the left hand side using the
steady state solution to their DDFT and comparing with the right hand side
computed with the equilibrium probability density
$P(\mathbf{r}_{1},\mathbf{r}_{2})$ obtained from the the Fokker–Planck
equation for the two-point correlation function for a system of exactly two
particles. It is likely that the local steady state approximation of the
inertial DDFT presented here also breaks down for dense non–equilibrium steady
states. To get around this, one could consider closing the hierarchy and
obtaining a steady state relation for the three–point correlation function
Tschopp and Brader (2022); Tschopp, Vuijk, and Brader (2023), or by using the
methods of power functional theory Schmidt (2022).

For the remainder of the derivation, we rely on the usual approximations: (i)
The non-steady contributions from the one-body distribution which are not
captured by the local-steady approximation are small or may be neglected. (ii)
For two-body HI, it is assumed that the widely used Enskog approximation is
compatible with externally flowing solvents which are subject to the
aforementioned assumptions. Further, by neglecting HI and external flows we
recover the DDFT derived in Ref. Archer, 2009, and by neglecting external
flows only we recover the DDFT derived in Ref. Goddard  _et al._ , 2012.

Via the numerical solutions presented in Section V, we have demonstrated that
the inclusion of inertia as well as general, inhomogeneous external flow
fields may be combined with the usual approximations made in deriving a DDFT
to produce accurate models for colloidal fluids undergoing external forcing,
applicable to a wide range of non-trivial systems. Promising future extensions
to the presented DDFT, both theoretical and numerical, include the extension
to multiple–particle species, anisotropic particles, self-propelled particles,
as well as a more thorough investigation of the effects of confined
geometries.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Author Declarations

The authors have no conflicts to disclose.

###### Acknowledgements.

B. D. Goddard would like to acknowledge support from EPSRC EP/L025159/1. R. D.
Mills-Williams is grateful to EPSRC for PhD funding. Both B. D. Goddard and A.
J. Archer gratefully acknowledge support for parts of this work from the
London Mathematical Society and the Loughborough University Institute of
Advanced Studies.

## References

  * Archer (2005) A. J. Archer, “Dynamical density functional theory: phase separation in a cavity and the influence of symmetry,” J. Phys.: Condens. Matter 17, S3253 (2005).
  * Goddard  _et al._ (2012) B. D. Goddard, A. Nold, N. Savva, P. Yatsyshin, and S. Kalliadasis, “Unification of dynamic density functional theory for colloidal fluids to include inertia and hydrodynamic interactions: derivation and numerical experiments,” J. Phys.: Condens. Matter 25, 035101 (2012).
  * Archer (2009) A. J. Archer, “Dynamical density functional theory for molecular and colloidal fluids: A microscopic approach to fluid mechanics,” J. Chem. Phys. 130, 014509 (2009).
  * Wittkowski and Löwen (2011) R. Wittkowski and H. Löwen, “Dynamical density functional theory for colloidal particles with arbitrary shape,” Mol. Phys. 109, 2935–2943 (2011).
  * Durán-Olivencia, Goddard, and Kalliadasis (2016) M. A. Durán-Olivencia, B. D. Goddard, and S. Kalliadasis, “Dynamical density functional theory for orientable colloids including inertia and hydrodynamic interactions,” J. Stat. Phys. 164, 785–809 (2016).
  * Durán-Olivencia  _et al._ (2017) M. A. Durán-Olivencia, P. Yatsyshin, B. D. Goddard, and S. Kalliadasis, “General framework for fluctuating dynamic density functional theory,” New Journal of Physics 19, 123022 (2017).
  * Hansen and McDonald (2013) J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids: with applications to soft matter_ (Academic press, 2013).
  * te Vrugt, Löwen, and Wittkowski (2020) M. te Vrugt, H. Löwen, and R. Wittkowski, “Classical dynamical density functional theory: from fundamentals to applications,” Adv. Physics 69, 121–247 (2020).
  * Callen and Welton (1951) H. B. Callen and T. A. Welton, “Irreversibility and generalized noise,” Phys. Rev. 83, 34 (1951).
  * Dobson  _et al._ (2013) M. Dobson, F. Legoll, T. Lelièvre, and G. Stoltz, “Derivation of Langevin dynamics in a nonzero background flow field,” ESAIM: Math. Model. Num. 47, 1583–1626 (2013).
  * Evans and Morriss (2007) D. J. Evans and G. P. Morriss,  _Statistical mechanics of nonequilbrium liquids_ (ANU Press, 2007).
  * Todd and Daivis (1999) B. D. Todd and P. J. Daivis, “A new algorithm for unrestricted duration nonequilibrium molecular dynamics simulations of planar elongational flow,” Comput. Phys. Commun. 117, 191–199 (1999).
  * Tuckerman  _et al._ (1997) M. E. Tuckerman, C. J. Mundy, S. Balasubramanian, and M. L. Klein, “Modified nonequilibrium molecular dynamics for fluid flows with energy conservation,” J. Chem. Phys. 106, 5615–5621 (1997).
  * Todd and Daivis (2007) B. D. Todd and P. J. Daivis, “Homogeneous non-equilibrium molecular dynamics simulations of viscous flow: techniques and applications,” Mol. Simulat. 33, 189–229 (2007).
  * Legoll, Luskin, and Moeckel (2007) F. Legoll, M. Luskin, and R. Moeckel, “Non-ergodicity of the Nosé–Hoover thermostatted harmonic oscillator,” Arch. Ration. Mech. An. 184, 449–463 (2007).
  * Penna, Dzubiella, and Tarazona (2003) F. Penna, J. Dzubiella, and P. Tarazona, “Dynamic density functional study of a driven colloidal particle in polymer solutions,” Phys. Rev. E 68, 061407 (2003).
  * Dzubiella, Löwen, and Likos (2003) J. Dzubiella, H. Löwen, and C. N. Likos, “Depletion forces in nonequilibrium,” Phys. Rev. Lett. 91, 248301 (2003).
  * Krüger and Rauscher (2007) M. Krüger and M. Rauscher, “Colloid-colloid and colloid-wall interactions in driven suspensions,” J. Chem. Phys. 127, 034905 (2007).
  * Rauscher  _et al._ (2007) M. Rauscher, A. Domínguez, M. Krüger, and F. Penna, “A dynamic density functional theory for particles in a flowing solvent,” J. Chem. Phys. 127, 244906 (2007).
  * Rauscher (2010) M. Rauscher, “DDFT for Brownian particles and hydrodynamics,” J. Phys.: Condens. Matter 22, 364109 (2010).
  * Almenar and Rauscher (2011) L. Almenar and M. Rauscher, “Dynamics of colloids in confined geometries,” J. Phys.: Condens. Matter 23, 184115 (2011).
  * Gutsche  _et al._ (2008) C. Gutsche, F. Kremer, M. Krüger, M. Rauscher, R. Weeber, and J. Harting, “Colloids dragged through a polymer solution: Experiment, theory, and simulation,” J. Chem. Phys. 129, 084902 (2008).
  * Reinhardt, Scacchi, and Brader (2014) J. Reinhardt, A. Scacchi, and J. M. Brader, “Microrheology close to an equilibrium phase transition,” J. Chem. Phys. 140, 144901 (2014).
  * Yu  _et al._ (2017) H.-Y. Yu, Z. Jabeen, D. M. Eckmann, P. S. Ayyaswamy, and R. Radhakrishnan, “Microstructure of flow-driven suspension of hardspheres in cylindrical confinement: a dynamical density functional theory and monte carlo study,” Langmuir 33, 11332–11344 (2017).
  * Menzel  _et al._ (2016) A. M. Menzel, A. Saha, C. Hoell, and H. Löwen, “Dynamical density functional theory for microswimmers,” J. Chem. Phys. 144, 024115 (2016).
  * Brader and Krüger (2011) J. M. Brader and M. Krüger, “Density profiles of a colloidal liquid at a wall under shear flow,” Mol. Phys. 109, 1029–1041 (2011).
  * Scacchi, Krüger, and Brader (2016) A. Scacchi, M. Krüger, and J. M. Brader, “Driven colloidal fluids: construction of dynamical density functional theories from exactly solvable limits,” J. Phys.: Condens. Matter 28, 244023 (2016).
  * Aerov and Krüger (2015) A. A. Aerov and M. Krüger, “Theory of rheology in confinement,” Phys. Rev. E 92, 042301 (2015).
  * Scacchi, Archer, and Brader (2017) A. Scacchi, A. J. Archer, and J. M. Brader, “Dynamical density functional theory analysis of the laning instability in sheared soft matter,” Phys. Rev. E 96, 062616 (2017).
  * Stopper and Roth (2018) D. Stopper and R. Roth, “Nonequilibrium phase transitions of sheared colloidal microphases: Results from dynamical density functional theory,” Phys. Rev. E 97, 062602 (2018).
  * Goddard, Pavliotis, and Kalliadasis (2012) B. D. Goddard, G. A. Pavliotis, and S. Kalliadasis, “The overdamped limit of dynamic density functional theory: Rigorous results,” Multiscale Model. Sim. 10, 633–663 (2012).
  * Marconi and Tarazona (1999) U. M. B. Marconi and P. Tarazona, “Dynamic density functional theory of fluids,” J. Chem. Phys. 110, 8032–8044 (1999).
  * Schmidt (2022) M. Schmidt, “Power functional theory for many-body dynamics,” Rev. Mod. Phys. 94, 015007 (2022).
  * Evans (1979) R. Evans, “The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids,” Adv. Phys. 28, 143–200 (1979).
  * Archer and Evans (2004) A. J. Archer and R. Evans, “Dynamical density functional theory and its application to spinodal decomposition,” J. Chem. Phys. 121, 4246–4254 (2004).
  * Rex and Löwen (2008) M. Rex and H. Löwen, “Dynamical density functional theory with hydrodynamic interactions and colloids in unstable traps,” Phys. Rev. Lett. 101, 148302 (2008).
  * Rex and Löwen (2009) M. Rex and H. Löwen, “Dynamical density functional theory for colloidal dispersions including hydrodynamic interactions,” Eur. Phys. J. E. 28, 139–146 (2009).
  * Jeffrey and Onishi (1984) D. J. Jeffrey and Y. Onishi, “Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow,” J. Fluid Mech. 139, 261–290 (1984).
  * Happel and Brenner (2012) J. Happel and H. Brenner,  _Low Reynolds number hydrodynamics: with special applications to particulate media_ , Vol. 1 (Springer Science & Business Media, 2012).
  * Goddard, Nold, and Kalliadasis (2017) B. D. Goddard, A. Nold, and S. Kalliadasis, “2DChebClass [Software],” http://dx.doi.org/10.7488/ds/1991 (2017).
  * Sherma and Fried (2003) J. Sherma and B. Fried,  _Handbook of thin-layer chromatography_ (CRC press, 2003).
  * Rosenfeld (1989) Y. Rosenfeld, “Free-energy model for the inhomogeneous hard-sphere fluid mixture and density-functional theory of freezing,” Phys. Rev. Lett. 63, 980 (1989).
  * Faxén (1922) H. Faxén, “Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist,” Ann. Phys. 373, 89–119 (1922).
  * Wu and Li (2007) J. Wu and Z. Li, “Density-functional theory for complex fluids,” Annu. Rev. Phys. Chem. 58, 85–112 (2007).
  * Roth (2010) R. Roth, “Fundamental measure theory for hard-sphere mixtures: a review,” J. Phys.: Condens. Matter 22, 063102 (2010).
  * Tschopp and Brader (2022) S. M. Tschopp and J. M. Brader, “First-principles superadiabatic theory for the dynamics of inhomogeneous fluids,” J. Chem. Phys. 157, 234108 (2022).
  * Tschopp, Vuijk, and Brader (2023) S. M. Tschopp, H. D. Vuijk, and J. M. Brader, “Superadiabatic dynamical density functional study of brownian hard-spheres in time-dependent external potentials,” J. Chem. Phys. 158, 234904 (2023).