Custom flow in overdamped Brownian Dynamics

TL;DR

This paper introduces a general iterative simulation scheme to determine the unique external force field needed to produce a specified flow in overdamped Brownian systems, facilitating systematic analysis of complex colloidal flows.

Contribution

It presents a novel computational method based on the force balance equation that can tailor flow fields and densities, extending the power functional framework for non-equilibrium systems.

Findings

The method converges to a unique external force field for prescribed flows.

It can generate conservative forces for target density profiles in equilibrium.

The approach provides detailed physical insights into complex flow behaviors.

Abstract

When an external field drives a colloidal system out of equilibrium, the ensuing colloidal response can be very complex and obtaining a detailed physical understanding often requires case-by-case considerations. In order to facilitate systematic analysis, here we present a general iterative scheme for the determination of the unique external force field that yields a prescribed inhomogeneous stationary or time-dependent flow in an overdamped Brownian many-body system. The computer simulation method is based on the exact one-body force balance equation and allows to specifically tailor both gradient and rotational velocity contributions, as well as to freely control the one-body density distribution. Hence compressibility of the flow field can be fully adjusted. The practical convergence to a unique external force field demonstrates the existence of a functional map from both velocity and density to external force field, as predicted by the power functional variational framework. In equilibrium, the method allows to find the conservative force field that generates a prescribed target density profile, and hence implements the Mermin-Evans classical density functional map from density distribution to external potential. The conceptual tools developed here enable one to gain detailed physical insight into complex flow behaviour, as we demonstrate in prototypical situations.

Figures

Equations

$$\displaystyle\rho({\mathbf{r}},t)=\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\Big{\rangle},$$

$$⟨\cdot ⟩$$

$$\gamma {\mathbf{v}}_i$$

$$\mathbf{J}(\mathbf{r},t)$$

$$\frac{\partial \rho (\mathbf{r},t)}{\partial t}$$

$$\displaystyle{\bf F}_{\rm int}({\mathbf{r}},t)=-\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla_{i}u({\mathbf{r}}^{N})\Big{\rangle}.$$

$$\gamma \mathbf{J}=-k_{B}T\nabla \rho +{\mathbf{F}}_{\mathrm{int}}+\rho {\mathbf{f}}_{\mathrm{ext}},$$

$$\mathbf{v}$$

$$\gamma \mathbf{v}$$

$${\mathbf{f}}_{\text{ext}}=k_{B}T\nabla \ln \rho -{\mathbf{f}}_{\text{int}}+\gamma \mathbf{v},$$

$${\mathbf{f}}_{\text{ext}}^{(k)}=k_{B}T\nabla \ln \rho -{\mathbf{f}}_{\text{int}}^{(k-1)}+\gamma \mathbf{v},$$

$${\mathbf{f}}_{\text{ext}}^{(0)}$$

$$\displaystyle\Delta\rho=\int d{\mathbf{r}}\big{(}\rho({\mathbf{r}})-\rho^{(k)}({\mathbf{r}})\big{)}^{2}/V<\epsilon_{1},$$

$$\displaystyle\Delta J=\int d{\mathbf{r}}\big{\lvert}{\bf J}({\mathbf{r}})-{\bf J}^{(k)}({\mathbf{r}})\big{\rvert}^{2}/V<\epsilon_{2},$$

$${\mathbf{f}}_{\mathrm{int}}([\rho ,\mathbf{v}],\mathbf{r},t)$$

$${\mathbf{f}}_{\mathrm{int}}(\mathbf{r})$$

$${\mathbf{f}}_{\mathrm{ext}}$$

$${\mathbf{f}}_{\mathrm{ext}}^{(k)}(\mathbf{r})$$

$${\mathbf{f}}_{\mathrm{ext}}^{(0)}(\mathbf{r})$$

$$v_{\text{ext}}^{(k)}(x)=-k_{B}T\ln \rho (x)+\int dx{f}_{\text{int}}^{(k-1)}(x).$$

$${\mathbf{f}}_{\mathrm{int}}([\rho ],\mathbf{r})$$

$${\mathbf{f}}_{\mathrm{ext}}$$

$${\mathbf{f}}_{\mathrm{int}}^{(k-1)}=k_{B}T\nabla \ln ({\rho }^{(k-1)})-{\mathbf{f}}_{\mathrm{ext}}^{(k-1)},$$

$${\mathbf{f}}_{\mathrm{ext}}^{(k)}-{\mathbf{f}}_{\mathrm{ext}}^{(k-1)}=k_{B}T\nabla \ln (\rho /{\rho }^{(k-1)}),$$

$$\ln (\rho /{\rho }^{(k-1)})=-({\rho }^{(k-1)}-\rho )/\rho +\dots$$

$${\mathbf{f}}_{\mathrm{ext}}^{(k)}(\mathbf{r},t_{\text{cg}})$$

$${\mathbf{f}}_{\text{ext}}^{(0)}(\mathbf{r},t_{\text{cg}})={\mathbf{f}}_{\text{ext}}(\mathbf{r},t_{\text{cg}}^{\prime }),$$

$$\phi(r)=\left\{\begin{array}[]{ll}4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}+\frac{1}{4}\right]&\text{if }r<r_{c}\\ 0&\text{otherwise,}\\ \end{array}\right.$$

$${\mathbf{r}}_i(t+\Delta t)={\mathbf{r}}_i(t)+\frac{\Delta t}{\gamma }[-{\nabla }_{i}u({\mathbf{r}}^N)+{\mathbf{f}}_{\mathrm{ext}}({\mathbf{r}}_i,t)]+{\mathbf{\eta }}_i(t),$$

$$\rho (x)=c_1{\sin }^2(\pi x/L)+c_2,$$

Full text

I Introduction

