Custom flow in overdamped Brownian Dynamics

Custom flow in overdamped Brownian Dynamics

Daniel de las Heras Email: Theoretische Physik II, Physikalisches Institut, Universität Bayreuth, D-95440 Bayreuth, Germany    Johannes Renner Theoretische Physik II, Physikalisches Institut, Universität Bayreuth, D-95440 Bayreuth, Germany    Matthias Schmidt Email: Theoretische Physik II, Physikalisches Institut, Universität Bayreuth, D-95440 Bayreuth, Germany
5 August 2018, revised version: 12 October 2018

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.

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 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 and time is given by


where the angles denote a statistical average, the expression inside the angles is the microscopic density operator, with indicating the Dirac distribution, and denoting the position of particle . In the Fokker-Planck picture, the information required for carrying out the average is encoded in the many-body probability distribution of finding microstate at time . The average is then defined as


where the integral runs over configuration space, i.e. each is integrated over the system volume. The time evolution of is determined by the Smoluchowski equation . Here the configuration space velocity of particle  is given on the many-body level by the instantaneous relation


where is the friction constant against the implicit solvent, is the Boltzmann constant, denotes absolute temperature, is the partial derivative with respect to , the interparticle interaction potential is denoted by , and 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 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


where , at time , 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 in (4).

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


The many-body coupling in (3) arises due to the presence of the internal interaction potential . On the one-body level, it is hence natural to define a corresponding internal force density field via


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

By multiplying (3) with , summing over , averaging according to (2), and identifying the one-body fields (1), (4), and (6), it is straightforward to show that


which we use as a basis for the non-equilibrium inversion procedure. Defining the microscopic velocity profile simply as the ratio between current profile and density profile,


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


Here the internal force field is defined as the internal force density (6) normalized with the density profile, i.e. .

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 is readily obtained by (numerical) differentiation of the density profile. The internal force density field 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, in (6).

In typical applications, the external force field is prescribed and and emerge as a result of the coupled many-body dynamics. In the following, we address the inverse problem of prescribing and a priory and calculating the required form of 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 and 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 , cf. (8), must be divergence-free, , which follows from (5) and represents a necessary condition on the allowed set of target functions , . The external force field that makes the target profiles stationary is obtained from first rearranging (9) as


where the internal force field, , is the only unknown quantity on the right hand side, as and are known input quantities. Here is from (6) in steady state.

In order to determine , and hence via (10), we proceed in two steps. First we present a fixed-point iterative scheme to solve (10), in which the -th iteration step is defined via


where the targets, and , are kept fixed for all steps . Here , where is the internal force density sampled in the previous iteration step, . Data for was obtained from BD sampling under the prescribed external force field . In order to initialize the iteration scheme, we set the external force field at step simply as


which is the exact external force field for the case of an ideal gas. Prescribing (12) allows to sample in BD, and then use as the required input for iteration step 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, and ; details on how to sample the current in BD are provided in Appendix A. As a criterion for judging the eventual convergence of 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 , i.e.,


where is the system volume, and 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 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 of the one-body density and current or analogously, viz. (8), of the density and the velocity field. 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 (), an excess (over ideal gas) part due to the internal interactions () and an external power () contribution, according to .

PFT implies that is a unique functional of density and current distributions, or equivalently of density and velocity profile. In particular, can be expressed as a functional derivative of the intrinsic excess (over ideal gas) free power functional , as


where the density distribution is kept fixed upon the variation, at fixed time .

Typically, one would split further into adiabatic and superadiabatic contributions, , where is the time derivative of the equilibrium excess (over ideal gas) Helmholtz free energy functional, and 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 . Recall that this is directly accessible via Eq. (6) in BD simulations.

The fact that 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,


where the right hand side is the force field functional (15) evaluated at the target profiles and ; 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


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 ; no hidden dependence on occurs on the right hand side. Recall that from (15), the internal force field depends solely on the “kinematic” fields and , 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 and prescribe in order to find the corresponding external force field . The external force field is necessary of conservative, gradient form, , where 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


