Search

The last few years have witnessed an explosion of new numerical methods for filament hydrodynamics. Aside from their ubiquity in biology, physics, and engineering, filaments present unique challenges from an applied-mathematical point of view. Their slenderness, inextensibility, semiflexibility, and meso-scale nature all require numerical methods that can handle multiple lengthscales in the presence of constraints. Accounting for Brownian motion while keeping the dynamics in detailed balance and on the constraint is difficult, as is including a background solvent, which couples the dynamics of multiple filaments together in a suspension. In this paper, we present a simulation platform for deterministic and Brownian inextensible filament dynamics which includes nonlocal fluid dynamics and steric repulsion. For nonlocal hydrodynamics, we define the mobility on a single filament using line integrals of Rotne-Prager-Yamakawa regularized singularities, and numerically preserve the symmetric positive definite property by using a thicker regularization width for the nonlocal integrals than for the self term. For steric repulsion, we introduce a soft local repulsive potential defined as a double-integral over two filaments, then present a scheme to identify and evaluate the nonzero components of the integrand. Using a temporal integrator developed in previous work, we demonstrate that Langevin dynamics sample from the equilibrium distribution of free filament shapes, and that the modeling error in using the thicker regularization is small. We conclude with two examples, sedimenting filaments and cross-linked fiber networks, in which nonlocal hydrodynamics does and does not generate long-range flow fields, respectively. In the latter case, we show that the effect of hydrodynamics can be accounted for through steric repulsion.

Semiflexible slender filaments are ubiquitous in nature and cell biology, including in the cytoskeleton, where reorganization of actin filaments allows the cell to move and divide. Most methods for simulating semiflexible inextensible fibers/polymers are based on discrete (bead-link or blob-link) models, which become prohibitively expensive in the slender limit when hydrodynamics is accounted for. In this paper, we develop a novel coarse-grained approach for simulating fluctuating slender filaments with hydrodynamic interactions. Our approach is tailored to relatively stiff fibers whose persistence length is comparable to or larger than their length, and is based on three major contributions. First, we discretize the filament centerline using a coarse non-uniform Chebyshev grid, on which we formulate a discrete constrained Gibbs-Boltzmann equilibrium distribution and overdamped Langevin equation. Second, we define the hydrodynamic mobility at each point on the filament as an integral of the Rotne-Prager-Yamakawa kernel along the centerline, and apply a spectrally-accurate quadrature to accurately resolve the hydrodynamics. Third, we propose a novel midpoint temporal integrator which can correctly capture the Ito drift terms that arise in the overdamped Langevin equation. We verify that the equilibrium distribution for the Chebyshev grid is a good approximation of the blob-link one, and that our temporal integrator samples the equilibrium distribution for sufficiently small time steps. We also study the dynamics of relaxation of an initially straight filament, and find that as few as 12 Chebyshev nodes provides a good approximation to the dynamics while allowing a time step size two orders of magnitude larger than a resolved blob-link simulation. We conclude by studying how bending fluctuations aid the process of bundling in cross-linked networks of semiflexible fibers.