The controlled application of an external field is a powerful means to drive colloidal systems of mesoscopic suspended particles out of equilibrium Löwen (2008); Erbe et al. (2008). The resulting complex interplay of the equilibrium properties of the system with the external perturbation is already present in Perrin’s pioneering work on colloidal sedimentation Perrin (1909). Gravitationally driven colloidal suspensions Köppl et al. (2006); Simeonova and Kegel (2004); Palacci et al. (2010) remain to this day primary model systems for studying structure formation phenomena. There is a large spectrum of different types of further specific external influence on colloids, such as the response of charged colloids to external electric fields Velev and Bhatt (2006); O’Brien and White (1978), and magnetic field-induced transport of both diamagnetic and paramagnetic colloidal particles Tierno et al. (2007); Snezhko and Aranson (2011) across a substrate, where the colloidal motion was recently analysed in terms of the powerful concept of topological protection against perturbation Loehr et al. (2016). The external magnetic fields in these setups varied periodically in both space and time. Alternatively, exerting shear-like flow on colloidal dispersions provides in-depth insights into important effects in fundamental material science, such as shear thickening Bender and Wagner (1996); Wagner and Brady (2009) and glass formation Besseling et al. (2007). Furthermore, optical tweezers form powerful and flexible tools for the generation of complex colloidal flow patterns Ruffner and Grier (2012); Sun et al. (2009); Lutz et al. (2004).

In order to systematically study the response of a soft material to an external perturbation, one typically first fixes the external field and then studies the resulting colloidal motion under the effect of the field. Certainly, this concept is compatible with our understanding of a causal relationship between forces and the motion that they generate. If the external force field is static and conservative, then the system will in general reach a new equilibrium state. This response of a complex system to such an external influence might be highly non-trivial. As a result of the action of the external potential, the system will in general become spatially inhomogeneous. In seminal work, Mermin showed for quantum systems Mermin (1965), as Evans subsequently did for the classical case Evans (1979), that a functional inversion of the relationship between external potential and one-body density distribution applies. Hence, reversing the above “causal” relationship, a unique mathematical map exists from the one-body density distribution to the corresponding external potential. This is an important and fundamental result of modern Statistical Physics, which generalizes Hohenberg and Kohn’s earlier result for quantum ground states Hohenberg and Kohn (1964). The functional relationship forms the basis of classical density functional theory, which is used in practically all modern microscopic theoretical treatments of spatially inhomogeneous systems Evans et al. (2016). Once the external potential is specified, the Hamiltonian is known (the internal interactions remain unchanged) and hence all equilibrium properties of the system are determined, and become functional dependent on the density distribution.

In this work we address the functional map in non-equilibrium steady and time-dependent states in overdamped Brownian many-body systems. We first set the desired colloidal motion, as specified by both the velocity field (or, equivalently, the one-body current distribution) and the density profile, and then determine the specific external force field that creates the prescribed motion in steady state. We develop and validate a numerical iterative method that enables efficient and straightforward implementation of this task.

That the map from motion to external force field exists and that it is unique follows formally from the power functional variational principle in general time-dependent non-equilibrium Schmidt and Brader (2013). The functional relationship has not, however, been explicitly demonstrated in an actual many-body system. Here we provide the first such demonstration, both for steady states and for time-dependent non-equilibrium.

As a special case, we also apply the inversion method to equilibrium systems. Here it allows to find the specific conservative external force field that stabilizes a predefined density distribution. As flow is absent and kinetic energy contributions are trivial in equilibrium, this method also applies to inertial (i.e. Hamiltonian) systems. An alternative iterative numerical method that also implements the functional map in equilibrium is presented in Ref. Fortini et al. (2014).

The conceptual progress in demonstrating the non-equilibrium inversion map explicitly enables the practical solution to the problem of generating tailor-made flow in complex systems. In the method that we present, the sole requirement is, besides the ability to freely control the external force field, to be able to measure the internal one-body force density distribution. This is a readily available quantity in many-body simulations, and it is conceivably also accessible in experimental work. We envisage that tailoring freely flow on demand constitutes a powerful concept for the systematic study of non-equilibrium physics. Here we investigate as a concrete model problem, the task of collecting particles in a certain region of space via a potential trap. The system is initially a homogeneous fluid, and we (i) speed up the natural dynamics by a factor of 2 and (ii) demonstrate that any unwanted effects due to superimposed external flow (e.g. due to convection or sedimentation in the system) can be fully cancelled.

The paper is organized as follows. In Sec. II we describe how we obtain the iterative method, based on the exact force balance relationship in overdamped Brownian dynamic, covering both non-equilibrium steady states and equilibrium. In Sec. III we present results for several model situations in which we custom-tailor the flow in a many-body system of repulsive particles. Sec. IV contains a discussion and provides some conclusions. Details about particle current sampling and convergence properties are given in the appendix.

II Theory

II.1 Dynamical one-body force balance

We consider a system of $N$ interacting Brownian particles in the overdamped limit, where inertial effects are absent and we neglect hydrodynamic interactions. The one-body density distribution (“density profile”) at space point ${\mathbf{r}}$ and time $t$ is given by

[TABLE]

where the angles denote a statistical average, the expression inside the angles is the microscopic density operator, with $\delta(\cdot)$ indicating the Dirac distribution, and ${\mathbf{r}}_{i}$ denoting the position of particle $i=1\ldots N$. In the Fokker-Planck picture, the information required for carrying out the average is encoded in the many-body probability distribution $\Psi({\mathbf{r}}^{N},t)$ of finding microstate ${\mathbf{r}}^{N}\equiv{\mathbf{r}}_{1}\ldots{\mathbf{r}}_{N}$ at time $t$. The average is then defined as

[TABLE]

where the integral runs over configuration space, i.e. each ${\mathbf{r}}_{i}$ is integrated over the system volume. The time evolution of $\Psi$ is determined by the Smoluchowski equation $\partial\Psi/\partial t=-\sum_{i}\nabla_{i}\cdot{\bf v}_{i}\Psi$. Here the configuration space velocity ${\bf v}_{i}$ of particle $i$ is given on the many-body level by the instantaneous relation

[TABLE]