which constitutes an explicit expression for the specific external force field that generates the given target profile 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


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


We then sample in equilibrium, under the external force , 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 instead of the force, . In systems that effectively depend on only one coordinate, say , the potential at each iteration can be easily obtained from the internal force profile, by performing a one-dimensional spatial integral


In two- and three-dimensional systems a line integral or, more generally, the use of an inverse operator 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


Inserting this into (18) yields


which is an explicit expression for the external force field, given 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 , where , 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 , the internal force field may be written as


in terms of the density distribution at equilibrium with that external force. Then Eq. (19) may be written as the change in the force, along the iterative procedure,


that vanishes when the target density is achieved, . The integration of (25), to get the change in external potential, and the linear expansion


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 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 and velocity field 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 we iterate the value of the external field according to


where and are the target fields, which enter via their values at time . The time is kept fixed under the iteration 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:


where indicates the time step previous to .

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 is known, and as consistency check, can be input into a “bare” non-steady BD run, in order to validate that the targets and 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),


where the parameters and set the energy and length scales, respectively, and indicates the center-center distance of the particle pair. The cutoff distance 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 with periodic boundary conditions along both directions. Using the standard Euler algorithm, the Langevin equation of motion is integrated in time via


where is a delta-correlated Gaussian random displacement with standard deviation in accordance with the fluctuation-dissipation theorem. Here, is the integration time step that we set to with the unit of time; the friction constant is set to .

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 , and the flow direction is along the -axis (i.e. with no shear motion occurring). Then the steady state condition reduces to the requirement of the current being constant, , where and is the unit vector in the -direction. Hence from (8) the velocity and the density profile possess a reciprocal relationship: .

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


with and such that . See the density profile in Fig. 1(a). The temperature is set to .

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 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 . 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 (equilibrium), and . 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 () 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 that generates the conservative force contribution. As a convention, we have introduced an (irrelevant) shift of the energy scale, such that at for all four cases. As expected, in equilibrium 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 -direction, the external potential changes its shape very substantially. Both the minimum and the maximum move towards smaller values of , i.e. against the direction of the flow, upon increasing (note the periodicity in ). 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, , are required in order to overcome the homogenizing effect of the flow.

Figure 1: One-body density (a) and current profiles (b), external force field (b), and external potential (c) as a function of the -coordinate in a system with target density profile with and after iterations. The inset in (a) shows the difference between target and sampled density profiles. Results are shown for different values of the target current: (blue dotted line), (orange dashed line), (violet solid line), and (green dot-dashed line), which is in equilibrium. The inset in (b) is a close-view for . Two-dimensional system with particles in a periodic box of side at . Data obtained by averaging over BD realizations (MC realizations in equilibrium).

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


with . As above, a companion target density profile cannot be chosen arbitrarily, since the resulting current must satisfy the steady state condition, . Given that (32) is divergence-free, , the steady state condition reduces to . As an immediate first choice, we set


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


