Dynamic Density Functional Theory with Inertia and Background Flow

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.


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 fluids 1 , inertia, with 2 and without 3 hydrodynamic interactions (HI), orientation 4,5 , and thermal fluctuations on the hydrodynamic scale 6 .Such methods predict the evolution of the onebody density distribution out of equilibrium with remarkable accuracy when compared to the underlying stochastic (Langevin) equations 7,8 .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 9 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. 10 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 11 or Kraynik-Reinelt BCs in the case of elongation flow 12 .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 systems 13,14 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-shearflow 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 14 .Additionally, Nosé-Hoover dynamics can be non-ergodic 15 .
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 16 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 17 , and separately Krüger and Rauscher 18 , 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. 19 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 20 (including the formal details of the moment closure), and there have been many practical applications including: dynamics in confined geometries 21 , advected colloids through DNA 22 , phase transitions in microrheology 23 , as well as the investigation of the hard sphere interactions in the presence of flow induced hydrodynamic lift in confined geometries 24 .The influence of solvent flow is also thought to be important in the dynamics of the density in active systems such as micro swimmers 25 .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 26 .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 27 .Additionally, there have been many extensions to DDFT specific to shear flow, including: showing the importance of HI in computing the microscopic stress tensor 28 , the investigation of laning instabilities 29 , and the emergence of nonequilibrium phase transitions depending on shear rate 30 .
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.19-21 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. 10 .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. 10 , Equation (3.1), Section 3).In particular, we consider the system of interacting stochastic differential equations (SDEs) on R 3 , which govern the positions r i and momenta p i of i = 1, . . ., N hard, spherically symmetric colloidal particles in the presence of a flowing (incompressible) background field of many more and much smaller particles given by where m is the particle mass, V is a potential energy (including one body and multi-body interactions), u is the inhomogeneous solvent flow (which may depend on both space and time), γ is the friction coefficient, k B is Boltzmann's constant, T is temperature, and the constant σ satisfies the fluctuation-dissipation relation σ = (mk ⊤ is a Gaussian white noise term with mean ⟨ζ a i (t)⟩ = 0 and autocorrelation The two conditions on the well-posedness of equations 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. 10 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. 10 , Sec 2.4.3, the stochastic equations (1a)-(1b) possess the invariant measure f (1) (r 1 , p 1 ) ∝ exp(−(2mk B T ) −1 |p 1 −mu| 2 ) only with the inclusion of an acceleration term p • ∇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. 3).In contrast, an analytical expression for stationary states of (1a)-(1b), are not known for general choice of u (other than the trivial choice 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 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 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. 5.

Potential Energy
We use the standard many-body expansion for the potential where V 1 is a one-body external potential and V n are the nbody interparticle potentials.We assume that these manybody 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, dp i /dt = 0, such that which, after substituting the corresponding value of p i into (1a), gives the overdamped Langevin dynamics 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. 19 , Rauscher 20 but here without the assumption that 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. 31 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) where f (N) (r N , p N ,t) is the probability density of finding each particle i at position r i with momentum 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. 2 and 3 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. 2 for relevant definitions.Multiplying (5) by N and integrating over all p i and all but one particle position, r 1 , yields the continuity equation ( 6) for the one body density ρ ≡ ρ (1) given by where ρ(r 1 ,t) is the one body density (Eq.( 11) below, with n = 1) and the current j is defined by Multiplying ( 5) by Np 1 /m and once again integrating over all variables except for r 1 yields the current evolution equation where dr n−1 denotes dr 2 . . .dr n .To close the system ( 6)-( 8), equations are needed for the quantities involving ρ (n) (r n ,t) and f (1) (r 1 , p 1 ,t).For the former we employ the standard adiabatic approximation 2,3,7,8,32,33 , while for the latter we use a modified local equilibrium (more properly, local steady state) approximation, which takes into account the background flow.