We develop a linearly-scaling variant of the Force Coupling Method [K. Yeo and M. R. Maxey, J. Fluid Mech. 649, 205-231 (2010)] for computing hydrodynamic interactions among particles confined to a doubly-periodic geometry with either a single bottom wall or two walls (slit channel) in the aperiodic direction. Our spectrally-accurate Stokes solver uses the Fast Fourier Transform (FFT) in the periodic $xy$ plane and Chebyshev polynomials in the aperiodic $z$ direction normal to the wall(s). We decompose the problem into two problems. The first is a doubly-periodic subproblem in the presence of particles (source terms) with free-space boundary conditions in the $z$ direction, which we solve by borrowing ideas from a recent method for rapid evaluation of electrostatic interactions in doubly-periodic geometries [O. Maxian, R. P. Pel\'aez, L. Greengard and A. Donev, J. Chem. Phys. 154, 204107 (2021)]. The second is a correction subproblem to impose the boundary conditions on the wall(s). Instead of the traditional Gaussian kernel, we use the exponential of a semicircle kernel to model the source terms (body force) due to the presence of particles, and provide optimum values for the kernel parameters that ensure a given hydrodynamic radius with at least two digits of accuracy and rotational and translational invariance. The computation time of our solver, which is implemented in graphical processing units, scales linearly with the number of particles, and allows computations with about a million particles in less than a second for a sedimented layer of colloidal microrollers. We find that in a slit channel, a driven dense suspension of microrollers maintains the same two-layer structure as above a single wall, but moves at a substantially lower collective speed due to increased confinement.

Slender fibers are ubiquitous in biology, physics, and engineering, with prominent examples including bacterial flagella and cytoskeletal fibers. In this setting, slender body theories (SBTs), which give the resistance on the fiber asymptotically in its slenderness $\epsilon$, are useful tools for both analysis and computations. However, a difficulty arises when accounting for twist and cross-sectional rotation: because the angular velocity of a filament can vary depending on the order of magnitude of the applied torque, asymptotic theories must give accurate results for rotational dynamics over a range of angular velocities. In this paper, we first survey the challenges in applying existing SBTs, which are based on either singularity or full boundary integral representations, to rotating filaments, showing in particular that they fail to consistently treat rotation-translation coupling in curved filaments. We then provide an alternative approach which approximates the three-dimensional dynamics via a one-dimensional line integral of Rotne-Prager-Yamakawa regularized singularities. While unable to accurately resolve the flow field near the filament, this approach gives a grand mobility with symmetric rotation-translation and translation-rotation coupling, making it applicable to a broad range of angular velocities. To restore fidelity to the three-dimensional filament geometry, we use our regularized singularity model to inform a simple empirical equation which relates the mean force and torque along the filament centerline to the translational and rotational velocity of the cross section. The single unknown coefficient in the model is estimated numerically from three-dimensional boundary integral calculations on a rotating, curved filament., Comment: 42 pages, 4 figures

In swimming microorganisms and the cell cytoskeleton, inextensible fibers resist bending and twisting, and interact with the surrounding fluid to cause or resist large-scale fluid motion. In this paper, we develop a novel numerical method for the simulation of cylindrical fibers by extending our previous work on inextensible bending fibers [Maxian et al., Phys. Rev. Fluids 6 (1), 014102] to fibers with twist elasticity. In our "Euler" model, twist is a scalar function that measures the deviation of the fiber cross section relative to a twist-free frame, the fiber exerts only torque parallel to the centerline on the fluid, and the perpendicular components of the rotational fluid velocity are discarded in favor of the translational velocity. In the first part of this paper, we justify this model by comparing it to another commonly-used "Kirchhoff" formulation where the fiber exerts both perpendicular and parallel torque on the fluid, and the perpendicular angular fluid velocity is required to be consistent with the translational fluid velocity. We then develop a spectral numerical method for the hydrodynamics of the Euler model. We define hydrodynamic mobility operators using integrals of the Rotne-Prager-Yamakawa tensor, and evaluate these integrals through a novel slender-body quadrature, which requires on the order of 10 points along the fiber to obtain several digits of accuracy. We demonstrate that this choice of mobility removes the unphysical negative eigenvalues in the translation-translation mobility associated with asymptotic slender body theories, and ensures strong convergence of the fiber velocity and weak convergence of the fiber constraint forces. We pair the spatial discretization with a semi-implicit temporal integrator to confirm the negligible contribution of twist elasticity to the relaxation dynamics of a bent fiber and study the instability of a twirling fiber., Comment: 76 pages, 13 figures. See https://github.com/stochasticHydroTools/SlenderBody for MATLAB codes that accompany this paper and https://cims.nyu.edu/~om759/TwistVideos/index.html for supplementary animations

We develop a fast method for computing the electrostatic energy and forces for a collection of charges in doubly-periodic slabs with jumps in the dielectric permittivity at the slab boundaries. Our method achieves spectral accuracy by using Ewald splitting to replace the original Poisson equation for nearly-singular sources with a smooth far-field Poisson equation, combined with a localized near-field correction. Unlike existing spectral Ewald methods, which make use of the Fourier transform in the aperiodic direction, we recast the problem as a two-point boundary value problem in the aperiodic direction for each transverse Fourier mode, for which exact analytic boundary conditions are available. We solve each of these boundary value problems using a fast, well-conditioned Chebyshev method. In the presence of dielectric jumps, combining Ewald splitting with the classical method of images results in smoothed charge distributions which overlap the dielectric boundaries themselves. We show how to preserve spectral accuracy in this case through the use of a harmonic correction which involves solving a simple Laplace equation with smooth boundary data. We implement our method on Graphical Processing Units, and combine our doubly-periodic Poisson solver with Brownian Dynamics to study the equilibrium structure of double layers in binary electrolytes confined by dielectric boundaries. Consistent with prior studies, we find strong charge depletion near the interfaces due to repulsive interactions with image charges, which points to the need for incorporating polarization effects in understanding confined electrolytes, both theoretically and computationally., Comment: 60 pages, 7 figures

We study the driven collective dynamics of a colloidal monolayer sedimentating down an inclined plane. The action of the gravity force parallel to the bottom wall creates a flow around each colloid, and the hydrodynamic interactions among the colloids accelerate the sedimentation as the local density increases. This leads to the creation of a universal "triangular" inhomogeneous density profile, with a traveling density shock at the leading front moving in the downhill direction. Unlike density shocks in a colloidal monolayer driven by applied torques rather than forces [Phys. Rev. Fluids, 2(9):092301, 2017], the density front during sedimentation remains stable over long periods of time even though it develops a roughness on the order of tens of particle diameters. Through experimental measurements and particle-based computer simulations, we find that the Burgers equation can model the density profile along the sedimentation direction as a function of time remarkably well, with a modest improvement if the nonlinear conservation law accounts for the sub-linear dependence of the collective sedimentation velocity on density.

Every animal cell is filled with a cytoskeleton, a dynamic gel made of inextensible fibers, such as microtubules, actin fibers, and intermediate filaments, all suspended in a viscous fluid. Numerical simulation of this gel is challenging because the fiber aspect ratios can be as large as $10^4$. We describe a new method for rapidly computing the dynamics of inextensible slender filaments in periodically-sheared Stokes flow. The dynamics of the filaments are governed by a nonlocal slender body theory which we partially reformulate in terms of the Rotne-Prager-Yamakawa hydrodynamic tensor. To enforce inextensibility, we parameterize the space of inextensible fiber motions and strictly confine the dynamics to the manifold of inextensible configurations. To do this, we introduce a set of Lagrange multipliers for the tensile force densities on the filaments and impose the constraint of no virtual work in an $L^2$ weak sense. We augment this approach with a spectral discretization of the local and nonlocal slender body theory operators which is linear in the number of unknowns and gives improved spatial accuracy over approaches based on solving a line tension equation. For dynamics, we develop a second-order semi-implicit temporal integrator which requires at most a few evaluations of nonlocal hydrodynamics and a few block diagonal linear solves per time step. After demonstrating the improved accuracy and robustness of our approach through numerical examples, we apply our formulation to a permanently cross-linked actin mesh in a background oscillatory shear flow. We observe a characteristic frequency at which the network transitions from quasi-static, primarily elastic behavior to dynamic, primarily viscous behavior. We find that nonlocal hydrodynamics increases the viscous modulus by as much as 25%, even for semi-dilute fiber suspensions., Comment: 84 pages, 12 figures

In this paper we develop a methodology for the mesoscale simulation of strong electrolytes. The methodology is an extension of the Fluctuating Immersed Boundary (FIB) approach that treats a solute as discrete Lagrangian particles that interact with Eulerian hydrodynamic and electrostatic fields. In both cases the Immersed Boundary (IB) method of Peskin is used for particle-field coupling. Hydrodynamic interactions are taken to be overdamped, with thermal noise incorporated using the fluctuating Stokes equation, including a "dry diffusion" Brownian motion to account for scales not resolved by the coarse-grained model of the solvent. Long range electrostatic interactions are computed by solving the Poisson equation, with short range corrections included using a novel immersed-boundary variant of the classical Particle-Particle Particle-Mesh (P3M) technique. Also included is a short range repulsive force based on the Weeks-Chandler-Andersen (WCA) potential. The new methodology is validated by comparison to Debye-H{\"u}ckel theory for ion-ion pair correlation functions, and Debye-H{\"u}ckel-Onsager theory for conductivity, including the Wein effect for strong electric fields. In each case good agreement is observed, provided that hydrodynamic interactions at the typical ion-ion separation are resolved by the fluid grid., Comment: 30 pages, 12 figures, 2 tables

We perform detailed computational and experimental measurements of the driven dynamics of a dense, uniform suspension of sedimented microrollers driven by a magnetic field rotating around an axis parallel to the floor. We develop a lubrication-corrected Brownian Dynamics method for dense suspensions of driven colloids sedimented above a bottom wall. The numerical method adds lubrication friction between nearby pairs of particles, as well as particles and the bottom wall, to a minimally-resolved model of the far-field hydrodynamic interactions. Our experiments combine fluorescent labeling with particle tracking to trace the trajectories of individual particles in a dense suspension, and to measure their propulsion velocities. Previous computational studies [B. Sprinkle et al., J. Chem. Phys., 147, 244103, 2017] predicted that at sufficiently high densities a uniform suspension of microrollers separates into two layers, a slow monolayer right above the wall, and a fast layer on top of the bottom layer. Here we verify this prediction, showing good quantitative agreement between the bimodal distribution of particle velocities predicted by the lubrication-corrected Brownian Dynamics and those measured in the experiments. The computational method accurately predicts the rate at which particles are observed to switch between the slow and fast layers in the experiments. We also use our numerical method to demonstrate the important role that pairwise lubrication plays in motility-induced phase separation in dense monolayers of colloidal microrollers, as recently suggested for suspensions of Quincke rollers [D. Geyer et al., Physical Review X, 9(3), 031043, 2019]., Comment: Submitted to Soft Matter. See https://cims.nyu.edu/~donev/Publications.html for supplementary material animations

We introduce a Rigid-Body Fluctuating Immersed Boundary (RB-FIB) method to perform large-scale Brownian dynamics simulations of suspensions of rigid particles in fully confined domains, without any need to explicitly construct Green's functions or mobility operators. In the RB-FIB approach, discretized fluctuating Stokes equations are solved with prescribed boundary conditions in conjunction with a rigid-body immersed boundary method to discretize arbitrarily-shaped colloidal particles with no-slip or active-slip prescribed on their surface. We design a specialized Split--Euler--Maruyama temporal integrator that uses a combination of random finite differences to capture the stochastic drift appearing in the overdamped Langevin equation. The RB-FIB method presented in this work only solves mobility problems in each time step using a preconditioned iterative solver, and has a computational complexity that scales linearly in the number of particles and fluid grid cells. We demonstrate that the RB-FIB method correctly reproduces the Gibbs-Boltzmann equilibrium distribution, and use the method to examine the time correlation functions for two spheres tightly confined in a cuboid. We model a quasi--two-dimensional colloidal crystal confined in a narrow microchannel and hydrodynamically driven across a commensurate periodic substrate potential mimicking the effect of a corrugated wall. We observe partial and full depinning of the colloidal monolayer from the substrate potential above a certain wall speed, consistent with a transition from static to kinetic friction through propagating kink solitons. Unexpectedly, we find that particles nearest the boundaries of the domain are the first to be displaced, followed by particles in the middle of the domain.

At mesoscopic scales electrolyte solutions are modeled by the fluctuating generalized Poisson-Nernst-Planck (PNP) equations [J.-P. P\'eraud et al., Phys. Rev. F, 1(7):074103, 2016]. However, at length and time scales larger than the Debye scales, electrolytes are effectively electroneutral, and the charged-fluid PNP equations become too stiff to solve numerically. Here we formulate the isothermal incompressible equations of fluctuating hydrodynamics for reactive multispecies mixtures involving charged species in the electroneutral limit, and design a numerical algorithm to solve these equations. Our model does not assume a dilute electrolyte solution but rather treats all species on an equal footing, accounting for cross-diffusion and non-ideality using Maxwell-Stefan theory. By enforcing local electroneutrality as a constraint, we obtain an elliptic equation for the electric potential that replaces the Poisson equation in the fluctuating PNP equations. We develop a second-order midpoint predictor-corrector algorithm to solve either the charged-fluid or electroneutral equations with only a change of the elliptic solver. We use the electroneutral algorithm to study a gravitational fingering instability, triggered by thermal fluctuations, at an interface where an acid and base react to neutralize each other. Our results demonstrate that, because the four ions diffuse with very different coefficients, one must treat each ion as an individual species, and cannot treat the acid, base, and salt as neutral species. This emphasizes the differences between electrodiffusion and classical Fickian diffusion, even at electroneutral scales., Comment: To appear in Phys. Rev. Fluids

We apply fluctuating hydrodynamics to strong electrolyte mixtures to compute the concentration corrections for chemical potential, diffusivity, and conductivity. We show these corrections to be in agreement with the limiting laws of Debye, H\"uckel, and Onsager. We compute explicit corrections for a symmetric ternary mixture and find that the co-ion Maxwell-Stefan diffusion coefficients can be negative, in agreement with experimental findings., Comment: Submitted to Current Opinion in Electrochemistry

Fluctuating hydrodynamics (FHD) provides a framework for modeling microscopic fluctuations in a manner consistent with statistical mechanics and nonequilibrium thermodynamics. This paper presents an FHD formulation for isothermal reactive incompressible liquid mixtures with stochastic chemistry. Fluctuating multispecies mass diffusion is formulated using a Maxwell-Stefan description without assuming a dilute solution, and momentum dynamics is described by a stochastic Navier-Stokes equation for the fluid velocity. We consider a thermodynamically consistent generalization for the law of mass action for non-dilute mixtures and use it in the chemical master equation (CME) to model reactions as a Poisson process. The FHD approach provides remarkable computational efficiency over traditional reaction-diffusion master equation methods when the number of reactive molecules is large, while also retaining accuracy even when there are as few as ten reactive molecules per hydrodynamic cell. We present a numerical algorithm to solve the coupled FHD and CME equations and validate it on both equilibrium and nonequilibrium problems. We simulate a diffusively-driven gravitational instability in the presence of an acid-base neutralization reaction, starting from a perfectly flat interface. We demonstrate that the coupling between velocity and concentration fluctuations dominate the initial growth of the instability.

We develop a Split Reactive Brownian Dynamics (SRBD) algorithm for particle simulations of reaction-diffusion systems based on the Doi or volume reactivity model, in which pairs of particles react with a specified Poisson rate if they are closer than a chosen reactive distance. In our Doi model, we ensure that the microscopic reaction rules for various association and disassociation reactions are consistent with detailed balance (time reversibility) at thermodynamic equilibrium. The SRBD algorithm uses Strang splitting in time to separate reaction and diffusion, and solves both the diffusion-only and reaction-only subproblems exactly, even at high packing densities. To efficiently process reactions without uncontrolled approximations, SRBD employs an event-driven algorithm that processes reactions in a time-ordered sequence over the duration of the time step. A grid of cells with size larger than all of the reactive distances is used to schedule and process the reactions, but unlike traditional grid-based methods such as Reaction-Diffusion Master Equation (RDME) algorithms, the results of SRBD are statistically independent of the size of the grid used to accelerate the processing of reactions. We use the SRBD algorithm to compute the effective macroscopic reaction rate for both reaction- and diffusion-limited irreversible association in three dimensions. We also study long-time tails in the time correlation functions for reversible association at thermodynamic equilibrium. Finally, we compare different particle and continuum methods on a model exhibiting a Turing-like instability and pattern formation. We find that for models in which particles diffuse off lattice, such as the Doi model, reactions lead to a spurious enhancement of the effective diffusion coefficients., Comment: To appear in J. Chem. Phys

We present a fluctuating boundary integral method (FBIM) for overdamped Brownian Dynamics (BD) of two-dimensional periodic suspensions of rigid particles of complex shape immersed in a Stokes fluid. We develop a novel approach for generating Brownian displacements that arise in response to the thermal fluctuations in the fluid. Our approach relies on a first-kind boundary integral formulation of a mobility problem in which a random surface velocity is prescribed on the particle surface, with zero mean and covariance proportional to the Green's function for Stokes flow (Stokeslet). This approach yields an algorithm that scales linearly in the number of particles for both deterministic and stochastic dynamics, handles particles of complex shape, achieves high order of accuracy, and can be generalized to three dimensions and other boundary conditions. We show that Brownian displacements generated by our method obey the discrete fluctuation-dissipation balance relation (DFDB). Based on a recently-developed Positively Split Ewald method [A. M. Fiore, F. Balboa Usabiaga, A. Donev and J. W. Swan, J. Chem. Phys., 146, 124116, 2017], near-field contributions to the Brownian displacements are efficiently approximated by iterative methods in real space, while far-field contributions are rapidly generated by fast Fourier-space methods based on fluctuating hydrodynamics. FBIM provides the key ingredient for time integration of the overdamped Langevin equations for Brownian suspensions of rigid particles. We demonstrate that FBIM obeys DFDB by performing equilibrium BD simulations of suspensions of starfish-shaped bodies using a random finite difference temporal integrator., Comment: Submitted to J. Comp. Phys

We derive a minimal continuum model to investigate the hydrodynamic mechanism behind the fingering instability recently discovered in a suspension of microrollers near a floor [Driscoll et al. Nature Physics, 2016]. Our model, consisting of two continuous lines of rotlets, exhibits a linear instability driven only by hydrodynamics interactions, and reproduces the lengthscale selection observed in large scale particle simulations and in experiments. By adjusting only one parameter, the distance between the two lines, our dispersion relation exhibits quantitative agreement with the simulations and qualitative agreement with experimental measurements. Our linear stability analysis indicate that this instability is caused by the combination of the advective and transverse flows generated by the microrollers near a no-slip surface. Our simple model offers an interesting formalism to characterize other hydrodynamic instabilities that have not been yet well understood, such as size scale selection in suspensions of particles sedimenting adjacent to a wall, or the recently observed formations of traveling phonons in systems of confined driven particles.

In this letter we analyze the effects of an externally applied electric field on thermal fluctuations for a fluid containing charged species. We show in particular that the fluctuating Poisson-Nernst-Planck equations for charged multispecies diffusion coupled with the fluctuating fluid momentum equation, result in enhanced charge transport. Although this transport is advective in nature, it can macroscopically be represented as electrodiffusion with renormalized electric conductivity. We calculate the renormalized electric conductivity by deriving and integrating the structure factor coefficients of the fluctuating quantities and show that the renormalized electric conductivity and diffusion coefficients are consistent although they originate from different noise terms. In addition, the fluctuating hydrodynamics approach recovers the electrophoretic and relaxation corrections obtained by Debye-Huckel-Onsager theory, and provides a quantitative theory that predicts a non-zero cross-diffusion Maxwell-Stefan coefficient that agrees well with experimental measurements. Finally, we show that strong applied electric fields result in anisotropically enhanced velocity fluctuations and reduced fluctuations of salt concentrations., Comment: 12 pages, 1 figure

We combine experiments, large scale simulations and continuum models to study the emergence of coherent structures in a suspension of magnetically driven microrollers sedimented near a floor. Collective hydrodynamic effects are predominant in this system, leading to strong density-velocity coupling. We characterize a uniform suspension and show that density waves propagate freely in all directions in a dispersive fashion. When sharp density gradients are introduced in the suspension, we observe the formation of a shock. Unlike Burgers' shock-like structures observed in other active and driven confined hydrodynamic systems, the shock front in our system has a well-defined finite width and moves rapidly compared to the mean suspension velocity. We introduce a continuum model demonstrating that the finite width of the front is due to far-field nonlocal hydrodynamic interactions and governed by a geometric parameter: the average particle height above the floor.

The Immersed Boundary (IB) method is a mathematical framework for constructing robust numerical methods to study fluid-structure interaction in problems involving an elastic structure immersed in a viscous fluid. The IB formulation uses an Eulerian representation of the fluid and a Lagrangian representation of the structure. The Lagrangian and Eulerian frames are coupled by integral transforms with delta function kernels. The discretized IB equations use approximations to these transforms with regularized delta function kernels to interpolate the fluid velocity to the structure, and to spread structural forces to the fluid. It is well-known that the conventional IB method can suffer from poor volume conservation since the interpolated Lagrangian velocity field is not generally divergence-free, and so this can cause spurious volume changes. In practice, the lack of volume conservation is especially pronounced for cases where there are large pressure differences across thin structural boundaries. The aim of this paper is to greatly reduce the volume error of the IB method by introducing velocity-interpolation and force-spreading schemes with the properties that the interpolated velocity field in which the structure moves is at least C1 and satisfies a continuous divergence-free condition, and that the force-spreading operator is the adjoint of the velocity-interpolation operator. We confirm through numerical experiments in two and three spatial dimensions that this new IB method is able to achieve substantial improvement in volume conservation compared to other existing IB methods, at the expense of a modest increase in the computational cost. Further, the new method provides smoother Lagrangian forces (tractions) than traditional IB methods. The method presented here is restricted to periodic computational domains. Its generalization to non-periodic domains is important future work., Comment: 49 pages, 13 figures

We develop numerical methods for reaction-diffusion systems based on the equations of fluctuating hydrodynamics (FHD). While the FHD formulation is formally described by stochastic partial differential equations (SPDEs), it becomes similar to the reaction-diffusion master equation (RDME) description when those SPDEs are spatially discretized and reactions are modeled as a source term having Poisson fluctuations. However, unlike the RDME, the FHD description naturally extends from the regime where fluctuations are strong, i.e., each hydrodynamic cell has few (reactive) molecules, to regimes with moderate or weak fluctuations, and ultimately to the deterministic limit. By treating diffusion implicitly, we avoid the severe restriction on time step size that limits all methods based on explicit treatments of diffusion, and construct numerical methods that are more efficient than RDME methods, without compromising accuracy. We find that an implicit midpoint tau leaping scheme attains second-order weak accuracy for the linearized FHD, and gives an accurate and stable structure factor for a time step size an order of magnitude larger than the hopping time scale of diffusing molecules. We study the numerical accuracy of our methods for the Schl\"ogl reaction-diffusion model both in and out of thermodynamic equilibrium. We demonstrate and quantify the importance of thermodynamic fluctuations to the formation of a two-dimensional Turing-like pattern, and examine the effect of fluctuations on three-dimensional chemical front propagation. By comparing FHD simulations to deterministic reaction-diffusion simulations, we show that fluctuations accelerate pattern formation in spatially homogeneous systems, and lead to a qualitatively-different disordered pattern behind a traveling wave., Comment: To appear in J. Chem. Phys., 2017

We develop efficient numerical methods for performing many-body Brownian dynamics simulations of a recently-observed fingering instability in an active suspension of colloidal rollers sedimented above a wall [M. Driscoll, B. Delmotte, M. Youssef, S. Sacanna, A. Donev and P. Chaikin, Nature Physics, 2016, doi:10.1038/nphys3970]. We present a stochastic Adams-Bashforth integrator for the equations of Brownian dynamics, which has the same cost as but is more accurate than the widely-used Euler-Maruyama scheme, and uses a random finite difference to capture the stochastic drift proportional to the divergence of the configuration-dependent mobility matrix. We generate the Brownian increments using a Krylov method, and show that for particles confined to remain in the vicinity of a no-slip wall by gravity or active flows the number of iterations is independent of the number of particles. Our numerical experiments with active rollers show that the thermal fluctuations set the characteristic height of the colloids above the wall, both in the initial condition and the subsequent evolution dominated by active flows. The characteristic height in turn controls the timescale and wavelength for the development of the fingering instability., Comment: Revised and resubmitted to J. Chem. Phys

Condensation of objects into stable clusters occurs naturally in equilibrium and driven systems. It is commonly held that potential interactions, depletion forces, or sensing are the only mechanisms which can create long-lived compact structures. Here we show that persistent motile structures can form spontaneously from hydrodynamic interactions alone with no sensing or potential interactions. We study this structure formation in a system of colloidal rollers suspended and translating above a floor, using both experiments and large-scale 3D simulations. In this system, clusters originate from a previously unreported fingering instability, where fingers pinch off from an unstable front to form autonomous "critters", whose size is selected by the height of the particles above the floor. These critters are a stable state of the system, move much faster than individual particles, and quickly respond to a changing drive. With speed and direction set by a rotating magnetic field, these active structures offer interesting possibilities for guided transport, flow generation, and mixing at the microscale.

We formulate and study computationally the low Mach number fluctuating hydrodynamic equations for electrolyte solutions. We are interested in studying transport in mixtures of charged species at the mesoscale, down to scales below the Debye length, where thermal fluctuations have a significant impact on the dynamics. Continuing our previous work on fluctuating hydrodynamics of multicomponent mixtures of incompressible isothermal miscible liquids (A. Donev, et al., Physics of Fluids, 27, 3, 2015), we now include the effect of charged species using a quasielectrostatic approximation. Localized charges create an electric field, which in turn provides additional forcing in the mass and momentum equations. Our low Mach number formulation eliminates sound waves from the fully compressible formulation and leads to a more computationally efficient quasi-incompressible formulation. We demonstrate our ability to model saltwater (NaCl) solutions in both equilibrium and nonequilibrium settings. We show that our algorithm is second-order in the deterministic setting, and for length scales much greater than the Debye length gives results consistent with an electroneutral/ambipolar approximation. In the stochastic setting, our model captures the predicted dynamics of equilibrium and nonequilibrium fluctuations. We also identify and model an instability that appears when diffusive mixing occurs in the presence of an applied electric field., Comment: 37 pages, 5 figures

We derive a coarse-grained description of the dynamics of a nanoparticle immersed in an isothermal simple fluid by performing a systematic coarse graining of the underlying microscopic dynamics. As coarse-grained or relevant variables we select the position of the nanoparticle and the \emph{total} mass and momentum density field of the fluid, which are locally conserved slow variables because they are defined to include the contribution of the nanoparticle. The theory of coarse graining based on the Zwanzing projection operator leads us to a system of stochastic \emph{ordinary} differential equations (SODEs) that are closed in the relevant variables. We demonstrate that our discrete coarse-grained equations are consistent with a Petrov-Galerkin finite-element discretization of a system of formal stochastic \emph{partial} differential equations (SPDEs) which resemble previously-used phenomenological models based on fluctuating hydrodynamics. Under suitable approximations we obtain \emph{closed} approximations of the coarse-grained dynamics in a manner which gives them a clear physical interpretation, and provides \emph{explicit} microscopic expressions for all of the coefficients appearing in the closure. Our work leads to a model for dilute nanocolloidal suspensions that can be simulated effectively using feasibly short molecular dynamics simulations as input to a FEM fluctuating hydrodynamic solver., Comment: Submitted to J. Chem. Phys

We introduce numerical methods for simulating the diffusive motion of rigid bodies of arbitrary shape immersed in a viscous fluid. We parameterize the orientation of the bodies using normalized quaternions, which are numerically robust, space efficient, and easy to accumulate. We construct a system of overdamped Langevin equations in the quaternion representation that accounts for hydrodynamic effects, preserves the unit-norm constraint on the quaternion, and is time reversible with respect to the Gibbs-Boltzmann distribution at equilibrium. We introduce two schemes for temporal integration of the overdamped Langevin equations of motion, one based on the Fixman midpoint method and the other based on a random finite difference approach, both of which ensure the correct stochastic drift term is captured in a computationally efficient way. We study several examples of rigid colloidal particles diffusing near a no-slip boundary, and demonstrate the importance of the choice of tracking point on the measured translational mean square displacement (MSD). We examine the average short-time as well as the long-time quasi-two-dimensional diffusion coefficient of a rigid particle sedimented near a bottom wall due to gravity. For several particle shapes we find a choice of tracking point that makes the MSD essentially linear with time, allowing us to estimate the long-time diffusion coefficient efficiently using a Monte Carlo method. However, in general such a special choice of tracking point does not exist, and numerical techniques for simulating long trajectories, such as the ones we introduce here, are necessary to study diffusion on long timescales., Comment: To appear in J. Chem. Phys, revised somewhat in v2

Diffusion processes are widespread in biological and chemical systems, where they play a fundamental role in the exchange of substances at the cellular level and in determining the rate of chemical reactions. Recently, the classical picture that portrays diffusion as random uncorrelated motion of molecules has been revised, when it was shown that giant non-equilibrium fluctuations develop during diffusion processes. Under microgravity conditions and at steady-state, non-equilibrium fluctuations exhibit scale invariance and their size is only limited by the boundaries of the system. In this work, we investigate the onset of non-equilibrium concentration fluctuations induced by thermophoretic diffusion in microgravity, a regime not accessible to analytical calculations but of great relevance for the understanding of several natural and technological processes. A combination of state of the art simulations and experiments allows us to attain a fully quantitative description of the development of fluctuations during transient diffusion in microgravity. Both experiments and simulations show that during the onset the fluctuations exhibit scale invariance at large wave vectors. In a broader range of wave vectors simulations predict a spinodal-like growth of fluctuations, where the amplitude and length-scale of the dominant mode are determined by the thickness of the diffuse layer., Comment: To appear in Scientific Reports

Fluctuations in a fluid are strongly affected by the presence of a macroscopic gradient making them long-ranged and enhancing their amplitude. While small-scale fluctuations exhibit diffusive lifetimes, larger-scale fluctuations live shorter because of gravity, as theoretically and experimentally well-known. We explore here fluctuations of even larger size, comparable to the extent of the system in the direction of the gradient, and find experimental evidence of a dramatic slowing-down in their dynamics. We recover diffusive behaviour for these strongly-confined fluctuations, but with a diffusion coefficient that depends on the solutal Rayleigh number. Results from dynamic shadowgraph experiments are complemented by theoretical calculations and numerical simulations based on fluctuating hydrodynamics, and excellent agreement is found. The study of the dynamics of non-equilibrium fluctuations allows to probe and measure the competition of physical processes such as diffusion, buoyancy and confinement., Comment: Includes see Supplementary Material. Submitted to PRL

We present a finite element discretization of a non-linear diffusion equation used in the field of critical phenomena and, more recently, in the context of Dynamic Density Functional Theory. The discretized equation preserves the structure of the continuum equation. Specifically, it conserves the total number of particles and fulfills an H-theorem as the original partial differential equation. Guided by the Theory of Coarse-Graining, we discuss the inclusion of thermal fluctuations in the non-linear diffusion equation. This sheds light on the meaning of such a fluctuating hydrodynamics equation and to the limitations of the approximations involved. The methodology proposed for the introduction of thermal fluctuations in finite element methods is general and valid for both regular and irregular grids in arbitrary dimensions. We focus here on simulations of the Ginzburg-Landau free energy functional using both regular and irregular 1D grids. Convergence of the numerical results is obtained for the static and dynamic structure factors as the resolution of the grid is increased., Comment: Updated version after review. Accepted for publication in the Journal of Chemical Physics

Fluctuating hydrodynamics provides a model for fluids at mesoscopic scales where thermal fluctuations can have a significant impact on the behavior of the system. Here we investigate a model for fluctuating hydrodynamics of a single component, multiphase flow in the neighborhood of the critical point. The system is modeled using a compressible flow formulation with a van der Waals equation of state, incorporating a Korteweg stress term to treat interfacial tension. We present a numerical algorithm for modeling this system based on an extension of algorithms developed for fluctuating hydrodynamics for ideal fluids. The scheme is validated by comparison of measured structure factors and capillary wave spectra with equilibrium theory. We also present several non-equilibrium examples to illustrate the capability of the algorithm to model multi-phase fluid phenomena in a neighborhood of the critical point. These examples include a study of the impact of fluctuations on the spinodal decomposition following a rapid quench, as well as the piston effect in a cavity with supercooled walls. The conclusion in both cases is that thermal fluctuations affect the size and growth of the domains in off-critical quenches.

In this paper we discuss the formulation of the fuctuating Navier-Stokes (FNS) equations for multi-species, non-reactive fluids. In particular, we establish a form suitable for numerical solution of the resulting stochastic partial differential equations. An accurate and efficient numerical scheme, based on our previous methods for single species and binary mixtures, is presented and tested at equilibrium as well as for a variety of non-equilibrium problems. These include the study of giant nonequilibrium concentration fluctuations in a ternary mixture in the presence of a diffusion barrier, the triggering of a Rayleigh-Taylor instability by diffusion in a four-species mixture, as well as reverse diffusion in a ternary mixture. Good agreement with theory and experiment demonstrates that the formulation is robust and can serve as a useful tool in the study of thermal fluctuations for multi-species fluids. The extension to include chemical reactions will be treated in a sequel paper.

We present explicit methods for simulating diffusions whose generator is self-adjoint with respect to a known (but possibly not normalizable) density. These methods exploit this property and combine an optimized Runge-Kutta algorithm with a Metropolis-Hastings Monte-Carlo scheme. The resulting numerical integration scheme is shown to be weakly accurate at finite noise and to gain higher order accuracy in the small noise limit. It also permits to avoid computing explicitly certain terms in the equation, such as the divergence of the mobility tensor, which can be tedious to calculate. Finally, the scheme is shown to be ergodic with respect to the exact equilibrium probability distribution of the diffusion when it exists. These results are illustrated on several examples including a Brownian dynamics simulation of DNA in a solvent. In this example, the proposed scheme is able to accurately compute dynamics at time step sizes that are an order of magnitude (or more) larger than those permitted with commonly used explicit predictor-corrector schemes., Comment: 54 pages, 8 figures, To appear in MMS

A novel stochastic fluid model is proposed with non-ideal structure factor consistent with compressibility, and adjustable transport coefficients. This Stochastic Hard Sphere Dynamics (SHSD) algorithm is a modification of the Direct Simulation Monte Carlo (DSMC) algorithm and has several computational advantages over event-driven hard-sphere molecular dynamics. Surprisingly, SHSD results in an equation of state and pair correlation function identical to that of a deterministic Hamiltonian system of penetrable spheres interacting with linear core pair potentials. The fluctuating hydrodynamic behavior of the SHSD fluid is verified for the Brownian motion of a nano-particle suspended in a compressible solvent., Comment: This work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 (LLNL-JRNL-401745). To appear in Phys. Rev. Lett. 2008

A novel Stochastic Event-Driven Molecular Dynamics (SEDMD) algorithm is developed for the simulation of polymer chains suspended in a solvent. The polymers are represented as chains of hard spheres tethered by square wells and interact with the solvent particles with hard core potentials. The algorithm uses Event-Driven Molecular Dynamics (EDMD) for the simulation of the polymer chain and the interactions between the chain beads and the surrounding solvent particles. The interactions between the solvent particles themselves are not treated deterministically as in event-driven algorithms, rather, the momentum and energy exchange in the solvent is determined stochastically using the Direct Simulation Monte Carlo (DSMC) method. The coupling between the solvent and the solute is consistently represented at the particle level, however, unlike full MD simulations of both the solvent and the solute, the spatial structure of the solvent is ignored. The algorithm is described in detail and applied to the study of the dynamics of a polymer chain tethered to a hard wall subjected to uniform shear. The algorithm closely reproduces full MD simulations with two orders of magnitude greater efficiency. Results do not confirm the existence of periodic (cycling) motion of the polymer chain., Comment: This work was performed under the auspices of the U.S. Department of Energy by the University of California Lawrence Livermore National Laboratory under Contract No. W-7405-Eng-48 (UCRL-ABS-233205)

We present, in a unifying way, the main components of three asynchronous event-driven algorithms for simulating physical systems of interacting particles. The first example, hard-particle molecular dynamics, is well-known. We also present a recently-developed diffusion kinetic Monte Carlo algorithm, as well as a novel stochastic molecular-dynamics algorithm that builds on the Direct Simulation Monte Carlo. We explain how to effectively combine asynchronous event-driven with classical time-driven or with synchronous event-driven handling. Finally, we discuss some promises and challenges for event-driven simulation of realistic physical systems., Comment: To appear in Simulation: Transactions of the Society for Modeling and Simulation International, 2008

In the first part of a series of two papers, we present in considerable detail a collision-driven molecular dynamics algorithm for a system of nonspherical particles, within a parallelepiped simulation domain, under both periodic or hard-wall boundary conditions. The algorithm extends previous event-driven molecular dynamics algorithms for spheres. We present a novel partial-update near-neighbor list (NNL) algorithm that is superior to previous algorithms at high densities, without compromising the correctness of the algorithm. This efficiency of the algorithm is further increased for systems of very aspherical particles by using bounding sphere complexes (BSC). In the second part of this series of papers we apply the algorithm presented in the first part of this series of papers to systems of hard ellipses and ellipsoids. The theoretical machinery needed to treat such particles, including the overlap potentials, is developed in full detail. We describe an algorithm for predicting the time of collision for two moving ellipses or ellipsoids. We present performance results for our implementation of the algorithm. The practical utility of the algorithm is demonstrated by presenting several interesting physical applications, including the generation of jammed packings inside spherical containers, the study of contact force chains in jammed packings, and melting the densest-known equilibrium crystals of prolate spheroids., Comment: Two-paper series, 45 pages per paper, 90 pages total. To be submitted to JCP