where $\gamma$ is the friction constant against the implicit solvent, $k_{B}$ is the Boltzmann constant, $T$ denotes absolute temperature, $\nabla_{i}$ is the partial derivative with respect to ${\mathbf{r}}_{i}$, the interparticle interaction potential is denoted by $u({\mathbf{r}}^{N})$, and ${\bf f}_{\rm ext}({\mathbf{r}},t)$ is the external force field, which in general is position- and time-dependent. The three contributions on the right hand side of (3) correspond to thermal diffusion (first term), deterministic motion due to interparticle interactions (second term), and the externally imposed force (third term). This formulation of the dynamics is analogous to the Langevin picture, where instead of (2), the average is taken over a set of stochastic particle trajectories for which a random (position) noise provides the effects of thermal motion. The corresponding discretized version is readily implementable in Brownian Dynamics (BD) computer simulations (details of our implementation are given in Sec III). Note that the configuration space velocity ${\bf v}_{i}$ defined in (3) is different from the average over the fluctuating velocity over realization of the noise, as represented in BD simulations.

We next supplement (1) by a corresponding one-body current distribution, defined as

[TABLE]

where ${\bf v}_{i}$, at time $t$, is given via (3). We show in detail in appendix A how a forward-backward symmetrical time derivative of the particle positions can be used in BD to represent ${\bf v}_{i}$ in (4).

The density distribution (1) and the current profile (4) are linked via the continuity equation,

[TABLE]

The many-body coupling in (3) arises due to the presence of the internal interaction potential $u({\mathbf{r}}^{N})$. On the one-body level, it is hence natural to define a corresponding internal force density field via

[TABLE]

Here contributions to the average occur due to two effects, namely (i) due to the bare value of the internal force field $-\nabla_{i}u({\mathbf{r}}^{N})$, but also (ii) due to the probability of finding particle $i$ at the considered space point ${\mathbf{r}}$, as measured by the delta function. Averages such as (6) hence constitute microscopically resolved force densities.

By multiplying (3) with $\delta({\mathbf{r}}-{\mathbf{r}}_{i})$, summing over $i$, averaging according to (2), and identifying the one-body fields (1), (4), and (6), it is straightforward to show that

[TABLE]

which we use as a basis for the non-equilibrium inversion procedure. Defining the microscopic velocity profile ${\bf v}({\mathbf{r}},t)$ simply as the ratio between current profile and density profile,

[TABLE]

allows us to rewrite (7), after division by the density profile $\rho$, as

[TABLE]

Here the internal force field ${\bf f}_{\rm int}$ is defined as the internal force density (6) normalized with the density profile, i.e. ${\bf f}_{\rm int}({\mathbf{r}},t)\equiv{\bf F}_{\rm int}({\mathbf{r}},t)/\rho({\mathbf{r}},t)$.

BD computer simulations allow straightforward access to the individual contributions to the force balance relationship (9). Sampling the density profile is straightforward either using the traditional counting method or more advanced techniques Borgis et al. (2013); de las Heras and Schmidt (2018). The ideal (diffusive) force field $-k_{B}T\nabla\ln\rho$ is readily obtained by (numerical) differentiation of the density profile. The internal force density field ${\bf F}_{\text{int}}$ can be sampled as an average over BD realizations of the time evolution, or as an average over time when investigating steady states; note that in BD one has direct access to the internal force on the many-body level, $-\nabla_{i}u({\mathbf{r}}^{N})$ in (6).

In typical applications, the external force field ${\bf f}_{\rm ext}({\mathbf{r}},t)$ is prescribed and $\rho({\mathbf{r}},t)$ and ${\bf v}({\mathbf{r}},t)$ emerge as a result of the coupled many-body dynamics. In the following, we address the inverse problem of prescribing $\rho$ and ${\bf v}$ a priory and calculating the required form of ${\bf f}_{\rm ext}$ that makes these fields stationary, such that the prescribed field values are identical to the true dynamical averages of density (1) and velocity (4), (8).

II.2 Inversion in nonequilibrium steady states

Let $\rho({\mathbf{r}})$ and ${\bf v}({\mathbf{r}})$ be the predefined stationary (i.e. time-independent) “target” profiles for density and velocity. In order to represent a valid steady state, the resulting target current profile $\rho{\bf v}$, cf. (8), must be divergence-free, $\nabla\cdot\rho{\bf v}=0$, which follows from (5) and represents a necessary condition on the allowed set of target functions $\rho$, ${\bf v}$. The external force field ${\bf f}_{\rm ext}({\mathbf{r}})$ that makes the target profiles stationary is obtained from first rearranging (9) as

[TABLE]

where the internal force field, ${\bf f}_{\text{int}}({\mathbf{r}})={\bf F}_{\text{int}}({\mathbf{r}})/\rho({\mathbf{r}})$, is the only unknown quantity on the right hand side, as $\rho$ and ${\bf v}$ are known input quantities. Here ${\bf F}_{\rm int}({\mathbf{r}})$ is from (6) in steady state.

In order to determine ${\bf f}_{\rm int}$, and hence ${\bf f}_{\rm ext}$ via (10), we proceed in two steps. First we present a fixed-point iterative scheme to solve (10), in which the $k$-th iteration step is defined via

[TABLE]

where the targets, $\rho({\mathbf{r}})$ and ${\bf v}({\mathbf{r}})$, are kept fixed for all steps $k$. Here ${\bf f}_{\rm int}^{(k-1)}={\bf F}_{\rm int}^{(k-1)}/\rho$, where ${\bf F}_{\rm int}^{(k-1)}$ is the internal force density sampled in the previous iteration step, $k-1$. Data for ${\bf F}_{\rm int}^{(k-1)}$ was obtained from BD sampling under the prescribed external force field ${\bf f}_{\text{ext}}^{(k-1)}$. In order to initialize the iteration scheme, we set the external force field at step $k=0$ simply as

[TABLE]

which is the exact external force field for the case of an ideal gas. Prescribing (12) allows to sample ${\bf F}_{\rm int}^{(0)}$ in BD, and then use ${\bf f}_{\rm int}^{(0)}$ as the required input for iteration step $k=1$ in (11). This completes the description the functional inversion.