A. Adiabatic Approximation
We consider the free energy of the non-equilibrium system at constant temperature and introduce the Helmholtz free energy functional 7,8,34 In the above functional Λ is the de Broglie wavelength (which is superfluous) and F ex [ρ] 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 manybody parts of the potential.First we introduce the reduced phase space densities, defined by and the reduced configuration densities ρ (n) (r n ,t), given by where, again, ρ (1) ≡ ρ and dx N−n = dx n+1 . . .dx N .It is assumed that the ρ (n) (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 u ≡ 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 Under this assumption, we see that, as the momentum distribution equilibrates, the kinetic pressure term may be calculated by Hence, at steady flow, using equations ( 12) and ( 13), equation (8) becomes Now, by computing the gradient of the variational derivative of F , the Euler-Lagrange equation for the steady state density is By subtracting ( 14) from ( 15) one obtains the steady sum rule 34,35 ρ(r We assume that ( 16) holds out of equilibrium to complete the adiabatic approximation on the higher n-body densities.

B. Moment Closure and Local Equilibrium Approximation
Using the definition of F along with the (out of equilibrium) approximation (16), equation (8) becomes where Note that at steady state, 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 v(r 1 ,t).We write f (1) as the sum of local steady and non-steady parts where The first few moments of f (1) and we impose the integral restrictions By using the integral conditions ( 22) and ( 25), combined with the expansion (19), the continuity equation ( 6) becomes We may also compute A(r 1 ,t) explicitly as We insert this expression for A(r 1 ,t) into ( 17) to obtain Now consider the following identity, for an arbitrary, suitably well behaved vector field z(r 1 ,t) This result follows from the product rule of differentiation and, e.g., equation (27) of Ref. 3. By setting z = v(r 1 ,t), using the continuity equation (27), and neglecting the f ns term, (since it can be shown to be small, at least in the overdamped limit 31 γ → ∞), we obtain our equations of motion.
where w(r 1 ,t) := v(r 1 ,t) − u(r 1 ,t).The momentum equation in (31) should be compared with the inertial DDFT given in Eq. ( 30) of Ref. 3. 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 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. 3 to the moving reference frame, which is equivalent to shifting the transport of the velocity distribution v(r 1 ,t) by the external flow u(r 1 ,t).We see that the velocity advection transforms to which differs from the advection term in (31).We note that it is possible, via a careful choice of u(r 1 ,t), to make the steady state velocity and density distribution resulting from (31) and a frame-transformed Eq. ( 30) in Ref. 3 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 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.

C. 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 where A ls is given by (28) with f ns = 0. Differentiating (6)  we obtain Additionally, by taking the divergence of (33), we have Combining Eqs. ( 34) and ( 35), we obtain It remains to treat the term on the right hand side.We note that in Ref. 2 we have A ls = A le , and, by the results of Ref. 31, at least for a two-body potential truncation of F ex [ρ], terms containing A le are negligible in the high-friction limit, γ → ∞.
Hence, under a similar adiabatic elimination process, based on a multiscale analysis of the Fokker-Planck equation ( 5), we expect that A ls will also be negligible in the presence of the solvent flow u.Therefore, (36) maybe approximated as The second-order time derivative ∂ 2 t ρ 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., Krüger and Rauscher 18 , Rauscher et al. 19 , Rauscher 20 , Almenar and Rauscher 21 , are recovered.Assuming ρ is strictly positive, we may also establish that the free energy is a monotonically decreasing function of time along the solvent flow where we have assumed that the boundary terms vanish; this is the case, e.g., when ρ = 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) 2,36,37 .By introducing the friction tensor Γ ∈ R 3N×3N , the nonequilibrium Langevin dynamics (1a)-(1b) becomes where B = (mk B T Γ) 1/2 and Γ comprises N 2 positive definite resistance matrices Γ i j given by Γ i j = γ1 + γ Γi j , where 1 is the 3 × 3 identity matrix, and Γi j ∈ R 3×3 .Physically we typically regard the forces introduced by γ Γi j as short range HIs mediated by the bath, (also known as lubrication forces) and γ 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, Z 1 = Z 2 = 0 (defined below) and γ → ∞ respectively.After carefully repeating the derivation from Section III onwards, we use the Enskog approximation to resolve f (2) [see Eq. (10)] where g is the pair correlation function, which must be supplied by auxiliary means.Additionally we restrict the HIs to be two-body, thus where Z 1 and 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 38 , Happel and Brenner 39 ).The corresponding DDFT is then 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 bath 2 .In addition, setting the HI tensors to zero (Z 1 = Z 2 = 0) recovers the results of Ref. 3.

A. Overdamped Limit
We now consider the overdamped limit with the inclusion of HI.Assuming that inter-particle HI are weak, that is Z 2 ≈ 0, we obtain where D is the diffusion tensor, known rigorously to be positive definite (see, e.g., Ref. 31) given by Now by adding and subtracting the same term to the continuity equation and combining (43) and ( 45) we obtain which may be compared with the DDFT in Eq. ( 27) of Ref. 19.Equation ( 46) is the Smoluchowski equation, differing from the one found in Rauscher et al. 19 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 31 .