such that is a spatially modulating function of the given form, is a constant such that (in order to ensure that the , and . Since is perpendicular to for all , it is straightforward to show that (32), (34) and (35) also define a valid steady state.

Figure 2: component (a1) and component (a2) of the velocity profile sampled after iterations using BD simulations. (b1) Sampled density profile for the steady state with constant density profile. (b2) and (b3) components of the external force field that produces the steady state with constant density. (c1) Sampled density profile for the steady state with inhomogeneous density profile. (c2) and (c3) components of the external force field that generates the steady state with inhomogeneous density. In both steady states we set , , and . The bin size is set to in both directions and the origin or coordinates is located in the middle of the box. Results are averages over BD realizations.

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 , , and . For the target velocity profile, we set in Eq. (32). For the inhomogeneous density profile we set in Eq. (34).

The two Cartesian components of the velocity profile, obtained after 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 component of the external force field , see Fig. 2(b2), is very similar in shape and magnitude to the component of the velocity profile, Fig. 2(a1). Given that the friction coefficient is set to , this means that generates the flow in the direction (there are small differences between and related to the component of the internal force field). The 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 . 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 . The component of the external force is inhomogeneous but the density profile is constant. Hence, the internal force must cancel the action of . 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 and 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 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 , 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 , 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, , that speeds up the dynamics of the system by a factor . 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 times faster. Hence in the new “fast forward” system the density profile at time is the same as the density profile in the original system at time . The current in the new system must be times the current in the old system due to the continuity equation. That is, in the new system:


In order to find the external field that induces the desired target dynamics, we discretize the time evolution at intervals , i.e. on a scale that is times larger than the time step of the BD simulation . At each coarse-graining time 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 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 () are required at each . At each iteration we average over trajectories. Finally, we average the results over independent simulation runs.

Figure 3: Inversion in full nonequilibrium as demonstrated by dynamics of confinement. Density profile (left column), current profile (middle column), and external force field (right column) as a function of for three different situations. (a) Time evolution of the system, density (a1) and current (a2), after switching on a static external field (a3). At the system is in equilibrium with vanishing external field. Results at four times are shown: , , , and . At the system is near equilibrium with the applied external field (a3). Profiles are obtained by averaging over different trajectories. In panels (b) we show a system that evolves following the same dynamics as in (a) but two times faster. The current (b2) is therefore twice as large as the current in the original system, and the required external field (b3) is time dependent. A movie showing a system that evolves three times faster is presented in the Supplemental Material SuppMat (). Panels (c) show the time evolution in a system that reproduces the same dynamics as (a) but with a global motion towards the right (note that the spatial average of the current (c2) at any time is ). The required external field (c3) is time-dependent. The horizontal and vertical dashed lines in the plots of the external field are drawn to help the comparison between systems. In all cases we set , , and .

The dynamics of the system with fast forward factor is shown in Fig. 3(b). The external field that is required to speed up the dynamics by the chosen factor is presented in Fig. 3(b3). As expected, is now a time dependent field, the amplitude of which decreases monotonically during the time evolution. In the limit 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 () 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, ), except for a prescribed additive constant . 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 , where 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 . 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 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 ). 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 instead of 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.

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 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 that can be used to directly sample the current. We need to sample: (i) the internal force density field 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 . Finally, the external force density field can either be calculated using the external force and the density profile () 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 th particle, , via the numerical central derivative of the position vector


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 , Eq. (37), yields the one-body current profile:


where indicates an average over either many different realizations or time in the case of a steady state.

To sample more efficiently it is convenient to rewrite Eq. (37) as Krinninger et al. (2016)


where . Plugging Eq. (30) into  (39) results in


The spatially resolved average of vanishes at any space point since is a gaussian random force, and therefore this average correlates the random force at time with the position at the same time . In contrast, it is important to realize that the spatially resolved average of does not vanish in general, since it correlates the random force at time with the current position at time .

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 and , 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 at two times and separated by a time interval . In our experience a value with 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.

Figure 4: One-body density profile (a), one-body current (b), and external force field (c) as a function of the -coordinate for different number of iterations, , as indicate in the legend of (a). The target current is set to . The target density profile is with and , which is practically identical to the sampled density profile after iterations (violet solid line). Two-dimensional system with particles in a periodic box of side at temperature . The data has been obtained by averaging over BD realizations. The total simulation time of iteration of one BD realization is set to , with . The bin size is .

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 . Then, at each iteration , we sample the internal force density during a given sampling time . 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 , we fix the sampling simulation time to


i.e., we double the runlength every three iterations. The total time of the first iteration is set to . Finally, we average over several () 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 . Less than 10 iterations suffice to get a very good estimate of the external force, and after iterations the sampled density, and current profiles are almost indistinguishable from the corresponding target profiles. The results have been obtained by averaging over 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).

Figure 5: (a) Scaled error of the density profile (green circles) and of the current profile (yellow squares) as a function of the iteration number (bottom axis) and the scaled simulation time of the iteration (upper axis) in non-equilibrium steady state using BD simulations. (b) Scaled error of the density profile as a function of the iteration number (bottom axis) and the number of Monte Carlo steps at iteration (upper axis) in equilibrium using MC simulations. In both panels, the data has been obtained by averaging over realizations. Note the logarithmic scale of the vertical axis and of the upper horizontal axis.