At each iteration step we also sample both the one-body density and one-body current profiles, $\rho^{(k)}$ and ${\bf J}^{(k)}$; details on how to sample the current in BD are provided in Appendix A. As a criterion for judging the eventual convergence of $f_{\rm ext}^{(k)}$ to the real external force field that makes the target density and current profiles stationary, i.e. the solution of (10), we use the difference between the target and the sampled profiles at step $k$, i.e.,

[TABLE]

where $V=\int d{\mathbf{r}}$ is the system volume, and $\epsilon_{1},\epsilon_{2}>0$ are small tolerance parameters. Numerical details of our implementation are given in Appendix B.

We find that in practice the iteration method converges reliably in all cases considered; results are shown below in Sec. III. That a solution exists for ${\bf f}_{\rm ext}$ and that it is unique are nontrivial properties of our scheme. We expect existence and uniqueness to hold, however, based on the power functional variational framework Schmidt and Brader (2013), which is a novel approach for the statistical description of the dynamics of many-body systems. The central object of power functional theory (PFT) is a “free power” functional $R_{t}[\rho,{\bf J}]$ of the one-body density and current or analogously, viz. (8), of the density and the velocity field. $R_{t}$ has units of energy per time (power) and plays a role analogous to the free energy functional (as detailed below) in equilibrium. It consists of an ideal gas contribution ($W_{t}^{\rm id}$), an excess (over ideal gas) part due to the internal interactions ($W_{t}^{\rm exc}$) and an external power ($X_{t}$) contribution, according to $R_{t}=W_{t}^{\rm id}+W_{t}^{\rm exc}-X_{t}$.

PFT implies that ${\bf f}_{\rm int}$ is a unique functional of density and current distributions, or equivalently of density and velocity profile. In particular, ${\bf f}_{\rm int}$ can be expressed as a functional derivative of the intrinsic excess (over ideal gas) free power functional $W_{t}^{\rm exc}[\rho,{\bf J}]$, as

[TABLE]

where the density distribution $\rho$ is kept fixed upon the variation, at fixed time $t$.

Typically, one would split further into adiabatic and superadiabatic contributions, $W_{t}^{\rm exc}=\dot{F}_{\rm exc}[\rho]+P_{t}^{\rm exc}[\rho,{\bf J}]$, where $\dot{F}_{\rm exc}[\rho]$ is the time derivative of the equilibrium excess (over ideal gas) Helmholtz free energy functional, and $P_{t}^{\rm exc}[\rho,{\bf J}]$ is the superadiabatic contribution, which describes the genuine nonequilibrium effects. This splitting offers great advantages in terms of the classification of the different types of forces that occur, but it is not required for our present purposes. We rather work directly with ${\bf f}_{\rm int}$. Recall that this is directly accessible via Eq. (6) in BD simulations.

The fact that ${\bf f}_{\rm int}$ is generated from a current-density functional, via functional differentiation (15), implies that the force field itself is a functional of density and current (or velocity profile). Hence (15) constitutes a unique map from density and velocity to the internal force field,

[TABLE]

where the right hand side is the force field functional (15) evaluated at the target profiles $\rho$ and ${\bf v}$; the left hand side is the corresponding (hitherto unknown) specific form of the internal force on the left hand side of (10). Hence by inserting (15) into (10) we obtain the explicit form

[TABLE]

with the iteration procedure (11) and (12) being a practical scheme for evaluating the right hand side. Note that (17) is an explicit expression for ${\bf f}_{\rm ext}$; no hidden dependence on ${\bf f}_{\rm ext}$ occurs on the right hand side. Recall that from (15), the internal force field depends solely on the “kinematic” fields $\rho$ and ${\bf v}$, but not on the external force field. This completes the proof. Before presenting results, we revisit the equilibrium case.

II.3 Inversion in equilibrium

In equilibrium, the velocity profile is identically zero, and we therefore can simply set the target ${\bf v}({\mathbf{r}})\equiv 0$ and prescribe $\rho({\mathbf{r}})$ in order to find the corresponding external force field ${\bf f}_{\rm ext}({\mathbf{r}})$. The external force field is necessary of conservative, gradient form, ${\bf f}_{\rm ext}({\mathbf{r}})=-\nabla v_{\rm ext}({\mathbf{r}})$, where $v_{\rm ext}$ is the external potential energy. Clearly, there is no dependence on time, and, as we show, one only needs to carry out equilibrium averages. Hence the method also applies to Hamiltonian systems, as the kinetic contributions are trivial.

In equilibrium we can simplify (10) to obtain

[TABLE]

which constitutes an explicit expression for the specific external force field that generates the given target profile $\rho$ in equilibrium.

As it is the case in nonequilibrium steady state, the internal force is unknown, but it can be found iteratively. The iteration step is

[TABLE]

and the external force is initialized with the exact solution of an ideal gas,

[TABLE]

We then sample ${\bf f}_{\text{int}}^{(0)}$ in equilibrium, under the external force ${\bf f}_{\rm ext}^{(0)}$, and then iterate, on the basis of (19), until the difference between the target and the sampled density profiles is small, cf. (13). As only the internal force and the density profiles are required, the sampling can be performed using BD or molecular dynamics simulations. If one wishes to use the Monte Carlo method, then one needs the actual value of the potential $v_{\text{ext}}^{(k)}$ instead of the force, ${\bf f}_{\rm ext}^{(k)}=-\nabla v_{\text{ext}}^{(k)}$. In systems that effectively depend on only one coordinate, say $x$, the potential at each iteration can be easily obtained from the internal force profile, by performing a one-dimensional spatial integral

[TABLE]

In two- and three-dimensional systems a line integral or, more generally, the use of an inverse operator $\nabla^{-1}$ is required to obtain the potential from the force field. Hence the situation is similar to the one addressed in modern “force sampling” methods that yield the density profile Borgis et al. (2013); de las Heras and Schmidt (2018).

That the method converges to a unique external potential is guaranteed by the Mermin-Evans functional map from density profile to external potential Mermin (1965); Evans (1979). In particular, the internal force field is obtained as a functional derivative of the excess free energy functional via

[TABLE]

Inserting this into (18) yields

[TABLE]