V. RESULTS FROM NUMERICAL SOLUTION
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 F ex [ρ] given by fundamental measure theory (FMT), and compare the resulting dynamics to previously derived DDFTs, specifically the inertial DDFT of Ref. 3 and the overdamped with external flow theory of Ref. 19.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 F ex [ρ].For example, we have also tested our DDFT in conjunction with a simple mean-field approximation for F ex [ρ] 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 MAT-LAB 2DChebClass package 40 .
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 41 , whereas the addition of a local sink models the flow of colloids around a plughole-like exit.The third example concerns flow in a confined geometry (a periodic two-dimensional capillary) for a simple model system in haemodynamics 24 .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. 19 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 ς (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 ρ(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 loidal particles.This hard sphere exclusion is described by the excess free energy F ex [ρ]; here we use the three-dimensional FMT functional of Rosenfeld 42 , 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 43 type corrections to the solvent fluid (c.f., Rauscher et al. 19 , Rauscher 20 ) 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 21 .
We remark now on initial conditions: In each of the numerical examples presented we assume u = 0 at t = 0, but thereafter (t > 0) 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 DFT 34,44 problem (δ F [ρ]/δ ρ) = 0. We obtain the initial condition via Picard iteration, a standard technique within the DFT community, see e.g., Roth 45 , 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.

A. Line Source in a Uniform Flow
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 e x and 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 u(x, y) = (u(x, y), v(x, y)) ⊤ where u(x, y) and v(x, y) are the scalar velocity fields of the solvent in the e x and e y directions respectively, given by where U and Q are the magnitudes of the uniform and line flow sources, respectively.We note that u(x, y) is derived from the gradient of the potential φ (x, y) = Ux + Q/(2π) log 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π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) = λ (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 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 λ = 0.01, σ = 50, τ = 0.01 and γ = 2 for N = 50 hard spheres.In Figure 1 we show the transient densities ρ(r,t) and velocities v(r,t) obtained from (31), alongside the solutions obtained with external flow turned off (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πU), 0) ≈ (−0.56, 0), which generates a dip in the density profile in the vicinity of that point.

B. Flow in a Localised Vortex and Sink
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 where R = 0.9 is a scalar parameter controlling the strength and direction of the vortex, Q = −1.57is 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 → ±∞, so we include an additional multiplicative term to decay the solvent velocity to zero such that ũ(x, y) = u(x, y)e −(x 2 +y 2 )/α , with the choice α = 8.5.We have maintained the same external potential width σ , time scale τ, friction coefficient γ, as well as the number of particles N = 50 from the previous example, and once again we have used Rosenfeld 42 FMT as the excess over ideal gas free energy functional F ex .For this example we use a slowly decay- ing confining potential with λ = 0.001 and (x a , y a ) = (0, 4).Here, V 1 (x, y,t) → λ (x 2 + y 2 ) as t → ∞, 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 u ̸ = 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.

C. 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 ∈ [−2, 2]) and length 2π (between y ∈ [0, 2π]).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, j • e x = 0 on x = −2, x = 2 and t ≥ 0, and ρ(r,t) = ρ(r + 2πe y ,t) for (x, y) ∈ [−2, 2] × [0, 2π] and t ≥ 0. We impose a parabolic external flow of the form where W = 4 is the width of the channel, and A(t) controls the strength and direction of the flow via: where α = 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 ρ = Z(cos(y + y 0 ) + ρ 0 ), where y 0 = π and ρ 0 = 2.We choose Z such that ρ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 timedependent, 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 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 ρ(r,t) and v(r,t); (iii) the density of the particle systems is relatively low.For the latter, in the overdamped limit, Almenar and Rauscher 21 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(r 1 , r 2 ) obtained from the the Fokker-Planck equation for the twopoint 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 46,47 , or by using the methods of power functional theory 33 .
For the remainder of the derivation, we rely on the usual approximations: (i) The non-steady contributions from the onebody distribution which are not captured by the local-steady approximation are small or may be neglected.(ii) For twobody 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. 3, and by neglecting external flows only we recover the DDFT derived in Ref. 2.
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.external potentials," J. Chem.Phys.158, 234904 (2023).

FIG. 1 :
FIG. 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 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 u = 0, and the four on the right are for a line source, i.e., uniform flow with u given by Eqs.(47a)-(47b).

FIG. 2 :
FIG. 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 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 u = 0, and the four on the right are for the case where u given by (51a)-(51b).

5 FIG. 3 :
FIG. 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 (ρ o − ρ 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.