For the equilibrium situation of the effectively one-dimensional system shown in Fig. 1 () 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 MCS at iteration , 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 MCS. As in the non-equilibrium steady state case, we improve the statistics by averaging over realizations.


  • Löwen (2008) Hartmut Löwen, “Colloidal dispersions in external fields: recent developments,” J. Phys. Condens. Matter 20, 404201 (2008).
  • Erbe et al. (2008) A. Erbe, M. Zientara, L. Baraban, C. Kreidler,  and P. Leiderer, “Various driving mechanisms for generating motion of colloidal particles,” J. Phys. Condens. Matter 20, 404215 (2008).
  • Perrin (1909) F. Perrin, Ann. Chem. Phys. 18, 5 (1909).
  • Köppl et al. (2006) M. Köppl, P. Henseler, A. Erbe, P. Nielaba,  and P. Leiderer, “Layer reduction in driven 2d-colloidal systems through microchannels,” Phys. Rev. Lett. 97, 208302 (2006).
  • Simeonova and Kegel (2004) Nikoleta B. Simeonova and Willem K. Kegel, “Gravity-induced aging in glasses of colloidal hard spheres,” Phys. Rev. Lett. 93, 035701 (2004).
  • Palacci et al. (2010) Jérémie Palacci, Cécile Cottin-Bizonne, Christophe Ybert,  and Lydéric Bocquet, “Sedimentation and effective temperature of active colloidal suspensions,” Phys. Rev. Lett. 105, 088304 (2010).
  • Velev and Bhatt (2006) Orlin D. Velev and Ketan H. Bhatt, “On-chip micromanipulation and assembly of colloidal particles by electric fields,” Soft Matter 2, 738 (2006).
  • O’Brien and White (1978) Richard W. O’Brien and Lee R. White, ‘‘Electrophoretic mobility of a spherical colloidal particle,” J. Chem. Soc., Faraday Trans. 2 74, 1607 (1978).
  • Tierno et al. (2007) Pietro Tierno, Ramanathan Muruganathan,  and Thomas M. Fischer, “Viscoelasticity of dynamically self-assembled paramagnetic colloidal clusters,” Phys. Rev. Lett. 98, 028301 (2007).
  • Snezhko and Aranson (2011) Alexey Snezhko and Igor S. Aranson, “Magnetic manipulation of self-assembled colloidal asters,” Nat. Mater. 10, 698 (2011).
  • Loehr et al. (2016) Johannes Loehr, Michael Loenne, Adrian Ernst, Daniel de las Heras,  and Thomas M Fischer, “Topological protection of multiparticle dissipative transport,” Nat. Commun. 7, 11745 (2016).
  • Bender and Wagner (1996) Jonathan Bender and Norman J. Wagner, “Reversible shear thickening in monodisperse and bidisperse colloidal dispersions,” J. Rheol. 40, 899 (1996).
  • Wagner and Brady (2009) Norman J Wagner and John F Brady, “Shear thickening in colloidal dispersions,” Phys. Today 62, 27 (2009).
  • Besseling et al. (2007) R. Besseling, Eric R. Weeks, A. B. Schofield,  and W. C. K. Poon, “Three-dimensional imaging of colloidal glasses under steady shear,” Phys. Rev. Lett. 99, 028301 (2007).
  • Ruffner and Grier (2012) David B. Ruffner and David G. Grier, ‘‘Optical forces and torques in nonuniform beams of light,” Phys. Rev. Lett. 108, 173602 (2012).
  • Sun et al. (2009) Bo Sun, Jiayi Lin, Ellis Darby, Alexander Y. Grosberg,  and David G. Grier, “Brownian vortexes,” Phys. Rev. E 80, 010401 (2009).
  • Lutz et al. (2004) Christoph Lutz, Markus Kollmann,  and Clemens Bechinger, “Single-file diffusion of colloids in one-dimensional channels,” Phys. Rev. Lett. 93, 026001 (2004).
  • Mermin (1965) N. David Mermin, “Thermal properties of the inhomogeneous electron gas,” Phys. Rev. 137, A1441 (1965).
  • 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 (1979).
  • Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864 (1964).
  • Evans et al. (2016) Robert Evans, Martin Oettel, Roland Roth,  and Gerhard Kahl, “New developments in classical density functional theory,” J. Phys.: Condens. Matter 28, 240401 (2016).
  • Schmidt and Brader (2013) Matthias Schmidt and Joseph M. Brader, “Power functional theory for Brownian dynamics,” J. Chem. Phys. 138 (2013).
  • Fortini et al. (2014) Andrea Fortini, Daniel de las Heras, Joseph M. Brader,  and Matthias Schmidt, “Superadiabatic forces in Brownian many-body dynamics,” Phys. Rev. Lett. 113, 167801 (2014).
  • Borgis et al. (2013) Daniel Borgis, Roland Assaraf, Benjamin Rotenberg,  and Rodolphe Vuilleumier, ‘‘Computation of pair distribution functions and three-dimensional densities with a reduced variance principle,” Mol. Phys. 111, 3486 (2013).
  • de las Heras and Schmidt (2018) Daniel de las Heras and Matthias Schmidt, “Better than counting: Density profiles from force sampling,” Phys. Rev. Lett. 120, 218001 (2018).
  • Weeks et al. (1971) John D. Weeks, David Chandler,  and Hans C. Andersen, “Role of repulsive forces in determining the equilibrium structure of simple liquids,” J. Chem. Phys. 54, 5237 (1971).
  • Leighton and Acrivos (1987) David Leighton and Andreas Acrivos, “The shear-induced migration of particles in concentrated suspensions,” J. Fluid Mech. 181, 415–439 (1987).
  • Frank et al. (2003) Martin Frank, Douglas Anderson, Eric R. Weeks,  and Jeffrey F. Morris, “Particle migration in pressure-driven flow of a Brownian suspension,” J. Fluid Mech. 493, 363–378 (2003).
  • Stuhlmüller et al. (2018) Nico C. X. Stuhlmüller, Tobias Eckert, Daniel de las Heras,  and Matthias Schmidt, “Structural nonequilibrium forces in driven colloidal systems,” Phys. Rev. Lett. 121, 098002 (2018).
  • Marconi and Tarazona (1999) U. M. B. Marconi and P. Tarazona, “Dynamic density functional theory of fluids,” J. Chem. Phys. 110, 8032 (1999).
  • Brader and Krüger (2011) Joseph M. Brader and Matthias Krüger, “Density profiles of a colloidal liquid at a wall under shear flow,” Mol. Phys. 109, 1029 (2011).
  • Scacchi et al. (2016) Alberto Scacchi, Matthias Krüger,  and Joseph M. Brader, “Driven colloidal fluids: construction of dynamical density functional theories from exactly solvable limits,” J. Phys. Condens. Matter 28, 244023 (2016).
  • (33) See Supplemental Material for a movie with an example in time-dependent situations
  • Thiele et al. (2008) M. Thiele, E. K. U. Gross,  and S. Kümmel, “Adiabatic approximation in nonperturbative time-dependent density-functional theory,” Phys. Rev. Lett. 100, 153004 (2008).
  • Krinninger et al. (2016) Philip Krinninger, Matthias Schmidt,  and Joseph M. Brader, “Nonequilibrium phase behavior from minimization of free power dissipation,” Phys. Rev. Lett. 117, 208003 (2016).
  • Zhang et al. (2017) J Zhang, PW Hess, A Kyprianidis, P Becker, A Lee, J Smith, G Pagano, I-D Potirniche, Andrew C Potter, A Vishwanath, et al., “Observation of a discrete time crystal,” Nature 543, 217 (2017).
  • Moessner and Sondhi (2017) Roderich Moessner and SL Sondhi, “Equilibration and order in quantum floquet matter,” Nat. Phys. 13, 424 (2017).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description