which is an explicit expression for the external force field, given $\rho$ as an input, in analogy to the nonequilibrium case, (15) and (17). For completeness, and briefly, (23) is formally obtained from the more general (15) and (17) by observing that in equilibrium $\delta W_{t}^{\rm exc}/\delta{\bf J}=\delta\dot{F}_{\rm exc}/\delta{\bf J}$, where $\dot{F}_{\rm exc}=\int d{\mathbf{r}}{\bf J}\cdot\nabla\delta F_{\rm exc}/\delta\rho$, cf. Schmidt and Brader (2013).

We next clarify the relationship to the method of Ref. Fortini et al. (2014). Note that at any step in the iteration, given the external force field ${\bf f}^{(k-1)}_{\rm ext}$, the internal force field may be written as

[TABLE]

in terms of the density distribution $\rho^{(k-1)({\mathbf{r}})}$ at equilibrium with that external force. Then Eq. (19) may be written as the change in the force, along the iterative procedure,

[TABLE]

that vanishes when the target density is achieved, $\rho=\rho^{(k-1)}$. The integration of (25), to get the change in external potential, and the linear expansion

[TABLE]

give precisely the method used in Ref. Fortini et al. (2014).

II.4 Inversion in time-dependent nonequilibrium

In order to perform the inversion in time-dependent nonequilibrium, we carry out the procedure of Sec. II.2 at a discretized sequence of (coarse-graining) times $t_{\text{cg}}$ during the time evolution. The method propagates the system forward in time, in sync with the target time evolution. At each coarse-graining time step the required external force field is obtained (via iteration) such that the prescribed target density $\rho({\mathbf{r}},t_{\text{cg}})$ and velocity field ${\bf v}({\mathbf{r}},t_{\text{cg}})$ are identical to their respective values in the target time evolution of the system. We interpolate linearly the values for the external force field between two consecutive times, which we find sufficient for the test cases presented below.

In detail, at each coarse-graining time step $t_{\text{cg}}$ we iterate the value of the external field according to

[TABLE]

where $\rho({\mathbf{r}},t)$ and ${\bf v}({\mathbf{r}},t)$ are the target fields, which enter via their values at time $t_{\text{cg}}$. The time $t_{\text{cg}}$ is kept fixed under the iteration $k\to k+1$ described by (27). For the first time step we initialize the external force using the exact ideal gas solution, Eq. (12). For the subsequent time steps we initialize the iterative scheme using the solution of the previous time step:

[TABLE]

where $t^{\prime}_{\text{cg}}$ indicates the time step previous to $t_{\text{cg}}$.

The iteration method proceeds forward in time. In order to correctly account for memory effects, the many-body dynamics evolves according to continuous, valid trajectories over the entire time-dependent dynamics. In the BD simulations, this requires to start a new coarse-graining time step using the many-body configuration(s) obtained at the end of previous coarse-graining time step. At the end of the process, the entire field ${\bf f}_{\rm ext}({\mathbf{r}},t)$ is known, and as consistency check, can be input into a “bare” non-steady BD run, in order to validate that the targets $\rho({\mathbf{r}},t)$ and ${\bf v}({\mathbf{r}},t)$ are met during the entire course of time.

III Results

In the following we demonstrate that the straightforward application of the method allows to cast new light on fundamental physical effects by studying a two-dimensional model fluid of Brownian particles interacting via the common Weeks-Chandler-Anderson potential Weeks et al. (1971), i.e. a purely repulsive, truncated-and-shifted Lennard-Jones (LJ) pair potential Weeks et al. (1971),

[TABLE]

where the parameters $\epsilon$ and $\sigma$ set the energy and length scales, respectively, and $r$ indicates the center-center distance of the particle pair. The cutoff distance $r_{c}=2^{1/6}\sigma$ is located at the minimum of the standard LJ potential, and hence the interaction is purely repulsive.

The particles are in a square box of length $L$ with periodic boundary conditions along both directions. Using the standard Euler algorithm, the Langevin equation of motion is integrated in time via

[TABLE]

where $\boldsymbol{\eta}_{i}$ is a delta-correlated Gaussian random displacement with standard deviation $\sqrt{2\Delta tk_{B}T/\gamma}$ in accordance with the fluctuation-dissipation theorem. Here, $\Delta t$ is the integration time step that we set to $\Delta t/\tau=10^{-4}$ with $\tau=\sigma^{2}\gamma/\epsilon$ the unit of time; the friction constant is set to $\gamma=1$.

III.1 Effective one-dimensional system

A considerably large class of nonequilibrium situations is effectively of one-dimensional nature, where both the density profile and the current distribution depend only on a single coordinate, say $x$, and the flow direction is along the $x$-axis (i.e. with no shear motion occurring). Then the steady state condition reduces to the requirement of the current being constant, ${\bf J}(x)=J_{0}{\bf e}_{x}$, where $J_{0}=\rm const$ and ${\bf e}_{x}$ is the unit vector in the $x$-direction. Hence from (8) the velocity and the density profile possess a reciprocal relationship: $v(x)=J_{0}/\rho(x)$.

We first study such an effective one-dimensional problem with $N=30$ particles in a square simulation box of side length $L/\sigma=10$. We choose the target density profile to contain a single nontrivial Fourier component, that modulates the homogeneous fluid,

[TABLE]

with $c_{1}\sigma^{2}=0.12$ and $c_{2}\sigma^{2}=0.24$ such that $\int dr\rho({\mathbf{r}})=N$. See the density profile in Fig. 1(a). The temperature is set to $k_{B}T/\epsilon=1$.

Our inversion method facilitates the study of fundamental aspects of driven systems. As an illustrative example, we construct a family of steady states that share the same density profile, cf. Eq. (31), but possess different values $J_{0}$ of the (constant) target current. In Appendix B we describe numerical details of the concrete implementation of the iterative procedure.

We show in Fig. 1 the external force field required to make the target density profile (31) stationary, for a range of different values of $J_{0}$. In Figs 1(a) and (b) the final converged density and current profiles are shown; these are indeed (numerically) identical to their targets. We consider four steady states with values of the current $J_{0}\sigma\tau=0$ (equilibrium), $0.1,0.5$ and $1$. The specific external force field required to produce each such steady states is depicted in Fig. 1(c) for all four cases. The force fields can be represented as the sum of a spatially constant force offset plus a conservative potential contribution. The constant offset drives the particle flow and it can be calculated as the spatial average of the total external force field. The conservative term generates the density modulation. As expected, in the equilibrium case ($J_{0}=0$) only the conservative term is present, and we find that the spatial average of the total external force vanishes. In Fig. 1(d) we show the external potential $v_{\text{ext}}(x)$ that generates the conservative force contribution. As a convention, we have introduced an (irrelevant) shift of the energy scale, such that $v_{\rm ext}=0$ at $x=0$ for all four cases. As expected, in equilibrium $v_{\text{ext}}$ possesses a minimum at the location of the density peak. It turns out that in order to keep the density profile unchanged upon imposing the constant flux of particles in the $x$-direction, the external potential changes its shape very substantially. Both the minimum and the maximum move towards smaller values of $x$, i.e. against the direction of the flow, upon increasing $J_{0}$ (note the periodicity in $x$). Clearly, this behaviour is a direct consequence of keeping the density profile constant while increasing the flow through this density “landscape”. In order to rationalize this effect, consider first the case where an external potential generates the the density profile, Fig. 1(a), in equilibrium. If we now switch on a an additional constant (positive) external force contribution, the result will be a particle flow and the density profile will respond by shifting the density peak in the direction of the flow (results not shown). In our system the density profile is rather kept constant and the shifting of the density peak needs to be cancelled by the external conservative field, which hence necessarily develops the observed shift in the direction opposite to the flow. Besides quantifying the positional shift, cf. Fig. 1(d), we also observe a marked increase in the amplitude of the external potential contribution; hence stronger “ordering” forces, $-\nabla v_{\rm ext}$, are required in order to overcome the homogenizing effect of the flow.

III.2 Two-dimensional system

The iteration scheme is general and it is not restricted to effectively one-dimensional inhomogeneous systems. As a proof of concept, we construct the external force field that makes a two-dimensional density profile stationary. We choose the target velocity field to be

[TABLE]

with $d_{1},d_{2}=\rm const$. As above, a companion target density profile cannot be chosen arbitrarily, since the resulting current must satisfy the steady state condition, $\nabla\cdot\rho{\bf v}=0$. Given that (32) is divergence-free, $\nabla\cdot{\bf v}=0$, the steady state condition reduces to ${\bf v}\cdot\nabla\rho=0$. As an immediate first choice, we set

[TABLE]

which trivially satisfies the steady state condition. Note that (32) and (33) represent a conceptually highly interesting case of a homogeneous, bulk-fluid-like one-body density distribution, with “superimposed” flow (32) that is fully inhomogeneous on microscopic length scales.

Furthermore, as a second choice together with (32), we consider the target density profile

[TABLE]

such that $Y(y)$ is a spatially modulating function of the given form, $a_{0}$ is a constant such that $a_{0}<N/L^{2}$ (in order to ensure that the $\rho>0$, and $d_{0}=d_{1}/d_{2}$. Since $\nabla\rho$ is perpendicular to ${\bf v}$ for all ${\mathbf{r}}$, it is straightforward to show that (32), (34) and (35) also define a valid steady state.

For both target states (constant and non-constant density profiles) we use the inversion method to find the external force fields that renders a stationary situation. We set $N=30$, $L/\sigma=10$, and $k_{B}T/\epsilon=0.5$. For the target velocity profile, we set $d_{1}=d_{2}=\tau/\sigma$ in Eq. (32). For the inhomogeneous density profile we set $a_{0}=0.5N/L^{2}$ in Eq. (34).

The two Cartesian components of the velocity profile, obtained after $40$ BD iterations of the inversion method, are shown in Fig. 2 (a1) and (a2). The sampled velocity and density profiles coincide with the target profiles within the imposed numerical accuracy. The sampled density profiles are shown in Fig. 2(b1) (constant density profile) and Fig. 2(c1) (inhomogeneous density profile). The corresponding external force fields are presented in panels (b2-b3) and (c2-c3) of Fig. 2.

For the case of constant density, the $x$ component of the external force field $f_{\text{ext}}^{(x)}$, see Fig. 2(b2), is very similar in shape and magnitude to the $x$ component of the velocity profile, Fig. 2(a1). Given that the friction coefficient is set to $\gamma=1$, this means that $f_{\text{ext}}^{(x)}$ generates the flow in the $x$ direction (there are small differences between $f_{\text{ext}}^{(x)}$ and $\gamma v_{x}$ related to the $x$ component of the internal force field). The $y$ component of the external force field, shown in Fig. 2(b3), shows a small deviation from an average value which is consistent with the value of the flow in $y$. This deviation is expected, since we have imposed a constant density profile, and hence the external force has to balance the migration force Leighton and Acrivos (1987); Frank et al. (2003) that results from the shear field imposed by $v_{x}$. The $y$ component of the external force is inhomogeneous but the density profile is constant. Hence, the internal force must cancel the action of $f_{\text{ext}}^{(y)}$. This is a purely superadiabatic effect Stuhlmüller et al. (2018), which is completely neglected in the widely-used dynamical density functional theory (DDFT)  Marconi and Tarazona (1999) which rather predicts internal forces to vanish for situations of constant density. Extended versions of DDFT have been recently proposed to try to overcome these limitations, see e.g. Brader and Krüger (2011); Scacchi et al. (2016).

The target velocity profile is effectively one-dimensional, and the target density profile is constant. As a result the external force is also effectively one-dimensional. This is not the case when the target density profile is inhomogeneous. Then the $x$ and $y$ components of the external force field depend on both coordinates, see Fig. 2, panels (c2) and (c3). Now the external force field generates the flow and also sustains the density gradient. Clearly, the required force field, which generates the fully inhomogeneous steady state, is very complex, and simple physical reasoning, such as we could rely on in the former two cases, is insufficient to obtain even a qualitative, let alone (semi-)quantitative rationalization of the occurring physics.

III.3 Dynamic confinement

While the above examples demonstrate custom flow for steady states, we next turn to its implementation for full (time-dependent) nonequilibrium situations, as laid out in Sec. II.4. We hence aim to show that the concept is general and valid even for complex dynamics.

As a prototypical situation, we address the time evolution of a system, which in its initial state is a homogeneous equilibrium fluid, (with no external field acting in this initial state). At time $t=0$ we switch on a conservative external field, which represents the potential trap shown Fig. 1 panels (c) and (d) for the equilibrium case (green dot-dashed-line). Hence, the external force induces migration of particles towards the center of the system, as shown in panels (a) of Fig. 3 for the density (a1) and the current (a2).

The external field is static for $t>0$, see panel (a3), and its influence evolves the system from the homogeneous state to a confined state that features a well-defined, peaked density modulation, Fig. 3(a1). After only few Brownian times, a new equilibrium state is reached. The particle current almost vanishes already at time $\tau_{4}/\tau=2.2$, see Fig. 3(a2).

Using the time-dependent version of the custom flow method described in Sec. II.4, we chose to determine the time- and position-dependent external force field, $f_{\text{ext}}(t,x)$, that speeds up the dynamics of the system by a factor $\alpha>1$. That is, we find a system that evolves through the same temporal sequence of density profiles as those in Fig. 3(a), but doing so at a rate which $\alpha$ times faster. Hence in the new “fast forward” system the density profile $\rho_{\alpha}$ at time $t$ is the same as the density profile in the original system at time $\alpha t$. The current in the new system must be $\alpha$ times the current in the old system due to the continuity equation. That is, in the new system:

[TABLE]

In order to find the external field that induces the desired target dynamics, we discretize the time evolution at intervals $\Delta t_{\text{cg}}/\tau=0.01$, i.e. on a scale that is $10^{2}$ times larger than the time step of the BD simulation $\Delta t/\tau=10^{-4}$. At each coarse-graining time $t_{\rm cg}$ we run an iterative process to find the desired external field at that time. We use a linear regression to approximate the external field at every time $t$ between two consecutive coarse-graining times. The imposed coarse-graining time is a good compromise between accuracy and computational time. Since the external field does not vary profusely during one time interval only a few iterations ($<10$) are required at each $t_{\text{cg}}$. At each iteration we average over $10^{6}$ trajectories. Finally, we average the results over $50$ independent simulation runs.

The dynamics of the system with fast forward factor $\alpha=2$ is shown in Fig. 3(b). The external field that is required to speed up the dynamics by the chosen factor $\alpha=2$ is presented in Fig. 3(b3). As expected, $f_{\text{ext}}$ is now a time dependent field, the amplitude of which decreases monotonically during the time evolution. In the limit $t\rightarrow\infty$ the external field converges to that of the original system, since the final equilibrium state are required to be the same in both cases. In the Supplemental Material SuppMat , we show a movie of the time evolution and the required external field in a system which moves three times faster ($\alpha=3$) than the original system.

As a further example, we conceive a system in which the current is the same as in the original system (no speed up, $\alpha=1$), except for a prescribed additive constant $J_{0}$. As the divergence of the constant vanishes, it has no effect on the dynamics of the density distribution via the continuity equation. Hence the density profiles of both systems are the same at any time and the current profile in the new system is $J(t,x)+J_{0}$, where $J(t,x)$ is the original current. The required external field that produces such a dynamical evolution is shown in Fig. 3(c3), in which we have set $J_{0}\sigma\tau=0.5$. The external force field is again time dependent. The extrema of the force field are shifted with respect to their original location in the case of the static force field. This was expected given our above results for the steady state, Fig. 1. The amplitude of the force and the magnitude of the shift vary in a nontrivial way in time, as a result of a delicate balance between memory effects and the amplitude of the density modulation. At $t\rightarrow\infty$ the system reaches the same steady state as that shown in Fig. 1 (orange dashed line).

IV Discussion and conclusions

We have presented a numerical iterative method to systematically construct the specific form of the external force field which is required to drive a prescribed steady state in an overdamped Brownian many-body system. The same scheme can be used to find the conservative potential for which a given density profile is in thermodynamic equilibrium. In equilibrium the method is not restricted to BD systems.

An alternative approach has been previously developed for the equilibrium case Fortini et al. (2014) (also in quantum systems Thiele et al. (2008)). Although we have not studied the relative performance of the two methods systematically against each other, preliminary tests suggest that the current approach applied to equilibrium is both faster and more reliable. Whether the present method can or cannot be extended to quantum systems is an open and interesting question.

In all cases that we have analysed so far, the iteration processes have reliably converged. Nevertheless, if the initial guess for the external force is very far from the actual force field, it might be necessary to improve the simple fixed-point iteration scheme presented here in order to avoid possible divergent trajectories (i.e. sequences of ${\bf f}_{\rm ext}^{(k)}$). Using e.g. Anderson acceleration-like methods should constitute a possible improvement of the method. Variations of the presented iterative scheme, such as e.g. defining ${\bf f}_{\rm int}^{(k-1)}={\bf F}_{\rm int}^{(k-1)}/\rho^{(k)}$ instead of ${\bf f}_{\rm int}^{(k-1)}={\bf F}_{\rm int}^{(k-1)}/\rho$ in (11), also converge to the desired external force and might be useful in cases where convergence issues occur, which might be the case e.g. in the vicinity of dynamical phase transitions.

As the method requires a discretization of the space coordinate, therefore yields a discretized external force field. The quality of the spatial discretization (e.g. size of the bins) is an important parameter of the method. Although we have shown only one- and two-dimensional mono-component examples, the extension to three-dimensional systems and/or mixtures is straightforward.

In all cases, whether time-dependent nonequilibrium, nonequilibrium steady state, or in equilibrium, the custom flow method requires to sample the internal force field. Therefore, the practical implementation for hard-body systems is not as straightforward as it is in the case of soft interparticle potentials. For steady-state hard-body systems it might be easier to extend the equilibrium approach of Ref. Fortini et al. (2014) to nonequilibrium conditions.

The custom flow method allows complete control of the dynamics of a given system in both steady state and full nonequilibrium. Possible future applications include the investigation of time crystals Zhang et al. (2017); Moessner and Sondhi (2017) in BD systems, removal of flow instabilities via the application of external fields in a controlled manner, and obtaining a better understanding of memory effects by e.g. a systematic analysis of the external fields required to speed up and/or slow down a given dynamical process.

Acknowledgements.

This work is supported by the German Research Foundation (DFG) via SCHM 2632/1-1.

Appendix A Sampling the current in Brownian dynamics simulations

We briefly comment on three different methods to sample the one-body current ${\bf J}({\mathbf{r}})$ in Brownian dynamics simulations.

Method 1: Force balance equation. First, we propose here a new simple method to measure the current based on the exact one-body force density balance equation, Eq. (7). This equation provides an expression for ${\bf J}$ that can be used to directly sample the current. We need to sample: (i) the internal force density field ${\bf F}_{\text{int}}$ as an average over time (in steady state) or over many realizations (in case of time dependent situations), and (ii) the density profile. Then, using the density profile one can calculate the thermal diffusive term $-k_{B}T\nabla\rho$. Finally, the external force density field can either be calculated using the external force and the density profile (${\bf F}_{\text{ext}}=\rho{\bf f}_{\text{ext}}$) or sampled directly during the simulation.

Method 2: Numerical derivative of the position vector. The second method, proposed in Refs. Fortini et al. (2014); Krinninger et al. (2016), is based on calculating the velocity of the $i$th particle, ${\bf v}_{i}(t)$, via the numerical central derivative of the position vector

[TABLE]

Due to the stochastic nature of the motion it is crucial to use the central derivative to properly compute the velocity of the particles, Eq. (37). Forward and backward derivatives give different results that are not consistent with the value of the current obtained by the alternative methods presented here.

A spatially resolved average of ${\bf v}_{i}$, Eq. (37), yields the one-body current profile:

[TABLE]

where $\langle\cdot\rangle$ indicates an average over either many different realizations or time in the case of a steady state.

To sample ${\bf v}_{i}$ more efficiently it is convenient to rewrite Eq. (37) as Krinninger et al. (2016)

[TABLE]

where $\Delta{\mathbf{r}}_{i}(t)={\mathbf{r}}_{i}(t+\Delta t)-{\mathbf{r}}_{i}(t)$. Plugging Eq. (30) into  (39) results in

[TABLE]

The spatially resolved average of $\boldsymbol{\eta}_{i}(t)$ vanishes at any space point since $\boldsymbol{\eta}_{i}$ is a gaussian random force, and therefore this average correlates the random force at time $t$ with the position at the same time $t$. In contrast, it is important to realize that the spatially resolved average of $\boldsymbol{\eta}_{i}(t-\Delta t)$ does not vanish in general, since it correlates the random force at time $t-\Delta t$ with the current position at time $t$.

Method 3: Continuity equation. The continuity equation, Eq. (5) provides an alternative route to compute the current in non-steady state situations, as shown in Ref. Fortini et al. (2014). Having sampled the density profile at different times $t$ and $t^{\prime}$, we can compute the numerical time derivative of the density profile, which must be equal to the divergence of the current. In effectively one-dimensional systems the result can be integrated in space and yields the one-body current profile (line integrals or other inversion methods are required in a higher dimensional space). This method yields the current profile up to an additive constant. If the actual value of the current at a given space point is known, then one can easily determine the missing additive constant. For instance, if the system is in contact with a hard wall the current at the hard wall must vanish. To use this method we need to sample $\rho$ at two times $t$ and $t^{\prime}$ separated by a time interval $\Delta_{t}$. In our experience a value $\Delta_{t}\approx 10^{2}\Delta t$ with $\Delta t$ the time step of the BD simulation provides good results.

We have checked the three methods presented above give the same one-body current profile within the inherent numerical accuracy of each procedure.

Appendix B Numerical details

We provide details of our precise implementation of the iterative scheme.

Each iteration step (11) of the nonequilibrium inversion method requires carrying out one BD simulation run for the given parameters and given external force field. Before acquiring the data, we let the system reach a steady state during $10^{2}\tau$. Then, at each iteration $k$, we sample the internal force density during a given sampling time $\tau_{k}$. The sampling time has a direct impact on both the statistical quality of the sampled internal force field (which is required for the next iteration) as well as on the performance of the method. Instead of using the same sampling time at each iteration, we find it preferable to start with short simulation runs and increase the run length at every iteration. Hence, at iteration $k$, we fix the sampling simulation time $\tau_{k}$ to

[TABLE]

i.e., we double the runlength every three iterations. The total time of the first iteration is set to $\tau_{0}/\tau=100$. Finally, we average over several ($10-100$) realizations of the iteration scheme, Eq. (11), to improve the statistics.

Fig. 4 illustrates the iterative process for the effectively one-dimensional system with target profile given by Eq. (31) and target current $J_{0}\tau\sigma=0.1$. Less than 10 iterations suffice to get a very good estimate of the external force, and after $k=40$ iterations the sampled density, and current profiles are almost indistinguishable from the corresponding target profiles. The results have been obtained by averaging over $25$ BD realizations of the iterative scheme. We show in Fig. 5(a) the evolution of the error of the density and the current profile during the iterative process, cf. Eqs. (13) and (14).

For the equilibrium situation of the effectively one-dimensional system shown in Fig. 1 ($J_{0}=0$) we have used Monte Carlo simulations to implement the iterative scheme, cf. (19) and (21). For completeness, we also show the efficiency of the method in equilibrium in Fig. 5(b), where we plot the difference between the sampled and the error in the density profile as a function of the number of iterations and the number of Monte Carlo sweeps (MCS). Each MCS is an attempt to sequentially and individually move all the particles in the system. We find it convenient to increase the number of MCS during the iterative process. We begin with $10^{4}$ MCS at iteration $k=1$, and increase the number of MCS after every iteration such that it doubles every three iterations. Before acquiring data we equilibrate the system by running $10^{4}$ MCS. As in the non-equilibrium steady state case, we improve the statistics by averaging over $25$ realizations.