Particle dynamics modeling methods for colloid suspensions

We present a review and critique of several methods for the simulation of the dynamics of colloidal suspensions at the mesoscale. We focus particularly on simulation techniques for hydrodynamic interactions, including implicit solvents (Fast Lubrication Dynamics, an approximation to Stokesian Dynamics) and explicit/particle-based solvents (Multi-Particle Collision Dynamics and Dissipative Particle Dynamics). Several variants of each method are compared quantitatively for the canonical system of monodisperse hard spheres, with a particular focus on diffusion characteristics, as well as shear rheology and microstructure. In all cases, we attempt to match the relevant properties of a well-characterized solvent, which turns out to be challenging for the explicit solvent models. Reasonable quantitative agreement is observed among all methods, but overall the Fast Lubrication Dynamics technique shows the best accuracy and performance. We also devote significant discussion to the extension of these methods to more complex situations of interest in industrial applications, including models for non-Newtonian solvent rheology, non-spherical particles, drying and curing of solvent and flows in complex geometries. This work identifies research challenges and motivates future efforts to develop techniques for quantitative, predictive simulations of industrially relevant colloidal suspension processes.

academic perspective and for their relevance in technological applications. Although the study of colloidal suspensions extends more broadly to different phases of the dispersed and dispersion media (e.g. foams, emulsions, gels, aerosols) [94,189], we focus exclusively on solid particles dispersed in a liquid medium (sols). For the purposes of this work, we define the colloidal regime based on a length scale that is much larger than molecular dimensions, but small enough that thermal fluctuations as well as colloidal forces (i.e. van der Waals dispersion and long-range electrostatic forces) are significant. For typical systems, this translates to effective particle diameters on the order of several tens of nanometers to several microns (∼0.1-2 μm). Traditional application areas for such colloidal suspensions include sedimentation of clays [54,166], oil recovery and transport [143,213], and rheological properties of paints, cosmetics, adhesives, food products, etc. [44,53]. More recently, there has been renewed interest in particle suspensions due to their presence in biological systems (e.g. red blood cells and platelets in blood [43,84,148,177], microorganisms in biological suspensions [61], proteins and nucleic acids in aqueous environments [108], drug encapsulation and delivery [157]) as well as their applications in nanotechnology and nanomanufacturing processes (e.g. nanoparticle-reinforced composites [76,215], battery electrodes [236], control of optoelectronic properties by specific arrangements of colloidal particles [195], sol-gel synthesis [28]).
From an applications perspective, the goal of studying colloidal suspensions is to understand, predict and control suspension rheological properties and particle microstructure. In applications where the desired product is a liquid suspension (e.g. cosmetics, adhesives, paints, food products), rheological properties are directly relevant to product performance; in cases where the solvent is cured or removed, the rheological properties of the suspension are critical to manufacturing processing conditions, whereas the final (dry) microstructure largely governs product properties and performance. For example, in injection-molded parts, particulate composites and thin film coatings, the spatial and orientational distribution of particles can drastically affect product properties.
However, both rheology and microstructure formation are complex, non-equilibrium phenomena. A more fundamental process that governs both of these is colloid diffusion. Although it may not be directly relevant in some applications, we argue that diffusion, which is a dynamic, equilibrium process, fundamentally underlies rheological properties and microstructure formation, and must be well-understood before these more complex phenomena can be addressed. Furthermore, diffusivity is in principle a readily measurable quantity, and diffusion is sufficiently complex to reflect many of the relevant underlying physics (hydrodynamic interactions, colloidal forces). Colloidal diffusion processes have also been studied extensively by theory, experiments and simulations. We therefore identify diffusion, rheological properties and microstructure formation as the three main properties of interest in the study of colloidal suspensions.
One of the difficulties in studying colloidal systems lies in obtaining accurate information on disparate time and length scales, which is essential to understanding macroscale behavior. Common experimental techniques include various forms of radiation scattering (e.g. X-ray, neutron, laser), confocal microscopy, fluorescence recovery after photobleaching (FRAP), pulsed-field-gradient nuclear magnetic resonance (PFG NMR), and a host of conventional techniques for probing rheology and imaging microstructures. With regard to dynamics, light scattering techniques provide extensive information on the spatial and temporal correlations between particles, typically in the form of structure factors and intermediate scattering functions [14,56,151,183]. However, with increasing complexity (e.g. particle aggregation, polydispersity, complex particle shapes and flow geometries), interpreting these quantities becomes more difficult, and significant limitations exist with regard to particle sizes, volume fractions and accessible time scales. Confocal microscopy techniques paired with image recognition algorithms have emerged as an important tool for directly probing the dynamics and structural evolution of colloidal systems [39,46,55,230], but also have limitations with regard to particle size, time scales and various limiting requirements for the optical properties of the particles and solvent (e.g. equal refractive indices). PFG NMR [158,209], FRAP [20,102] and rheology techniques [77,155,188,189] provide key information about average dynamical properties such as diffusion coefficients and frequency-dependent rheological data, but cannot directly probe dynamics at the individual particle scale, and the interpretation of microstructure evolution is often not straightforward. Computer simulations at the mesoscale, where colloidal particles are treated as dis-crete elements much larger than the atomic scale, have therefore emerged as an important additional tool for the study of colloidal suspensions, both to aid in the interpretation of experimental data and in some cases as standalone predictive techniques.
The focus of this work is to discuss and evaluate several computer simulation methods relevant to colloidal suspensions. We do not attempt a comprehensive literature review of all such techniques, as the field is quite large; instead, we focus on several methods that are primarily off-lattice and particle-based, where the colloidal particles are always treated explicitly (sometimes referred to as a Lagrangian approach), while the solvent is treated either as a separate set of particles (explicit) or incorporated in the colloid-colloid interaction (implicit). The methods we have selected are Fast Lubrication Dynamics (FLD) [31,124,126], which is a recently-developed expedient approximation to Stokesian dynamics (SD) [25][26][27]; multi-particle collision dynamics (MPCD) [92,164,169] originally known as stochastic rotation dynamics (SRD) [144]; and dissipative particle dynamics (DPD) [87,95]. We largely exclude several important classes of colloid suspension simulation techniques, including lattice-based treatments of the solvent (in particular, Lattice-Boltzmann techniques [1,37,57,128]; field-based or Eulerian treatments of the particles, where a spatially-varying continuous density variable is used rather than explicit particles to represent the solid phase [64]; and a host of continuum Navier-Stokes-based techniques for treatment of the solvent, including various finite element, finite difference, volume-of-fluid, boundary element and spectral methods [66,78,80,83,96,134,163]). That said, we do cite important work in these areas in Sect. 5 of this paper, as often they offer the only route to address challenging applications. We also exclude particle-based continuum treatments of the solvent (sometimes referred to as "meshless methods" [13,161]), such as smoothed particle hydrodynamics (SPH) [82,159], the element-free Galerkin method (EFGM) [12,73], the reproducing kernel particle method (RKPM) [137,138] or the finite point method (FPM) [167]. To our knowledge, these latter methods have not been applied to the study of colloidal suspensions, and their potential in this area remains largely unexplored. All of the aforementioned techniques have proved extremely useful in many cases, and we discuss them briefly where appropriate.
We note that many excellent reviews of computer simulation techniques for colloidal suspensions similar to those on which we are focusing already exist in the literature. Barnes et al. [9] provided an early survey of particle-based simulation methods for dense suspension rheology, but did not carry out a quantitative comparison of the various methods, and several currently popular techniques (DPD, MPCD) had not been developed at the time of their work. Brady and Bossis [27] and Foss and Brady [75] have summarized the development and extensive applications of the Stokesian Dynamics method, which has been highly successful for a broad range of problems. Harting et al. [91] reviewed the use of MPCD and Lattice-Boltzmann methods for particle suspensions, and provided useful guidelines as to which method is more appropriate for various physical situations. The review by Dunweg and Ladd [57] of Lattice-Boltzmann methods for colloidal suspensions also gives an excellent overview of DPD, MPCD and SD-based methods and the underlying theory. Van der Sman has presented an excellent survey of methods that include Lattice-Boltzmann, a variety of continuum-based treatments of the solvent as well as DPD and MPCD for suspension flows in confined geometries [207]. His work provides qualitative insights into matching these diverse methods with various length scales and flow regimes, and suggests that the Lattice-Boltzmann method may be the most versatile for complex geometries over multiple length scales. More recently, Dickinson [52] has reviewed the state of the art for a number of particle-based simulation methods and the insights they provide on colloid aggregation and phase behavior.
However, the reviews listed above rarely report quantitative comparisons of simulation methods for the same system in a manner that allows for direct critical evaluation of their accuracy. Several authors have recently presented rigorous quantitative comparisons of some of these methods: Batôt et al. [11] have compared colloid diffusion and conductivity for neutral and charged particles using MPCD and Brownian dynamics (BD). They employed two different collisional coupling schemes for MPCD, and BD with no hydrodynamic interactions as well as using a Felderhof hydrodynamic tensor [70]. Tomilov et al. [217] quantified various colloid aggregation properties using BD, BD with a Yamakawa-Rotne-Prager hydrodynamic tensor [234] and MPCD. Schlauch et al. [193] compared finite element, Lattice-Boltzmann, and Stokesian dynamics methods for their ability to reproduce forces and torques on assemblies of stationary colloid particles, but did not investigate dynamic suspensions. In previous work [197], we reported a comparison of MPCD, FLD and DPD in the context of diffusion and rheological properties for a system of colloidal polystyrene particles in aqueous solution. The simulation methods were tested briefly against a system of hard-sphere colloids, but the emphasis was on their ability to reproduce the experimentally measured properties of the polystyrene-water system. In this work, we delve into greater detail for the hard-sphere colloid system with an emphasis on colloid diffusion characteristics, and explore additional methodological details for all of these simulation techniques. All the methods that we discuss have been incorporated into the same parallel molecular dynamics software package (LAMMPS [176]), which will allow for an even more direct comparison, including computational performance. Furthermore, we devote substantial discussion in Sect. 5 to the feasibility of extending these techniques for predictive simulations of more complex physics and realistic manufacturing processes. Overall, it is hoped that this work will help identify areas of deficiency for the various techniques discussed and motivate additional method development efforts.
A high-level summary of the key findings of this work is presented in Table. 1. The table shows that while implicit solvent methods are more accurate for the simple systems addressed in this paper, explicit methods offer far greater promise for more advanced applications. The underpinning reasons for this will become clearer in the sections that follow. Specifically, in Sect. 5 we discuss the various extensions of these techniques to more complex physical situations. We include this table here as motivation for what follows, as well as a convenient summary of outstanding research challenges in the field.
The remainder of the paper is organized as follows: Sect. 2 presents a historic perspective on key developments in colloid suspension theory and an overview of the dominant physics involved in colloidal suspension modeling. Section 3 discusses the details of the different simulation techniques that we focus on in this work (FLD, MPCD and DPD), including a section devoted to the efficient implementation of simulation algorithms for colloidal systems. In Sect. 4, we present simulation results from our work as well as from the literature. We emphasize diffusion of hard-sphere colloids as a key metric for the accuracy of various simulation methods, but also discuss predictions of rheological behavior and equilibrium colloid particle microstructure. In all cases, we provide an evaluation of the different simulation methods, including accuracy, flexibility and computational speed. Section 5 is devoted to advanced capabilities standing in the way of extension of these methods to complex applications. Finally, we conclude with a summary of our findings, some of the outstanding challenges, and future directions for particle-based simulations of colloids.

Historical perspective
Seminal work in the theory of colloidal suspensions dates back to investigations of single-particle Brownian motion by Einstein [62] and Smoluchowski [227], and subsequent experimental validation by Perrin [171], all of which were instrumental in confirming the atomic nature of matter. Subsequent theoretical studies on the equilibrium properties of colloidal suspensions by Derjaguin and Landau [49] as well as Verwey and Overbeek [224] focused on understanding the origins of interactions between colloidal particles. This resulted in the Derjaguin-Landau-Verwey-Overbeek Table 1 Summary of the current capabilities of simulation methods treated in this work. Coloring indicates the relative state of development for various capabilities (green-most mature; yellow-moderately mature; red-least mature). (Color table online) (DLVO) theory, which describes forces of molecular origin between bodies immersed in an electrolyte, and can be used to quantify several key features of colloid aggregation. With regard to dynamic properties, extensive theoretical developments towards Fokker-Planck formulations of single-and multi-particle colloidal systems appear in several sources [34,50,147]. While the direct use of such formulations is generally unwieldy, they provide a rigorous starting point and a justification for the Langevin and generalized Langevin equations that form the basis of many current simulation algorithms. Additionally, several approximations for colloid diffusion coefficients and suspension viscosity have been derived theoretically for simple systems in the dilute limit [10]. Unfortunately, for most practical applications of interest, the complexity of realistic colloidal suspensions precludes the direct use of a purely theoretical treatment.
Computer-aided simulation has greatly accelerated the understanding of colloidal systems. The typical computational framework consists of an N -body Newtonian dynamics solver to treat the colloidal particles (i.e. a discrete element model, or DEM framework), coupled with a solver for the hydrodynamic effects of the suspending fluid. Origi-nal algorithms for DEM simulations date back to molecular dynamics (MD) simulations of simple Lennard-Jones fluids in the late 1950s [2], but a fully atomistic treatment of a colloidal suspension still remains largely intractable due to the vast separation in length and time scales between solvent molecules and colloidal particles. Typically, the solvent degrees of freedom are coarse-grained in some fashion. Seminal work by Ermak and McCammon resulted in a Brownian dynamics algorithm for the simulation of a collection of spherical particles, with various approximations for hydrodynamic interactions [65]. A significant refinement of the Ermak and McCammon algorithm came in the form of Stokesian dynamics (SD) [25][26][27]59]. SD includes a near-field correction to the interparticle hydrodynamic interaction as well as improvements to the far-field interaction by accounting for torques and stresslets. Both the Ermak-McCammon algorithm and Stokesian dynamics algorithms treat hydrodynamics based on multipole expansions to solve the Stokes equation in the presence of moving spheres. Despite often prohibitive computational costs, these methods have proved to be remarkably accurate and extremely useful, and more recent accelerated versions [8,154,200] as well as simplified versions such as Fast Lubrication Dynamics (FLD) [126] show even greater promise. However, due to their underlying assumptions, they are often limited in terms of flexibility to different geometries (particle shapes as well as overall domain geometry) and complicating physical phenomena (e.g. complex solvent rheology, drying-see Table. 1).
In parallel, an increasing number of lattice-based (e.g. lattice gas automata and Lattice-Boltzmann; various gridbased Navier-Stokes solvers) as well as off-lattice particlebased treatments of the solvent have emerged. Some of the significant developments aimed at colloidal suspensions came in the form of Lattice-Boltzmann algorithms which include thermal fluctuations [130] and the development of coarse-grained particle-based solvents methods like dissipative particle dynamics (DPD) [95] and smoothed particle hydrodynamics (SPH) [82,140]. More recently, a hybrid lattice-off-lattice algorithm known as stochastic rotation dynamics (SRD) [144] or Multi-Particle Collision Dynamics (MPCD) [92,164] has gained traction in the soft matter simulation community. The SRD algorithm was derived from techniques developed for rarefied gas flow simulations, e.g. direct simulation Monte Carlo (DSMC) [18]. The DPD and MPCD methods are discussed in greater detail in Sect. 3, following a description of the key physics involved in modeling colloidal suspensions.

Governing physical principles
We describe a suspension of colloidal particles as a collection of N interacting, moving particles in carrier liquid, or "solvent". In three dimensions, each colloidal particle i has in general three translational degrees of freedom r i and three orientational degrees of freedom ω i . We group these into the 6N -dimensional vector r = {r 1 , r 2 , . . . , r N , ω 1 , ω 2 , . . . , ω N }. Similarly, the translational and angular velocities of all particles are lumped into the vector U, the forces and torques into the vector F, and the masses and moments of inertia into the matrix M.
In the general framework of a collection of N colloidal particles submerged in a bath of N s solvent particles, the equations of motion for the colloidal particles can be written as: Here, F C is the resultant (conservative) force acting on the particles due to interactions with other colloidal particles, and F S is the resultant force due to interactions with solvent particles. When the solvent is treated atomistically, the nature and form of the two interactions are often similar, and an additional N s equations similar to (1) are required; when the solvent particles are coarse-grained representations of actual solvent molecules, F S is given by a modified pairwise interaction or a suitable solvent-colloid collision rule. The choice of these rules and interactions depends on the method. We will return to the details of the various methods in subsequent sections.
An alternative approach to the simulation of additional solvent particles relies on a continuum treatment of the solvent. As a result, the solvent degrees of freedom can be discarded, but the forces experienced by the colloid particles must be modified to reproduce the key effects of the solvent. This typically assumes that the momentum relaxation timescale in the solvent is much smaller than the timescale of colloid diffusion (i.e. steady flow). The motion of a given colloidal particle induces a flow field in the solvent, which in turn affects its own motion (drag) as well as the motion of the remaining colloidal particles (hydrodynamic interactions). Collectively, these effects are lumped into the hydrodynamic forces, F H . In addition, colloidal particles experience a large number of collisions with solvent molecules, which give rise to Brownian motion. We denote the resulting forces as F B , the total Brownian force. The equations of motion for such a model are therefore given by: We refer to all such models of colloidal suspensions as implicit solvent models. While the term "implicit solvent" often refers to models that also account for the thermodynamic effects of solvents on conservative inter-particle interactions [187] (e.g. electrostatic screening, hydrophobic effects, etc.), we use the term to refer exclusively to the hydrodynamic and dissipative effects of the solvent. The aforementioned thermodynamic properties of the solvent are accounted for in the conservative interparticle force term F C . In both implicit and explicit solvent models, F C is a conservative force that depends only on particle positions (and orientations), while the hydrodynamic force F H in implicit solvent models depends additionally on particle velocities, and the Brownian force is constrained to satisfy the fluctuationdissipation theorem.
In the remainder of this section, we provide a discussion of the general physical principles underlying each of these types of forces. The details pertaining to the treatment and implementation of the solvent models are largely deferred to Sect. 3.

Hydrodynamic interactions
Hydrodynamic interactions and Brownian motion of colloidal particles arise naturally when using explicit representations of the solvent, so long as several relatively simple conditions are met. First, the solvent dynamics should conserve mass, linear and angular momentum and energy on the time and length scales relevant to the motion of colloidal particles. When a constant temperature must be maintained in the solvent, the thermostatting algorithm must be carefully designed so as to minimize interference with the flow characteristics. Second, the coupling between colloidal particles and solvent must be done in a manner consistent with particle-scale flow characteristics; in the canonical case of a hard-sphere colloidal suspension, this typically entails enforcing no-slip boundary conditions at the particle surface. Several other features are additionally desirable for explicit solvents, such as Galilean invariance, incompressibility, realistic solvent rheological and thermodynamic properties and realistic Schmidt number (ratio of momentum to mass diffusivity). In many cases, the solvent coarse-graining process compromises on several of these features in order to keep the methods computationally expedient (cf. Sect. 3).
Implicit solvent methods augment the particle forces with hydrodynamic interactions and fluctuating Brownian forces, as shown in Eq. (2). The most common starting point for implicit solvent models is a continuum treatment via the Navier-Stokes equations. Due to the small size of colloidal particles, inertial effects are insignificant (low Reynolds number), which leads to hydrodynamics well-described by the Stokes equations: where v is the fluid velocity, p is the pressure and μ is the fluid dynamic viscosity. We note that in writing Eq. (3), we have also assumed a Newtonian solvent. A Green's function (equivalent to the Oseen tensor) can be obtained for the Stokes equations, and a boundary integral equation representation can be used to reduce the three-dimensional system of Eq. (3) to a two-dimensional integral equation system [116]. This is the basis of the Boundary Element Method (BEM), which can be used to solve the flow field for arbitrary particle geometries, and as such represents a highly general approach to the treatment of Stokes hydrodynamics [116]. Unfortunately, the BEM is computationally too expensive for the time scales of interest in the present work, so we do not pursue it here. Several approaches have been devised to approximate the hydrodynamic interactions for systems that exhibit certain symmetries in particle shapes and relatively simple computational domains. A key simplifying feature of the Stokes equations is their linearity, which makes the generalized hydrodynamic force/torque vector F H a linear function of the particle translational/angular velocities U: Here, R is known as the hydrodynamic resistance tensor, which in the current formulation is a 6N -dimensional matrix that describes the effects of the velocity of any particle on the hydrodynamic force experienced by any other particle in the system. Additionally, we have implicitly assumed that no bulk flow fields are imposed. The Langevin equation for the colloid particles then becomes: Both hydrodynamic and Brownian forces originate from interactions between colloidal and solvent particles, so it is not surprising that they are closely related. This relationship can be derived based on the well-known fluctuationdissipation theorem, and yields the following requirement for the properties of the 6-N -dimensional Brownian force vector F B : Here, angular brackets denote an ensemble average, t and τ represent time, k B and T are Boltzmann's constant and the system temperature, respectively, and δ D is the Dirac delta function. Clearly, the many-body resistance tensor R is central for both hydrodynamic and Brownian forces. Much of the research effort in developing implicit solvent methods has been devoted to finding accurate and computationally efficient approximation schemes to this resistance tensor. The work of Ermak and McCammon [65] set the groundwork for many subsequent implicit solvent methods. The basic premise of their algorithm is to integrate Eq. (5) twice in time and determine the particle displacements over a time step Δt much larger than the inertial time scale of the colloid particles τ C = m/6πμa, where μ is the solvent viscosity, m is the particle mass and a is the particle radius. The interested reader is referred to the original work for details [65]. We note that this approach effectively ignores the inertial time scale of the colloid particles, which is not always a good assumption. In the current notation, the resulting expression for the evolution of particle positions r is given by: The inverse of the resistance tensor R −1 is known as the mobility tensor, M. The constraint from fluctuationdissipation theorem on the magnitude of the Brownian force F B in this formulation is: The simple forward time-stepping scheme in Eq. (7) can be problematic due to difficulties associated with computing the divergence of the mobility tensor (∇ · R −1 ). Several tractable alternatives exist, including the use of a midpoint time-stepping scheme [71,85]. The entries in the mobility tensor (or equivalently, the resistance tensor) can be derived from various analytical treatments of the Stokes equations (3). In the simplest case inter-particle hydrodynamic interactions are ignored, leaving only isotropic drag terms corresponding to the dilute limit. This yields resistance tensor entries R i j = 6πμRδ i j for all force-translational velocity couplings, and R i j = 8πμR 3 δ i j for all torqueangular velocity pairs. The accuracy of this simple treatment, often termed Brownian dynamics (BD), is expected to deteriorate when colloid particles are in close proximity. However, its convenient diagonal form is computationally expedient owing to the simple inversion of the resistance tensor. More sophisticated mobility tensors have been derived based on a multipole expansion treatment of the Stokes equations. In the original formulation of Ermak and McCammon, the Oseen tensor [90,235] as well as the Rotne-Prager tensor [186] were implemented. The resulting algorithms are also classified as Brownian dynamics, but modified based on various approximations for the hydrodynamic tensor. Such methods have met with some success and can often reproduce qualitative trends in suspension diffusion and rheological properties, but perform poorly for dense suspensions or when particles are in close contact. In near-contact regimes, lubrication forces can only be accounted for with high-order terms in the moment expansion, which is impractical. As it stands, the computational cost of these methods is fairly significant, due to the calculation of the resistance tensor, which scales as N 2 , and the inversion of the mobility matrix required for the Brownian displacement terms, which scales as N 3 [7,57].
The more rigorous Stokesian dynamics (SD) technique [25][26][27]59,75] has addressed many of the shortcomings in the accuracy of these techniques, albeit at even higher computational cost. The discussion so far has centered on a quiescent suspension of monodisperse spherical particles in an infinite medium. In its original form, the SD method efficiently accounts for near-field lubrication terms, considers the possibility of an imposed linear shear flow, and improves the far-field hydrodynamics with the inclusion of stresslet (S H ) to rate-of-strain (E ∞ ) coupling. For more details, the interested reader is referred to the original literature [27,59,75]. In Sect. 5, we discuss the application of SD to non-spherical particles and general geometries.
The Stokesian dynamics technique has proved extremely powerful, and remains one of the most accurate methods available for simulating multi-body hydrodynamics in the context of colloidal suspensions. However, due to its high computational cost it is currently extremely difficult to reach the time scales of interest for long-time diffusion. While the recently-developed accelerated Stokesian dynamics (ASD) techniques [8,200] have extended the timescales accessible to SD, the algorithm still scales poorly on parallel architectures. Variants and further improvements to ASD are available [48,72,107,125,153,200,225]. Additionally, several simpler pair-drag models [7] that are based primarily on pairwise lubrication terms have emerged as potential alternatives, particularly for high colloid volume fractions. The most recent variant of ASD, now commonly known as fast lubrication dynamics (FLD) [31,124,126], represents a significant simplification of SD and has shown considerable potential for use in predictive simulations of realistic systems. This technique is very closely related to an approximate version of ASD introduced by Banchio and Brady [8], which they dubbed ASDB-near-field, or ASDB-nf. We adopt the FLD nomenclature simply because more details are available on its implementation and performance [31,124,126] and it has been tested more extensively for the types of systems and the time scales that we are interested in here [126,197].
The FLD technique is based on splitting the resistance tensor R into an isotropic part that accounts for far-field multi-body effects, R 0 and a part that accounts for shortrange pairwise lubrication effects, R δ : The key simplifying assumption in FLD is choosing R 0 to be a diagonal tensor: The scalars f 0 FU , f 0 T and f 0 SE are functions only of the volume fraction of colloid particles φ, and are empirically fitted to match the short-time diffusivity and viscosity obtained from SD simulations. The details of the fitting procedure as well as other aspects of the method are described in greater detail in works by Kumar, Bybee and Higdon [31,124,126].
The lubrication component of the resistance tensor R δ is computed using only pairwise frame-invariant interactions [7] based on lubrication theory solutions of two-particle interactions [106,117], similar to what is done for the nearfield component of the SD algorithm. This gives rise to interactions that include terms of order δ as well as δ log(δ), where δ is the inverse of the distance between surfaces at the point of nearest approach for a pair of particles (δ = 1/(r i j −a i −a j ); a i and a j are the particle radii, and r i j is the distance between their centres). Following previous work [31,124,197], we have tested the algorithm both with and without the log(δ) terms.
The FLD method thus combines the pairwise near-field lubrication interactions that are dominant at high volume fractions and small inter-particle separations with a reasonable approximation for the far-field resistance tensor, which dominates at large particle separations. It inherits the flexibility of SD to account for imposed linear shear flows, but also its limitations with regard to more complex boundaries, flow conditions and solvent rheological properties. Although the discussion thus far has focused on monodisperse spherical particles, the extension to polydispersity is relatively straightforward, with pairwise lubrication expressions given by Kim and Karrila [116]. Extensions of SD and FLD to nonspherical particles will be discussed in Sect. 5, and more details pertaining to the implementation of FLD in this work are deferred to Sect. 3.

Conservative inter-particle interactions
The conservative force F c in Eqs. (1) and (2) is only a function of particle positions, and includes inter-particle forces as well as any effects from external forces. The physical basis of conservative inter-particle interactions is typically a very short-range repulsion due to steric effects (i.e. volume exclusion/particle collision), a short-range attractive force that has its physical origins in dispersion interactions, and a long-range screened electrostatic force, which may be either attractive or repulsive. For multicomponent solvents, attractive depletion forces due to the presence of additional solutes can also be significant. Typically, surface chemistries of colloidal particles and solvent properties (pH, ionic strength) are modulated in order to balance the attractive van der Waals force with a repulsive electrostatic force and promote particle dispersion. Conservative inter-particle forces can often be approximated as a sum of pairwise two-body interactions, each of which depend only on the separation r i j between particles i and j (as well as relative orientation for non-spherical particles). The total conservative force can thus be written as a sum of pairwise inter-particle forces and any external forces: The assumption of pairwise interactions ignores the possibility of multi-body effects, which can be problematic in some cases. For instance, colloidal particles coated with relatively large surface ligands (e.g. polymers) can interact in ways that alter the coating structure, which in turn affects their interactions with additional colloids. The discrete element framework can in principle be extended to include multi-body effects if their functional form can be determined, but the added complexity and expense are typically not justified.
The form of pairwise inter-particle forces F c i j (r i j ) is highly dependent on the details of the colloid particle material, colloid surface topology and chemistry, solvent properties and thermodynamic state variables. Typically, the potential energy U i j (r i j ) associated with the interaction of two particles is described, and the force is then trivially obtained from: Determining the functional form of U i j (r i j ) for two particles in a solvent requires an averaging (or coarse-graining) of all degrees of freedom other than the inter-particle separation. This may include molecular details of the solvent (e.g. counterion distribution, hydration layers), the colloid surface (e.g. passivating ligands/polymers, any condensed counterions) and the colloid interior particle structure. In a rigorous statistical mechanical framework, U i j (r i j ) is the potential of mean force between colloidal particles. Indeed, fully atomistic simulations of two colloidal particles can be used to compute the potential of mean force for a given system, which can then be tabulated and used in a larger scale simulation of a many-particle colloidal suspension [132]. This approach is computationally expensive, often suffers from statistical convergence issues, is only practical for systems with monodisperse, homogeneous particles and is only as accurate as the underlying atomistic force field. However, with improved force fields and increasingly powerful computing resources, it has the potential to become a key component in a multiscale modeling framework for colloidal suspensions.
The most common approaches to evaluating U i j (r i j ) rely on approximate analytical theories that describe the dominant physical processes responsible for inter-particle interactions. The classic DLVO theory treats screened electrostatic/electrical double layer interactions in the context of a linearized mean-field approximation for the counterion distribution in the solvent, and dispersion interactions in terms of summations of empirical forms of atomic dispersion potentials. A detailed derivation is not given here for either case, since these are available in many other works [49,63,68,104,224]. Several forms of the electrostatic interaction can be found in literature, and all contain a term that decays exponentially with inter-particle separation. An excellent discussion of these electrostatic potentials is given by Elimelech et al. [63]. In previous work [197], we have used the following form for the electrical double layer interaction potential for two particles with radius a separated by a center-to-center distance r : ρ ∞ is the bulk electrolyte concentration (assumed here to be a 1:1 salt, but other cases can easily be incorporated), κ is the inverse Debye screening length and ψ 0 is the electric potential at the particle surface, which can be related to experimentally measured zeta potentials. The interested reader is referred to the work of Schunk et al. [197], which includes a more detailed description of this potential and a discussion of how parameters can be selected to match experimental conditions. The colloid dispersion interaction components of the DLVO model are treated by summing the pairwise dispersion interactions between the constituent atoms/molecules in the particles. Once again, a variety of physical assumptions and approximation schemes have been used to carry out this summation [63]. Most formulations of DLVO theory treat only the attractive dispersion terms, which scale as r −6 . In the original formulation of Hamaker [89], the resulting attractive potential for two spheres of radius a separated by a center-to-center distance r yields the following for the attractive potential: Here and in previous work [197], we use the form derived by Everaers and Ejtehadi [68], which treats colloidal particles as collections of Lennard-Jones (LJ) particles, i.e. particles that interact according to: The parameters and σ are the usual LJ parameters representing the depth of the potential well and the characteristic dimensions of the LJ particles, respectively. The resulting integrated colloid-colloid potential U cc , which includes the r −12 repulsive component of the LJ potential, has the following form: The overall potential that we advocate for DLVO-type interactions is then the sum of the electrostatic component, Eq. (13), and the attractive and repulsive dispersion interactions, Eq. (16). These are implemented in the LAMMPS software package [176] as 'pair_style yukawa/colloid' and 'pair_style colloid', respectively.
While the qualitative features of DLVO theory have provided invaluable insights into the behavior of colloidal suspensions, its approximations can often be problematic. Real systems typically consist of particles that are heterogeneous in their composition, structure and surface characteristics. As such, even simple parameters like the effective particle radius a are not straightforward to assign for the purposes of calibrating the potential. Based on previous work [197], we have found the Hamaker constant A cc to be particularly difficult to estimate. Typical values based on the derivation discussed above [180] can vary by an order of magnitude depending on some of the underlying assumptions. Typically, several key parameters can be assigned based on experimental measurements of directly related properties, such as the zeta-potential and particle diameter; however, the Hamaker constant is often treated as a fitting parameter, which requires some level of calibration based on experimentally measured macroscopic properties. For the idealized systems in the present work, colloid particles are assumed to be made up of small Lennard-Jones particles with parameters σ and . The Hamaker constant can then be calculated analytically as A cc = 4π 2 ρ 1 ρ 2 σ 6 , where ρ 1 and ρ 2 are the number denisities of constituent LJ particles in colloid particles 1 and 2.
Most analytical treatments of dispersion potentials, including Eq. (16) lead to expressions that diverge at particle contact or near-contact. For dense suspensions under high shear or compressive forces, DLVO-like potentials must be replaced by short-range interactions for small particle separations. This is typically done using potentials derived in the context of granular flow simulations [33,105,202]. Various forms of such granular potentials have been proposed depending on various assumptions relating to the particle contact mechanics [51]. For most colloidal simulations in a liquid solvent, friction effects can be ignored, and simple frictionless Hookean or Hertzian granular potentials [45] are adequate. More complex potentials that also include friction effects [201] and inter-particle adhesion surface forces [110,145] can also be included.
Attempts have also been made at deriving closed-form potentials that include effects of grafted polymer ligands. For instance, Vincent et al. [226] derived an interaction potential between particles with grafted polymer chains in a solvent containing additional polymer molecules, which includes parameters such as the polymer χ -parameter, polymer adsorbed densities, molar volumes, etc. Despite the added complexity, such approaches can be useful in simulations for predicting general trends as a function of changes in various system parameters. However, the additional parameters can be difficult to evaluate, making quantitative predictions for such systems even more challenging.

Methods
As previously stated, the goal of this work is to directly compare several techniques for the simulation of colloidal suspensions. To make the comparison as straightforward as possible, we have selected the canonical system of hard-sphere monodisperse particles immersed in a Newtonian solvent.
In this section, we provide methodological details specific to the simulations carried out in this work. We therefore discuss our implementation of the MPCD and DPD methods, as well additional details pertaining to the FLD method described in Sect. 2.2.1. All of the simulations have been carried out in the LAMMPS software package [176], which is freely distributed as open source code (http://lammps.sandia.gov). We make reference to various modules within LAMMPS (known as 'styles') for the reader interested in applying some of these methods. The basic molecular dynamics algorithm and its parallel implementation in LAMMPS are not discussed here; it is covered extensively in other works [2,4,176]. However, the last portion of this section provides more details of algorithmic considerations specific to colloidal suspensions.
An important distinction in this work compared to previous similar studies is that we attempt to match quantitatively the physical properties of a particular solvent, rather than the usual approach of using convenient values for various parameters in the solvent models, and reporting all results in terms of dimensionless quantities [11,22,36,169,181,217]. For all simulations presented here, we select the solvent properties to match those of a Lennard-Jones fluid of density 0.66 σ −3 at temperature kT = and pressure P = 0 σ 3 / with a cutoff distance of r c = 3.0σ and particle mass m L J . In previous work, it was shown that this results in a dynamic viscosity of μ = 1.01 ± 0.03σ 2 /τ [172], where τ is the Lennard-Jones time unit (τ = m L J σ 2 / ). Note that we are not simulating Lennard-Jones particles at any point in this work, only using these solvent properties and Lennard-Jones units for convenience. The well-characterized Lennard-Jones system shares many of the same fundamental characteristics as real fluids, and the mapping procedure is in principle equivalent. This choice of solvent in combination with the hard-sphere-like interaction of colloid particles allows for a straightforward comparison of the different hydrodynamic treatments. The interested reader is referred to previous work [197] for discussions of the additional complications that arise from mapping to the properties of a more realistic system (polystyrene particles in water).

FLD implementation
The physical basis and mathematical treatment of hydrodynamic interactions and Brownian forces underpinning implicit solvent methods such as SD and FLD were discussed in Sect. 2.2.1. When implementing these methods in a molecular dynamics framework, the time evolution of colloid particles can be carried out without inertial terms, as was largely assumed throughout the development in Sect. 2.2.1. Alternatively, the inertial terms in Eq. (5) can be retained, and the equations of motion can be explicitly integrated in the usual manner. In previous work [197], we referred to these two schemes as implicit and explicit integration, respectively; here, we refer to them as 'non-inertial' (e.g. non-intertial FLD, or nFLD for short) and 'inertial' (iFLD), respectively. This is to avoid confusion with the terms implicit and explicit as they are applied to the solvent. As the name suggests, inertial schemes can resolve timescales shorter than the inertial timescale of the colloid, τ C = m/6πμa, but naturally require a time step significantly smaller than this. On the other hand, non-inertial methods require a time step significantly larger than τ C . The time step size is otherwise limited only by system-specific numerical considerations (e.g. conservative forces). We also note that both schemes assume that the underlying momentum relaxation in the solvent is instantaneous (i.e., the underlying hydrodynamic expressions are steady flow; see Eq. (3)).
We also present results for Brownian dynamics simulations (BD), which are implemented as a straightforward simplification of FLD. In particular, if all pariwise lubrication terms R δ in Eq. (9) are discarded, and volume fraction corrections to the isotropic terms are not carried out (i.e. f 0 FU and f 0 T are set to unity in Eq. (10)), the resulting simulation corresponds to Brownian dynamics (BD). The inertial version of BD (iBD) is then equivalent to Langevin dynamics (LD). All of these simplifications can be carried out in LAMMPS using various settings to the 'pair lubricate' styles.
The implementation of the FLD and BD schemes in LAMMPS follows the works of Kumar, Bybee and Higdon [31,124,126]. The discussion of Sect. 2.2.1 largely assumed quiescent flow conditions for simplicity, to which we now add the possibility of an imposed shear flow. In this case, the generalized velocities of colloid particles U are expressed relative to the velocity/angular velocity of the bulk fluid evaluated at the center of the particles, U ∞ . Also, the stresslet (S H )/rate-of-strain(E ∞ ) coupling introduced in the original Stokesian dynamics formulation is retained. The rate of strain tensor is a simple function of the shear rateγ and the flow geometry [126]. The hydrodynamic force and stresslet are given by the following linear relationship: We discussed the origins of the resistance tensor R in Sect. 2.2.1. The additional Brownian force is calculated based on Eq. (8), where the resistance tensor does not include the R F E component. This requires taking the square root of the force-velocity resistance tensor (R FU ), which is trivial for the isotropic, diagonal portion R 0 (see Eq. (10)), and can be carried out in a pairwise fashion for the lubrication component R δ [7]. The interested reader is referred to the relevant literature for additional details [7,31,124,126]. In the inertial version of FLD (iFLD), the Brownian force, hydrodynamic force and conservative inter-particle forces are summed, and the particle positions are updated using standard molecular dynamics schemes. In LAMMPS this is accomplished using "hybrid" pairwise interactions that include pair styles 'lubricate', 'brownian' and additional conservative pair styles (e.g. 'colloid' or 'yukawa/colloid'), in conjunction with a standard extended-body integrator (for translation and rotation) such as 'fix nve/sphere'.
For nFLD, the forces sum to zero, which leads to the following matrix problem for the particle velocities: Since the resistance tensor R FU is symmetric and positive definite, a conjugate-gradient algorithm can be used to solve for the particle velocities. The particle positions are then updated accordingly. In LAMMPS, this is accomplished using a similar hybrid pair style scheme, but with the use of pair style 'lubricateU' instead of 'lubricate', in conjunction with 'fix nve/noforce' to update particle positions.
The rate-of-strain is imposed for both iFLD and nFLD and the resulting stress which includes the solvent component can be measured. Once the particle velocities are known, the stresslet S H and the total stress in the solution τ can be computed readily in a post-processing step for a given configuration of particles (see Eqs. (13) and (14) in the work of Kumar and Higdon [126]). The viscosity μ r of the suspension relative to the viscosity of the solvent μ 0 can be computed as a function of the appropriate measured stress component τ xy and the imposed shear rateγ :

MPCD implementation
The MPCD method and its implementation in LAMMPS have been detailed elsewhere for a pure fluid [172] as well as for forced flow in the presence of walls [24]. A summary is given here for purposes of comparison, with additional details relevant to colloidal suspensions given in Sect. 3.5. In MPCD, the solvent is represented as point particles with no pairwise particle/particle or long-range interactions. In order to propagate momentum through the fluid, MPCD particles are first streamed at a constant velocity, and their positions x i are updated accordingly: At prescribed time intervals, particles are grouped into evenly-spaced cubic regions (bins), and various schemes are used to swap momentum among particles in order to simulate solvent collisions. In the original MPCD (or SRD) scheme proposed by Malevanets and Kapral [144], the components of the velocities of particles relative to the centre of mass velocity of all particles in a given bin ξ are rotated around a randomly selected orthogonal direction: Here, R s is a stochastic rotation matrix and u ξ is the center of mass velocity of all N ξ particles located in bin ξ : In our implementation, we select randomly one of six directions corresponding to positive and negative orthogonal axes, and always rotate by 90 • . Alternative rotation schemes have been proposed, and the rotation angle can be modulated to control the properties of the solvent [3,219]. It can be easily verified that this collision scheme conserves linear momentum and energy; as a result, the correct hydrodynamics are reproduced [101,144]. We refer to simulations based on this collision scheme as SRD.
An alternative collisional scheme known as multi-particle collision/Andersen thermostat (MPC-AT) [164,165] that does not entail rotation has also gained traction, due its inherent ability to thermostat the fluid. In this scheme, particle velocities are updated according to: Here, v i,rand represent velocities drawn randomly out of a Maxwell-Boltzmann distribution, and u ξ is the mean velocity for a particular bin, defined as in Eq. (22). A similar scheme can be used to conserve angular momentum (MPC-AT+a) [92,164], but we do not include it here as we are interested primarily in translational aspects of colloidal motion. In both SRD [92,219] and MPC-AT [164] collision schemes, the viscosity and diffusion coefficients of the pure fluid can be calculated analytically. The kinematic viscosities ν for the two schemes are given by: Here, M is the mean number of MPCD particles per bin, Δt is the time between collisions, Δx is the bin size and ρ is the mass density of the fluid. Values of ν and ρ are chosen to match the desired physical properties of the fluid, while  (24) and (25) yield two distinct solutions for the collision time step Δt [24]. The small time step solution corresponds to a short mean free path λ relative to the bin size, which requires random shifting of collision bins to remove inter-particle correlations and restore Gallilean invariance [100,101]. In contrast, the larger time step solution leads to a larger mean free path, and bin shifting is not required. However, the short mean free path solution also corresponds to higher, more realistic Schmidt Here we present results for several sets of parameters for both methods. In all cases we attempt to match the properties of the same Lennard-Jones solvent discussed earlier, following previous work [24]. The resulting parameters for both collision schemes are summarized in Table 2. The Schmidt number in Lennard-Jones fluids is typically of order 10-100 [152], comparable to values obtained using the small collision time step, but significantly larger than that of the large time steps. More realistic fluids have Schmidt numbers of order 10 3 , which is difficult to attain with the MPCD method. Several schemes have been proposed for coupling MPCD fluids to colloid particles. The simplest of these includes colloid particles in the multi-particle collision scheme as if they were another MPCD particle (located at the colloid center) [11,109], with the option of weighting the colloid particle influence by the colloid mass. This is referred to as 'collisional coupling' by Batôt et al. [11]. While this scheme yields qualitatively reasonable results, it does not resolve the particle surface in any way, and cannot be expected to compare quantitatively to the other methods tested here; the hydrodynamic radius of the colloid particles in this approach is a function of the MPCD parameters, rather than being set to the desired value [109]. A more sophisticated approach relies on introducing a pairwise potential between MPCD particles and colloid particles, typically with a short repulsive character [11]. However, this introduces additional parameters that must be carefully selected for a given set of MPCD fluid parameters [11], and it is not clear that no-slip boundaries can easily be enforced in this manner. Finally, the approach that we adopt here involves detecting collisions between MPCD particles and colloid particles that take place during the MPCD streaming step (Eq. 20), re-directing the MPCD particles, and adjusting the forces/torques on the colloids appropriately. These collisions are not to be confused with the bin-wise multi-particle collisions, in which only MPCD particles participate (Eqs. 21 and 23). The physical basis of MPCD particle-colloid collisions is not easily interpreted, since MPCD particles do not correspond directly to solvent molecules or fluid elements. One of the key requirements of such collisions is that they enforce no-slip boundary conditions at the colloid particle surface, i.e. the average tangential component of the MPCD fluid velocity must vanish at the colloid surface. We have tested two such collision schemes.
The first collision scheme, which we refer to as stochastic boundary conditions, is based on the work of Inoue et al. [103]. Briefly, an MPCD particle that collides with a colloidal particle is returned to the point of collision, after which it is assigned a random velocity with a component along the outward normal to the colloid particle surface, and two components tangential to this. The normal and tangential components (v n and v t , respectively) are drawn out of the following distributions: Here, β = m f /(2k B T ), with m f the mass of the MPCD fluid particles (taken here to be unity). The colloid particle surface velocity at the point of collision (Eq. 29) is then added to the randomly assigned MPCD particle velocity, and the MPCD particle is propagated with the new velocity for the remainder of the streaming time step. Forces are applied to the colloid at the collision point based on the momentum change of the MPCD particle. In previous work [24], we showed that this boundary condition can be problematic for situations with forced flow, but we expect it to be adequate for quiescent conditions. In the second collision scheme ('reverse'), the MPCD particle is returned to the point of collision, and the magnitude of each of its outgoing velocity components is calculated assuming a momentum-conserving (perfectly elastic) collision with the point on the colloid surface. In the limit of the colloid mass being much larger than the MPCD particle mass, this yields: Here, v f 0 and v f 1 are the velocities of the MPCD particle before and after the collision, respectively, and v 0 cs is the velocity at the colloid surface before the collision: v 0 c and c are the translational and angular velocities of the colloid particle prior to collision, and r s and r c are the locations of the surface collision point and the colloid particle centre, respectively.
In order to avoid artificially low local viscosities in MPCD bins that are partially occupied by colloid particles [131,232], we use the "virtual particle" method first suggested by Lamura et al. [131]. Although several versions have been proposed [24,232], we only implement one here, which has been shown to be adequate for forced flow between parallel walls [24]. In this scheme, "virtual" MPCD particles (VPs) are added at random locations to the interior of colloid particles where MPCD bins partially overlap colloids. Next, the velocities of the VPs are randomly assigned out of a Maxwell-Boltzmann distribution, and augmented by the velocity of the interior point in the colloid corresponding to their location. VPs then participate in multi-particle collisions with actual MPCD particles, and resulting changes in their momenta are transferred as forces/torques to the colloids. Note that VPs only participate in the multi-particle collisions, and new VPs are generated for each such event. The total number of VPs assigned to each colloid particle is selected to yield the same overall density of VPs as actual MPCD particles in the bulk fluid. Additional details are provided in our previous work [24], where this method was referred to as "VP dens " and "VP multi " .

DPD implementation
The DPD method introduced by Hoogerbrugge and Koelman [95] treats a fluid using a collection of particles with 'soft' interactions, allowing for a much larger time step than atomistic MD. The characteristic dimension of DPD particles (i.e. the cutoff of the DPD interparticle interactions) is much larger than the molecular dimension of fluid molecules, but should be smaller than the relevant flow characteristic size (in the case of suspensions, the colloid particle diameter). Additionally, thermostatting in a DPD fluid is carried out using a pairwise scheme that preserves hydrodynamics. The physical basis of the DPD method can be justified if DPD particles are thought of as 'packets' of fluid, but as with MPCD, rigorous mapping to the properties of a given fluid and resolving geometric features comparable to the size of DPD particles is problematic [175].
From a molecular dynamics perspective, DPD is a traditional pairwise interaction model with a short-range cutoff and can be implemented similar to a Lennard-Jones potential. The implementation in LAMMPS follows that of Groot and Warren [87], where the force on particle i due to particle j separated by a distance r < r c is given as: is the relative velocity of the two particles, α is a Gaussian random number with zero mean and unit variance, Δt is the timestep size, and w(r ) is a weighting factor that varies between 0 and 1.
Like all explicit solvent models, mapping the physical properties of a real fluid to a DPD model is challenging. In particular, the lack of analytical expressions for various dynamic properties such as viscosity complicates matters (i.e. equations analogous to (24) and (25) for the MPCD method). Approximate expressions have been proposed [87], which provide a useful starting point for the mapping procedure we carry out herein. As previously discussed, we aim to reproduce the viscosity of the carrier fluid, rather than simply the key dimensionless parameters for a given flow situation. The combination of DPD parameters that gives the desired viscosity of 1.01m L J /σ τ is not unique. We therefore select a value of A = 25 /σ , which, following previous work [170,181,197], results in a realistic solvent compressibility; the DPD interaction cutoff r c is chosen to be somewhat smaller than the colloid particle diameter [181], r c = 3.0σ ; and the DPD particle number density is set to 3.0r −3 c , or 1/9σ −3 based on computational considerations. The DPD particle mass is set to 5.94m L J to match the target mass density of 0.66m L J /σ 3 . Finally, the parameter γ is set to 86 τ/σ 2 , which yields a viscosity of ∼ 1.01m L J /σ τ , as measured using a trial-and-error approach and the Müller-Plathe viscosity calculation method (see below). The nonlinear dependence of the viscosity on the damping parameter γ is well-known for the DPD method [114,218]. For this set of parameters, the system pressure is ∼ 3 /σ 3 , and the Schmidt number is ∼ 64.
Several approaches have been proposed for coupling a DPD fluid to colloidal particles. In early works based on the DPD method, suspended objects were simulated by constraining clusters of DPD particles to move as rigid bodies [21][22][23]118]. While this approach is flexible with regard to representing various object shapes (see also Sect. 5.2), it leads to significant solvent penetration into colloid par-  ticles due to the soft DPD interaction, and in the case of spherical colloids, it entails unnecessary computational expense. Instead, we model colloids as single spherical particles, and introduce an additional, relatively stiff conservative interaction between colloids and DPD particles. Similar approaches have been proposed by others and shown to be more conducive to producing correct hydrodynamic interactions [36,60,170]. Given the non-physical nature of DPD particles, the choice of interaction potential is somewhat arbitrary, and serves only to prevent solvent penetration into the colloid and satisfy the required boundary conditions. As such, we use a simple shifted Lennard-Jones potential: Here, C S and σ C S are Lennard-Jones parameters specific to the colloid-DPD interaction, r is the distance between the colloid and DPD particle centers, and δ C S is a parameter that effectively shifts the potential so that it diverges at a non-zero particle separation. As these parameters are not straightforward to select, we have tested several combinations to achieve the desired diffusion characteristics. These parameters were selected to approximately match the short-time diffusivity of colloidal suspensions at a volume fraction of 0.1 based on iFLD simulations; short DPD simulations were therefore carried out with several parameter sets to ensure that reasonable agreement was attained. We have tested three combinations of parameters for DPD-colloid coupling, which we denote as DPD-v1, DPD-v2 and DPD-v3. The values of the relevant parameters are summarized in Table 3.
Additional modifications of the DPD method have been proposed in literature that have resulted in improved accuracy. The centro-symmetric colloid-solvent interaction potential in Eq. (31) can only transfer linear momentum to colloid particles, and all rotational aspects of the flow are lost. To address this, Espanol [67] has proposed a generalization of the DPD method known as the fluid particle model (FPM), which adds non-central shear components to the dissipative forces. Several variations of this method have been successfully applied to simulations of colloidal suspensions [170,175,181]. In order to achieve higher Schmidt numbers, Fan et al. [69] have suggested a modified form of the DPD conservative interaction, w(r ) = (1 − r/r c ) s , where non-unity values of the exponent s have been shown to yield higher Schmidt numbers. For high-viscosity fluids, the use of a Lowe thermostat [139] to replace the DPD random force has been advocated by several workers [35,36]. For additional discussion, the interested reader is referred to the excellent review by Pivkin et al. [175].

Colloid particles and simulation details
In all of the methods discussed here, the colloidal particles are treated as finite-sized particles with translational as well as orientational degrees of freedom (LAMMPS atom style 'sphere'). Parameters for the hard-sphere-like colloids are similar to those in previous work [197], with the particle radius set to a = 5σ . Inter-colloid interactions are captured using the integrated Lennard-Jones potential (Eq. (16); pair style 'colloid' in LAMMPS) with a cutoff distance of R c = 2a+30 −1/6 σ and a Hamaker constant of A cc = 4π 2 ∼ 39.478 . Note that this interaction is not infinitely hard, but is slightly softened to prevent colloid overlaps and allow for standard molecular dynamics time integration. Physically, this corresponds to adding a short surfactant coating to the colloid particles to avoid flocculation. Due to its short range and relatively stiff nature, the differences from a true hardsphere interaction are expected to be minor. Although lubrication forces prevent the close approach of the colloids, we also include the inter-particle potential with FLD, in order to maintain consistency with the other methods. In all cases, the inner and outer cutoff distances for FLD lubrication terms are set to 2.0002a and 3a, respectively (i.e. FLD lubrication terms are active for 2.0002a < r i j < 3a).
Based on simple energy conservation considerations, we found that the hard-sphere potential described above limits the time step to a value of Δt C ∼ 0.04τ . This is smaller than the inertial timescale of colloid motion (τ C = m/6πμa ∼ 3.63τ ), implying that nFLD cannot take full advantage of its large time-stepping ability; indeed, results for these methods for timescales not satisfying t >> τ C are physically suspect. In order to facilitate the comparison among various methods, we nonetheless retain the colloid potential and a small time step of δt C = 0.01τ C ∼ 0.0363τ for iFLD and δt C = 0.005τ C ∼ 0.0182τ for nFLD simulations, where colloid particle overlaps were occasionally problematic. We will return to this discussion when comparing computational performance of the different methods, and note that for other conservative inter-particle potentials, nFLD can be much more efficient. Similar values of the MD time step (∼ 0.03τ ) were used for MPCD, where the multi-particle collision time step (see Table 2) must be an exact multiple of the colloid/MD time step. For DPD, where no such constraints exist, we used a time step of 0.03τ for all cases.
The starting configurations for our simulations were generated so as to maintain consistency between the different methods. The MPCD method is the only one that imposes any constraint on the simulation domain size, which should be a multiple of the MPCD bin size (2.0σ ). As such, we set the simulation box size for all cases based on this criterion. Computational considerations for the explicit solvent methods (MPCD and DPD) at low volume fractions limit the system size to several hundred colloid particles, so we use a system size of 200 particles of radius a = 5σ for all methods and volume fractions. Based on these two constraints, we used cubic domains of side length 104, 80, 70 and 64σ , corresponding to volume fractions φ of 0.0931, 0.2045, 0.3053, and 0.3995, respectively; for simplicity, we refer to these as 0.1, 0.2, 0.3, and 0.4 volume fraction systems. The diffusion coefficient measured in simulations with periodic boundary conditions is known to have a strong dependence on system size [8,75,129]. We therefore use the corrections suggested by Ladd [129] for all volume fractions. Further details are provided in the Sect. 4.
The lowest volume fraction system (0.1) was constructed by random placement of colloidal particles in a cubic domain of size 104σ , ensuring that no particles overlap. This method is robust up to volume fractions of approximately 0.3, beyond which random placement of particles quickly exhausts all possibilities. As such, the 0.1 volume fraction system was slowly compressed in a simple Langevin dynamics (LD) simulation to the desired dimensions to generate the starting colloid particle configurations for the remaining volume fractions. In the case of MPCD simulations, additional solvent particles were added to achieve the desired number of MPCD particles per bin (M = 5.28), while ensuring that no MPCD particle started inside a colloid particle. For DPD simulations, simulation boxes of the same size as the desired colloid suspension systems were first constructed using only DPD particles with the requisite parameters (see Sect. 3.3), and the resulting pressures were measured. The colloidal particles were then added, and any DPD particles within a distance a of any colloidal particle center were removed. For each volume fraction, the pressure of the resulting colloid-DPD mixture resulting from this procedure was found to be in good agreement with the pressure for the pure DPD fluid; as such, the viscosity and other key parameters of the DPD fluid are consistent among volume fractions.
In all simulations, several hundred thousand equilibration steps were carried out prior to sampling for the purposes of computing diffusion coefficients or radial distribution functions. Viscosity for all FLD methods was determined using a non-equilbirium molecular dynamics (NEMD) approach, where the simulation box was shear-deformed at a constant rate (equivalent to Lees-Edwards boundary conditions; implemented in LAMMPS using the "fix deform" capability). The resulting stress was then measured as outlined near the end of Sect. 3.1, and viscosity values were computed for several shear rates. The Müller-Plathe velocity swapping algorithm [160] was used to compute viscosities for all explicit solvent methods. In this method, a momentum flux is imposed by random exchanges of particle velocities, and the resulting velocity profile (i.e. shear rate) is measured. The shear rate can be varied by controlling the frequency and number of particles undergoing velocity exchanges. For additional details, the interested reader is referred to the original work of Müller-Plathe [160] as well as the work of Petersen et al. [172]. The method is implemented in LAMMPS in the "fix viscosity" capability. For both MPCD and DPD, only solvent particles are included in the velocity exchanges; this minimizes interference with the colloid particle dynamics, and leads to a much faster convergence of the measured viscosities as compared to the NEMD method that must be used for FLD simulations. In both DPD and MPCD, simulations of the pure fluid were first carried out at the desired shear rates to ensure that no shear thinning of the solvent takes place in the shear range of interest.

Computational considerations
From a computational perspective, all the techniques described in this section (FLD, MPCD, DPD) share similar computational kernels: computing colloid/colloid interactions and (optionally) colloid/solvent or solvent/solvent interactions. This enables all the methods to be implemented in a discrete element (DEM) or molecular dynamics (MD) framework such as LAMMPS. Moreover, all the methods can exploit parallelism using the spatial-decomposition approach provided by LAMMPS (and many other MD codes) whereby the simulation box is partitioned across processors and processors compute forces on particles within their subdomain and communicate particle information to processors owning neighboring sub-domains [176]. Due to the shortrange nature of the computations, so long as there are sufficient particles per processor, the computational work outweighs the communication costs, and the parallel performance of a simulation can scale to large numbers of processors. Prototypical performance and scaling data are presented in Sect. 4.4.
Since FLD is an implicit solvent method, colloid/colloid interactions are the only pairwise particle computations to perform. As discussed above, the inertial FLD method is effectively an explicit time-stepping method, so the forces and torques on particles resulting from pairwise interactions are time-integrated in the usual manner. The non-inertial version of FLD is effectively an implicit time-stepping method, which requires solution of a linear matrix equation to obtain particle velocities at each time-step. This is a sparse matrix of order N = the number of particles, with non-zero elements for each particle/particle interaction within the cutoff distance. The matrix equation can thus be solved in a modest number of conjugate gradient iterations which invoke the same kind of local inter-processor communication needed to exchange particle information between neighboring processors in the three-dimensional spatial partitioning.
MPCD for hybrid colloidal/solvent systems has been implemented in the parallel MP2C code by Sutmann et al. [211], and has been used to model hydrodynamic effects on shearing polymers [97] and colloidal systems [205]. In LAMMPS, MPCD is implemented as a 'fix srd' style, which allows MPCD solvent particles to be added to a colloidal system. Colloid/solvent interactions are only computed when a streaming solvent particle collides with a colloid particle. As discussed above, the no-slip collision imparts force and torque to the colloid. The manner in which collisions are efficiently detected in LAMMPS is by binning the colloid particles onto a regular 3-D grid, which may but need not be the same resolution as the grid used for binning MPCD particles to perform multi-particle collisions. Since colloid particles have finite extent, each may overlap several such bins. After an MPCD particle streams to its final position, a loop over the small number of colloid particles in its bin is used to check for collisions.
MD codes such as LAMMPS use neighbor lists of nearby particles to efficiently enumerate particle/particle interactions [176], e.g. for colloid/colloid interactions in all the methods discussed here. MPCD particles do not interact directly with each other, and the number of solvent particles is orders of magnitude larger than the number of colloid particles in typical mixture models. For performance reasons it is important to be able to "skip" the solvent particles when building neighbors lists for colloid/colloid interactions or other operations such as time integration, which only involve colloid particles. LAMMPS does this by separating the two kinds of particles within the list of particles owned by each processor. Colloid particles are at the beginning of the list and can thus be looped over, when necessary, as if no solvent particles were present. Likewise, no MPCD particles need be communicated to neighboring processors to act as "ghost" particles when computing pairwise interactions.
As previously mentioned, DPD is straightforward to implement as a short-range pairwise interaction potential in an MD code. However, an additional complication arises when using small interacting particles as a background sol-vent for large colloidal particles. Such a model has multiple length scales, as evidenced by the disparate cutoffs for colloid/colloid, colloid/solvent, and solvent/solvent interactions. For example, if colloid particles are 20-times the size of solvent particles (diameter = σ ), then coarse-grained colloid/colloid potentials like that of Everaers and Ejtehadi [68] may be cut off at a distance ∼2.5× the colloid diameter = 50σ . Solvent/solvent interactions (like Lennard-Jones or DPD) are typically cut off at ∼2.5σ , and colloid/solvent interactions at some intermediate distance, ∼25σ . Standard parallel algorithms for building neighbor lists and communicating ghost atoms, designed for systems with a single cutoff length, can become very inefficient for such systems. To overcome this, LAMMPS has optional 'multi'-style neighborfinding and communication algorithms tailored for systems with multiple cutoff lengths. Benchmarks of prototypical solvated colloidal systems with size disparities up to 20× showed speed-ups of up to 100× in serial or parallel, using the algorithms presented by in't Veld et al. [222]. Although the disparities in particle sizes and cutoff distances are not as significant for the systems in the present work, these algorithms are nonetheless beneficial.
A final computational issue that can affect all of the methods, when running in parallel, is that of load-imbalance. Colloidal particles in dilute systems (low volume fraction) may aggregate, whether the solvent is implicit or explicit. This can lead to unequal numbers of particles per processor, if the simulation box is partitioned into equal-sized sub-domains. The problem can be exacerbated if explicit solvent is used, since both the large coarse-grained colloid and tiny solvent particles are represented as single particles. Thus the number of total particles per processor may vary widely. This typically becomes less of an issue as the colloid volume fraction rises, since the colloid density remains more spatially homogenous. LAMMPS has a 'balance' option to adjust processor sub-domain sizes to improve load balance, which can compensate for one-dimensional density variations, e.g. in an evaporation model with a liquid/vapor interface [38].
The use of GPUs for computational acceleration is becoming common place in molecular dynamics modeling. There is recent work by Wang et al. for GPU-enabled DPD models [228], and by Westphal et al. for GPU-enabled MPCD modeling [231] for pure MPCD fluids. We are not aware of GPU implementations of hybrid colloidal/MPCD models.

Diffusion coefficients
The equilibrium diffusion of colloid particles is fundamental to more complex phenomena of interest, such as shear viscosity and microstructure formation. If any simulation tech- nique is to be quantitatively predictive, it must first be able to predict diffusion properties accurately. As such, we present comparisons of the various methods discussed above, with a particular focus on diffusion. The details of each simulation method, the treatment of the colloidal interparticle interactions and the system construction have been discussed in the Sect.3.
In Fig. 1, we plot diffusion coefficients of hard-sphere colloid suspensions at various volume fractions as a function of time based on FLD simulations, using both nFLD and iFLD formulations discussed in Sect. 3.1. These plots were obtained by taking central finite difference numerical derivatives of the mean square displacement of colloid particles as a function of time: The angular brackets denote an average over particles as well as times τ 0 , i.e. using a standard moving time-origin analysis. Note that these values are not corrected for finite size effects, as these corrections are only applicable to plateau values of the diffusion coefficients, rather than the entire range of timescales.
Of particular interest are the early-and late-time plateau values of the diffusion coefficient, indicated on each plot by a cross or circle symbol, respectively. For the non-inertial FLD methods (nFLD), the early time regime is not physically real-istic, as this method is predicated on the assumption of having a much larger time step than the colloid inertial time scale τ C = m/6πμa; this point was discussed in greater detail in Sect. 3.1. Nevertheless, for sufficiently long times, the nFLD diffusion coefficients converge exactly to the iFLD values. The physical origins of these different diffusion regimes are well-known-the early-time plateau corresponds to the diffusion of colloids on length scales characteristic of interparticle separations, i.e. diffusion within a cage-like structure formed by neighboring colloid particles. The decrease in the diffusion coefficient at later times is associated with diffusing particles encountering neighboring particles and being unable to move past them. At sufficiently long times, particles diffuse past their neighbors, and a Fickian diffusion regime is recovered. Note that the very short-time ballistic regime encountered in atomistic systems is not accessible to either FLD variant, as this requires resolving a time scale shorter than colloid-solvent collisions, and the form of the Brownian force used in both variants of FLD assumes a much larger time scale. At the shortest times shown in Fig. 1, the diffusion coefficient is independent of volume fraction, as this corresponds to diffusion on length scales smaller than those at which colloid particles interact. With increasing volume fraction, inter-colloid collisions happen on shorter length (and time) scales, so that early-time values of diffusion are reached faster, and D early values are lower; similarly, late-time diffusivity values D late decrease, as higher colloid densities lead to greater hindrance to diffusion.
Qualitatively similar plots to those in Fig. 1 are obtained for all simulation methods. For ease of discussion, we present only D early and D late values as a function of volume fraction; while there are additional small differences in the behavior of the diffusion coefficient as a function of time, these plateau values suffice to provide quantitative comparisons of the different methods. Additionally, the values of the early-and latetime diffusion coefficients obtained in this manner are known to be sensitive to the number of colloid particles N in the system, with a dependence that scales as N −1/3 [8,75,129]. Laddy [129] has proposed an accurate correction for these finite-size effects, where the true infinite-size diffusion coefficient D ∞ is given by: Here, D(N ) is the diffusion coefficient measured based on the plateau values in Fig. 1, D 0 = k B T /6πμ 0 a is the Stokes-Einstein diffusion coefficient, μ 0 and μ are the viscosity of the pure solvent and the suspension, respectively and φ is the volume fraction. For the viscosity of the suspension μ, we use the shear viscosity measured at Pe = 0.1; see Sect. 4.2.
For simplicity, we do not retain the D ∞ notation, and leave it  to be understood that diffusion coefficients for FLD, MPCD and DPD in all subsequent discussions have been corrected in this manner. In Fig. 2, we present early-and late-time diffusivity values for the FLD and LD methods (see Sect. 3.1). The noninertial versions (nFLD and BD) yield the same results (data not shown). Results from theory [216,223], experiments [150,168] and other simulation works [74] are presented as reference values. We include Brownian dynamics results from the work of Foss and Brady [74], which, in agreement with our own LD simulations, show the need to account for hydrodynamic effects. Note that LD results need not be corrected for finite size effects. While some variation exists among theoretical, experimental and simulation reference values for both D early and D late , it is clear that the iFLD method yields diffusion coefficient values in excellent agreement with the expected values. In particular, the agreement in late-time diffusion coefficients between the iFLD method and the more rigorous Stokesian dynamics method is excellent. We therefore use the iFLD results as a benchmark for the remaining methods.
In Fig. 3, we plot the early-and late-time diffusion coefficients obtained using several variations on the MPCD method, all of which were discussed in Sect. 3.2. We have also tested the variants with the inclusion of "virtual particles" (see discussion near the end of Sect. 3.2), but found no appreciable effect on diffusion characteristics. Even though the key physical properties of the background MPCD fluid (viscosity, density) are the same in all cases, and the colloid particle boundaries are equally well-defined, there are clearly strong effects of all the methodological details tested.
While no clear systematic trends emerge from Fig. 3, several general comments can be made. First, with the exception of the SRD/large collision time step/reverse boundary conditions scheme, all other variations attempted here fail to reproduce the iFLD results. However, in most cases the discrepancy is not drastic (∼20-30 %). Second, stochastic boundary conditions (Eqs. 26 and 27) always lead to larger colloid diffusion coefficients as compared to reverse boundary conditions (Eq. 28). In all cases, this leads to an overestimate of the late-time diffusivity. We therefore advocate the use of reverse boundary conditions, and limit the discussion of viscosity that follows to simulations based on this treatment. Third, the small collision time step generally leads to larger values of the diffusion coefficient, and a better match to the iFLD results is obtained with the large collision time step. This is somewhat surprising, as the large collision time step leads to unrealistically low Schmidt numbers for both algorithms, but the effect appears to be minor. Fourth, the MPC-AT collision scheme tends to yield diffusion coefficient values that are mostly larger than those resulting from the SRD collision scheme.
Overall, the MPCD method yields qualitatively reasonable results, and in one case (SRD/large collision time step/reverse boundary conditions), a near-perfect match to iFLD data. However, given the differences observed as a result of the different variants, it may well be that the match in this case is fortuitous. It is not clear that the finite size corrections given in Eq. (33) apply equally to all MPCD variants, particularly at finite Schmidt number, where hydrodynamic correlations may decay much differently. Additional simulations with larger numbers of colloid particles may elucidate this, but these are beyond the scope of the present   Fig. 4, we summarize colloid diffusion results for DPD using three sets of colloid-solvent coupling parameters (see Table 3). We were unable to obtain late-time diffusion coefficient values for a volume fraction of 0.4 (the diffusivity drops to near zero at long times), and the short-time values at this volume fraction are likewise unreliable. This is a result of the finite size of DPD particles, which results in an artificial gel-like state of the system at high volume fractions. At lower volume fractions, the agreement with FLD values is quite good; both early-and late-time values are best reproduced by the parameter set DPD-v3, while DPD-v1 and DPD-v2 overestimate the diffusion coefficients. Presumably, with an alternate choice of parameters or an alternate functional form of the DPD-colloid interaction potential (viz. Eq. 31), even better agreement for both short-and late-time diffusivity can be attained. Additionally, some of the more sophisticated versions of DPD discussed in Sect. 3.3 (FPM method and Lowe-Andersen thermostat) can be expected to yield further improvements. However, in all cases, parameters and coupling schemes cannot be chosen directly, and must be calibrated using an ad hoc approach. We therefore identify the lack of a systematic procedure for the selection of DPD-colloid coupling parameters as a major shortcoming of the method.

Shear viscosity
Shear viscosity was computed for all simulation methods, using a NEMD/Lees-Edwards approach for FLD and the Müller-Plathe method [160] for DPD and MPCD. Additional details are provided in Sect. 3.4. In all cases, we express the shear rateγ in terms of the Peclet number, i.e. the ratio of the applied shear forces to diffusion forces (Pe = a 2γ /D 0 = 6πμ 0 a 3γ /k B T ). We note that this is only for numerical convenience, since all other quantities (μ 0 , a) are equal for all cases. In Fig. 5, we plot the suspension viscosity normalized by the solvent viscosity μ 0 as a function of Pe for various volume fractions, as calculated using iFLD and MPCD. In all cases, nFLD and iFLD yielded the same results to within statistical error, so we only show data for the latter.
At low volume fractions, the shear viscosity computed by MPCD is in excellent quantitative agreement with that produced by iFLD. As the volume fraction increases, the agreement remains strong, but with some notable discrepancies. The small collision time step for both the SRD and MPC-AT collision schemes results in a lower suspension viscosity in all cases as compared to the iFLD values; the agreement is better with the large time step. This is to be expected given that the small collision time step simulations tend to overestimate the diffusion coefficient (indicating a lower effective viscosity). In all cases, the well-known shear thinning effect is observed at higher volume fractions and moderately Early-time diffusivity Late-time diffusivity high shear rates [39,233]. This effect has traditionally been explained as a result of shear-induced layering of particles, although recent work suggests alternative mechanisms [233]. Overall, the MPCD method appears to be an excellent tool for the prediction of rheological properties of colloidal suspensions; however, it is limited to moderate shear rates such as those tested here (Pe < 100), as higher shear rates result in shear thinning of the MPCD fluid (data not shown); while many fluids of interest, including the LJ solvent tested here do exhibit shear thinning, it occurs at much higher Pe values. Figure 6 shows a similar comparison of shear viscosity for iFLD and DPD. For clarity, we only show results for DPD-v2 and DPD-v3, since DPD-v1 diffusion results are similar to DPD-v2, but consistently overestimate the diffusion coefficient. At low volume fractions, the agreement with iFLD is strong for both DPD-v2 and DPD-v3, as expected given the agreement in the diffusion coefficient values. However, at higher volume fractions, the agreement breaks down significantly, as was the case for diffusion characteristics. The higher volume fraction cases correspond to a regime in which the typical inter-colloid separation is on the order of the size of DPD particles, which leads to a highly structured gellike state. For sufficiently high shearing rates, the collective motion of colloid and DPD layered particles yields signifi- cant shear thinning. However, suspension viscosities at low Pe values cannot be reliably measured, as the statistical convergence of the velocity profiles requires unfeasible simulation time. At high volume fractions, we were not able to reach a sufficiently low Pe value to obtain the zero shear rate value of the viscosity. As such, DPD is an adequate method at low volume fractions, but breaks down at volume fractions above ∼ 0.3. The choice of DPD-colloid coupling method will not yield significant improvements in this area, as this effect is inherently related to the finite size of DPD particles. An obvious solution is to decrease the size of the DPD particles (r c ) relative to colloid particles, but this requires a higher density of DPD particles to maintain the same viscosity, and therefore a significant increase in computational costs.

Equilibrium colloid microstructure
As a final point of comparison, we also quantify the equilibrium microstructure of colloid suspensions using the radial distribution function g(r ). In simulation studies of colloid suspensions, it is often of interest to compare pair distribution functions in non-equilibrium conditions (i.e. under shear), particularly in different directions (e.g. velocity gradient and vorticity directions) [39,126,233]. We do not carry out such an analysis here, as it does not advance our methodscomparison appreciably. Equilibrium g(r ) plots are shown for the various methods at volume fractions of 0.1 and 0.4 in Fig. 7. The radial distribution function data corresponding to MPCD is only shown for one variation of the method (MPC-AT collision scheme, small collision time step, reverse boundary conditions), as other variations do not lead to appreciable differences in g(r ).
Clearly, the FLD and MPCD methods yield the same equilibrium colloid suspension structure at both low and high volume fractions; DPD on the other hand predicts a significantly different structure. With FLD, no solvent particles are present, and FLD pairwise and isotropic terms only modify the dynamics of the suspension, but do not affect structural properties. Similarly, the point-mass particles used to represent solvent in all MPCD methods do not have a significant effect on equilibrium structural properties; at sufficiently high MPCD number densities, depletion-like forces can be expected for particles in close contact, but they do not appear to be significant here. This is likely because the inter-colloid repulsive potential used here prevents such close particle approaches from occurring. With DPD, the finite size of the solvent particles and the conservative forces used to couple DPD and colloid particles clearly have a strong effect on the structural and thermodynamic properties of the suspension. The inter-colloid interaction potential is effectively modified by the presence of DPD particles, as shown by the resulting radial distribution functions. At high volume fractions, where the DPD particles and the range of DPD-colloid interactions approaches the inter-colloid separation, this leads to a highly-structured, gel-like suspension, which explains the breakdown of the DPD method for diffusion and shear viscosity calculations in this range. A related and more detailed discussion of the effects of explicit versus implicit solvents on inter-colloidal interactions is given by Grest et al. [86].
Overall, it is clear that the FLD method yields more accurate diffusion coefficients than any variant of the MPCD and DPD methods that we have tested. A significant challenge for the latter two methods is reproducing the desired physical characteristics of the fluid, in particular dynamic properties such as viscosity. Given the variations in diffusivities observed as a result of changes in the solvent-colloid coupling scheme for both MPCD and DPD, it appears that enforcing the correct colloid particle size and boundary conditions is also problematic. The precise nature of the difficulties in this regard is not clear from the results presented here for either MPCD or DPD, and requires additional investigation and development of these methods. With FLD, both the solvent viscosity and particle size are direct inputs to the simulation, which makes for a trivial parameterization procedure.
Despite discrepancies between explicit solvent methods and FLD with regard to diffusion coefficient values, viscosity data for certain variants of MPCD are in excellent agreement. It appears that for the system tested here, the unusually low Schmidt numbers (order unity) of the large collision time step variants of MPCD do not have a significant adverse effect; however, for more realistic fluids, this is bound to become a concern for all MPCD-based methods. With DPD, good agreement with iFLD viscosity results is obtained at low vol- ume fractions, but the method breaks down at high volume fractions due to the finite size of DPD particles, which is comparable to the inter-colloid separation distance. This can be overcome with a larger number of smaller DPD particles, but the computational advantage of the method is then diminished. Despite the ability to reproduce a higher Schmidt number using DPD, the method overall appears to be a less physically accurate model for colloid suspensions as compared to MPCD.

Computational performance
The computational performance of the different simulation methods is summarized in Table 4. All methods are implemented in the LAMMPS software package [176], as described in Sect. 3. The data in Table 4 are based on runs carried out on the Sandia National Laboratories Red Mesa system, a cluster of Intel Xeon 5500 processors with InfiniBand connectivity. Short equilibrium simulations were carried out for the system described in the preceding section, i.e. 200 colloidal particles interacting with a purely repulsive integrated Lennard-Jones potential. DPD data correspond to the DPD-v1 system, and MPCD data correspond to the MPC-AT collision scheme with a large time step and reverse boundary conditions. While there are some differences in speed for different variants of the MPCD method, they are small compared to differences among the three different solvent models.
As expected, the iFLD method provides by far the fastest performance of all the methods tested, while MPCD and DPD show comparable speed and parallel scaling. However, we note that for both volume fractions shown, there are close to ten times more MPCD solvent particles than DPD particles.
We have used a constant number of colloid particles (200) and varied the size of the simulation box to achieve different volume fractions. This results in the explicit solvent methods being significantly faster at higher volume fractions, where much fewer solvent particles are needed. In contrast, both iFLD and nFLD become slower at higher volume fractions, since the computational effort in these methods is dominated by the calculation of pairwise colloid-colloid interactions; at higher volume fractions, particles are in closer proximity, leading to more frequent near-field pairwise lubrication interactions. In both cases, the nFLD method is significantly slower, due to the matrix problem that must be solved at every time step (see Eq. (18)). As already mentioned, we have used a similar integration time step size for both methods, as the integration time step in this case is limited by the inter-colloid potential. This leads to unphysical diffusion results at short times for the nFLD method, which is predicated on the use of a much larger time step. For the situation considered here, the iFLD method is clearly the better choice; however, for cases where a much softer colloid interaction potential is used, the nFLD method can potentially allow for the use of a time step that is orders of magnitude larger than iFLD; as a result, nFLD can be a much more expedient method in terms of total simulation time achieved, even if each time step requires more wall time to compute.
As previously discussed, a plethora of solvent parameter choices are possible for both the MPCD and DPD methods, which can have significant impacts on computational performance. In this work, several parameters were selected for optimal computational speed; as shown in the previous three subsections, the relatively large size (and consequently lower number density) of DPD particles leads to some undesirable effects, particularly at high colloid volume fractions. For Values are wall times in seconds required to perform 10,000 time steps (lower number means better performance). Numbers in braces below each explicit solvent method in the title row indicate the total numbers of solvent particles in the system more realistic simulations, smaller DPD particles at higher number densities would likely be required, which could lead to a large increase in computational expense. In contrast, the MPCD results do not suggest any need for larger number densities of MPCD particles, except perhaps to achieve higher shear rates without shear thinning of the solvent, and higher Schmidt numbers. Although we have not attempted to do so, it is clear that in order to achieve the accuracy of MPCD using DPD, a much larger number of DPD particles would be required, leading to a drastically higher computational effort. We therefore consider MPCD to be the more computationally expedient method of the two, despite comparable performance data in Table 4. Even for a well-characterized system like the suspension of monodisperse hard spheres treated here, the quantitative prediction of seemingly simple properties such as equilibrium diffusion and shear viscosity is challenging. Based on the results presented, we find the FLD method to offer the best balance between accuracy and computational expedience. Presumably, more rigorous variants of the SD method would provide even more accurate results, albeit with some additional computational expense. However, the flexibility of the different simulation techniques beyond the simple system discussed above must also be considered when evaluating their overall use and future potential. The possibility of introducing more advanced capabilities in these and other similar simulation methods is therefore addressed in detail in the next section.

Advanced capabilities
Newtonian dynamics particle solvers coupled with any approach to hydrodynamic interactions all suffer from at least one fundamental underpinning that prevents general application to practical processing routes for integrating colloidal particles into useful materials. In this section we present the current state and outstanding challenges in what we believe are the key barriers to such general application: (1) non-Newtonian, complex solvent rheology (2) solvent drying/curing (3) non-spherical particles, and (4) complex flow geometries.
Despite these barriers, we want to stress that even in their current state, modeling and simulation tools can be quite useful in addressing practical problems. Processing colloidal dispersions into either highly-ordered films/structures or disordered particle compacts or composites involves complex flows (casting, coating, extrusion, mixing). Regarding the underpinning process flow, much can be learned from the work undertaken and reviewed in Sects. 3 and 4, as nearly all practical processes include colloidal diffusion and shearinduced (or viscometric) flow. Moreover, many systems of interest consist of Newtonian or nearly Newtonian solvent rheology and many colloidal particles are spherical or can be treated as such for practical purposes. Hence, computational studies such as the ones reviewed and taken up in Sect. 4 are useful for determining the effects of inter-particle potentials on flow rheology, volume fraction effects, and microstructural evolution tendencies. Moreover, these approaches can also be used to contrive and/or fit coarse-grained constitutive models for use in larger scale simulations.
Unfortunately solvents in real-world applications of nanocomposite fabrication are often non-Newtonian, and many particles are not spherical. New materials development is most often focused on the solid state (except for lubricants or liquid-crystals), and so drying and solidification is key to most processes. Finally, many processing flows (coating, extrusion, etc.) are not simply bulk shear, but a com-bination of shear and extensional in nature, together with wall-confinement effects. In this section we address these complications and provide some guidance based on existing literature and from our own experience on how best to proceed with existing technology, as well as research challenges that must be addressed.

Non-Newtonian solvents
The non-Newtonian response of a colloidal suspension to impressed shear, even when the solvent is Newtonian, is well known (see Sect. 4). Classic shear thinning is always observed in such systems. Even more extraordinary non-Newtonian responses are observed with large particle loadings and strong colloidal interactions (attractive and/or repulsive), including time-dependent relaxation, viscoelasticity, and even yield-stress behaviors. It is clear that the modeling and simulation community has made great strides in predicting such behaviors [75,91,197]. For processing flows associated with directed colloidal assembly (e.g. the work of Snyder et al. [208]) solvents are usually of low molecular weight (high vapor pressure) and demonstrate a Newtonian (constant viscosity) response to deformation. However, processing flows of this sort remain largely bench-top research, and are in rare instances practiced in industry for large-scale manufacturing. More often, dense suspensions are processed into particle compacts and composites for a number of applications, including the production of energy materials (e.g. batteries, fuel cells), catalytic materials (porous with high surface area), and other structural polymer, metal and ceramic materials. These suspensions are replete with binders and surfactants added for the express purpose of particle interaction control and binding/adhesion (through curing) upon solvent removal. Due to their higher molecular weight, these constituents result in non-Newtonian solvent response, especially during the curing process. Unfortunately, the modeling and simulation methods discussed above are not easily extended to account for such effects.
Extension of the hydrodynamics to accommodate timedependent relaxation and/or viscous response (the classic G and G quantities measured from oscillatory shear experiments [19]) has received the most attention. Clearly, any implicit solvent model based on the Stokes equation for a Newtonian fluid (see Eq. 3) is fundamentally problematic in this respect. This includes all pair-wise and isotropic terms for methods such as SD, FLD and BD. However, simple approximations can be used to account for certain types of non-Newtonian rheological behaviors in these models. Perhaps the simplest approach to this problem is to augment expressions for the hydrodynamic resistance tensor (see Eq. 10) with a frequency dependent solvent viscosity. To this end the equation of motion for a single colloid can be written as follows: The hydrodynamic term is written as a time convolution integral: The kernel κ(t) contains the non-Newtonian behavior of the solvent (note, in Eq. 1 the kernel is assumed to be κ(t) = δ D (t); hence, after evaluating the inegral one is left with a time-independent or quasi-steady drag). Finally the Brownian term must reflect the correlations due to the time dependence of the viscosity: In each of these expressions we have also assumed an FLD-type mean-field-like approximation for the volume fraction dependence of the viscosity via diagonal terms in the R 0 tensor (i.e., [R 0 ] ii ), while ignoring the lubrication terms R δ (see Eq. 9) for simplicity; although the lubrication terms can be included as well. This means that any non-Newtonian effects would not be accounted for in long-range hydrodynamic interactions. This approach of course assumes that the colloid particle size is significantly larger than the long-chain polymers which lead to the time-dependent solvent viscosity, or that, in accordance with assumptions of microrheology, the viscosity that the colloid feels is equivalent to the shear viscosity of the solvent [146]. One example of an algorithm to numerically simulate the so-called Generalized Langevin equation implemented in LAMMPS has been presented by Baczewski and Bond [6] who deployed a frequency (time) dependent viscosity with proper Brownian fluctuation terms (see also references therein).
Extending explicit solvent methods such as MPCD and DPD coupled to DEM colloidal solvers towards non-Newtonian solvents has also been the subject of research. Pryamitsyn and Ganesan [182] extended DPD to account for colloidal particle systems in viscoelastic solvents. Other groups have paid considerable attention to extending DPD to viscoelastic fluid mechanics, largely by connecting a portion of the particles with springs. This can either be interpreted as a direct representation of polymer additives [119,194,212] or as an abstract augmentation of the DPD algorithm that yields viscoelastic solvent behavior [210]. In a similar vein, Tao et al. [214] advanced the MPCD technique presented in Sect. 3.2 to viscoelastic fluids by connecting pairs of MPCD particles with harmonic spring-like potentials. The method is based on alternating streaming and collision steps, just like the Newtonian solvent implementation. They applied this approach to oscillatory shear flow and found the elasticviscous frequency response to be consistent with that of a Maxwell fluid. This work did not address interactions with other colloidal particles. Despite the ability of these types of models to qualitatively capture several features of viscoelastic behavior, no approach that we are aware of has been able to map the non-Newtonian rheology of a realistic solvent and yield quantitative predictions for real-life applications. The scaling and mapping of realistic solvent properties that were of significant concern even for the simple system in Sect. 4 are greatly compounded for these more complex solvent models.
Several groups have extended Stokesian dynamics (SD) as well as the more general boundary element method (BEM, see Sect. 2.2.1) to account for non-Newtonian and viscoelastic solvent behavior. The approach just described for DPD and MPCD, i.e. including an additional discrete element representation of particles with elastic character, was adopted by Binous and Phillips for SD [16,17]. By adding spherical dumbbells connected by finite-extension nonlinear-elastic (FENE) springs to the colloid suspension, they were able to recover the behavior of a Boger viscoelastic fluid-i.e. exhibiting elastic effects, but having a constant shear viscosity. In a different approach, an approximate analytical reformulation of SD to simulate a Maxwell fluid has been presented by Schaink et al. [192]. The more generalized Boundary Element Method has also been extended to viscoelastic solvents by Phan-Thien and Fan using an analytical treatment for an Oldroyd-B fluid [173], which in principle can be extended to other rheological models, but is limited to unbounded domains (or possibly domains with simple boundary conditions). A similar approach has been advanced in the continuum flow arena, viz. to use a microstate equation in the continuum which carries the response of elastic molecules to solvent deformation [98]. The continuum approach typically must be mesh-based (e.g. finite element, finite volume methods), or posed in a Lagrangian-particle framework, as is discussed below.
Immersed Boundary Methods (IBM) as recently reviewed by Lechman et al. [134] offer the most direct route to non-Newtonian solvent rheology, as the added, but expensive, benefit of a background or body-fitted mesh provides a route for several decades of mesh-based algorithm development to address non-Newtonian effects. Halin et al. [88] advanced the so-called Langrangian particle method for computing timedependent viscoelastic flows. Their method deploys either a differential constitutive equation (macroscopic approach) or a kinetic theory model (micro-macro approach) on a background mesh for the viscoelastic stress. The discrete Langrangian particles carry the stress state and its history through the interpolation, and they can serve in a dual role as suspension particles. Hwang et al. [99] pioneered an approach based on earlier work of Baaijens [5] to couple discrete particle motion with a uniform background FEM model of incom-pressible viscoelastic flow. They deployed Lagrange multiplier constraints to match the particle-fluid stress and solventmass-displacement boundary conditions on the particles. However, their application was restricted to two dimensions with no apparent effort to extend it to three. A recent method developed by Noble et al. [163] known as the conformaldecomposition finite element method (CDFEM) may offer the most general framework to address non-Newtonian solvent effects. They have made considerable progress towards scalability to parallel platforms and hence may be well positioned to incorporate such effects in flows of colloidal suspensions. A more recent IBM-like approach was advanced by Chrispell and Fauci [40] and applied to a complex flow geometry (refer to Sect. 5.4). Using a finite-difference, markerin-cell Navier-Stokes solver they developed a method to integrate an Oldroyd-B constitutive equation. They successfully applied their approach to peristaltic pumping of colloidal suspensions at reasonably high Weisenberg numbers. In summary, IBM methods coupled with colloidal DEM offer a straightforward route to solve this problem, but with the requirement of at least 100-1,000 particles in a threedimensional meso-scale simulation, computational expediency remains an outstanding issue with this approach.
Accounting for non-Newtonian solvents in colloidal dynamics solvers is clearly still an outstanding multiscale challenge. From the most general, straightforward approaches based on coupled FEM/FDM and DEM to the more esoteric extensions of FLD, SRD, and DPD, the research challenges seem not surprisingly to lead to the need for computational efficiency and rheological accuracy. FLD, with all of its merits vis-à-vis DPD and SRD, will face the same potentially unsolvable challenges, mainly due to its quasi-analytical origins, a potential pitfall that extends equally to SD and BEM.

Non-spherical particles
Colloidal particles are often spherical due to the manner in which they are synthesized. In both high-temperature gasphase reactions and solution-based nucleation growth, spherical shapes are thermodynamically favorable. As already discussed, the spherical particle assumption mathematically enables and is in fact the foundation of most DEM approaches that include contact and colloidal forces. On the other hand, particles larger than colloids are commonly aspherical as fabrication routes in this regime are often based on pulverization and milling. Larger particles made by spraying/atomization and drying also often end up ellipsoidal due to drying stresses. In the colloidal size regime, however, recent topical interest in so-called nano-materials has drawn attention to processing particles like nanotubes (e.g. carbon nanotubes), nanowires (e.g. silicon) and graphene flakes [136,198]  ated hydrodynamics to accommodate such shapes has come to the forefront. In addition to nano-materials, some creative solution-based approaches to producing mildly aspherical shapes like di-spherical colloids, or di-colloids, have emerged [111]. However, colloidal particles of this sort are rarely encountered in large-scale manufacturing, perhaps because the process of making them has not been scaled up. Highly aspherical colloidal particles can also be fabricated with nano-imprinting processes. Caldorera-Moore et al.
[32] used imprint lithography to make a variety of shapes of hydrogel particles. Again, such approaches hold promise for small batches of particles used for applications like drug delivery, but they represent a specialized case visà-vis modeling requirements for concentrated dispersions used in nanocomposites. While the development of modeling approaches that address generalized particle shapes is an important challenge to be met in the broader DEM space, most colloidal applications today deal with easily parameterized shapes, such as spheres, cylinders, and plates.
There are a number of ways to construct and parameterize models of non-spherical particles, as illustrated in Fig. 8. To first order, arbitrary shapes such as a di-colloid can be approximated with dumbbell assemblies, or if a smooth, single particle is required with a spheroid, as shown in Fig. 8a, d. Clearly, if the goal is to represent a hexahedral body, these simple shapes do not capture important features, such as sharp corners and flat surfaces, which can have a significant effect on microstructure and dynamics. A more generic and accurate representation can be accomplished with collections of spheres (overlapping or non-overlapping) [81,178], which are then constrained to move as rigid bodies, as in Fig. 8b, c. This generalization is also applicable to cylinders, rods and sheet-like objects. The corrugation that results from this approach is sometimes undesirable, as it can lead to artificial inter-particle stacking at close range. The most generic approach involves representing a shape using a surface triangular mesh, with the possibility of an interior tetrahedral volume mesh, as shown in Fig. 8e, f. The LAMMPS software package [176] provides the basic computational infrastructure to accommodate all of these representations. What remains a challenge, however, is a workflow which allows for building a collection of such non-spherical particles to be used for a mesoscale simulation. Several helpful algorithms for generating and using libraries of arbitrary shapes have been presented in literature [81,133,178], but the process is highly dependent on the problem at hand, and typically requires some effort on the part of the analyst.
Regardless of the representation, extensions for accommodating aspherical particles are difficult to implement in most simulation methods for two reasons. First, the method itself may be predicated on the spherical particle assumption, as is the case in quasi-analytical hydrodynamic treatments underpinning Brownian Dynamics, Stokesian Dynamics and FLD. Second, soft/long-range colloidal forces are difficult to reconcile without generalized and expensive surface-tosurface distance calculations required to determine the pairwise force. Significant work has been pursued to address both of these areas, as we discuss below.
Hydrodynamic interactions in implicit solvents pose a significant challenge. The Stokesian Dynamics (SD) method and its expedients, such as FLD, are limited to spheres due to their analytical underpinnings. However, dicolloid shapes in suspensions have been modeled with SD-related extensions [125,127]. Meng and Higdon [153,154] have generalized to plate-like particles using an extension to SD that approximates the hydrodynamics based on a planar assemblage of hard-spheres (viz. no colloidal interaction). An extension of their work added Brownian motion to the model [154]. They reported the effects of Peclet number (shear rate) and volume fraction on the effective suspension viscosity. Additionally, SD has been extended to prolate spheroids [41,42] and fibres [30,190]. Nevertheless, the SD approach and all related methods are fundamentally limited to relatively simple particle shapes that are amenable to analytical treatment. Unfortunately, generalizations of these methods, like the Boundary Element Method (BEM), which easily accommodate nonspherical shapes, are often too expensive to be practical for large systems. On the other hand, explicit solvents such as MPCD and DPD can in principle be used in conjunction with any arbitrary particle shapes. The challenge there lies in developing efficient solvent-colloid collision detection and surface interaction algorithms that enforce no-slip boundary conditions at the particle surface. This can largely be resolved for simple flow geometries with MPCD [24], and significant effort has been expended for DPD [58,174,185]. Further investigations for more complex flow situations involving the various particle representations discussed above are warranted for these methods, but the difficulties appear to be manageable.
Two alternative modeling routes exist that are feasible for general application to aspherical particles. First, atomistic models have no such limitations regarding particle shape, as both solvent and colloid can be represented at the atomistic scale. Presuming accurate atomistic potentials are available for all atom-atom (solvent and particle) interactions, addressing nonspherical colloidal systems is feasible [115] but limited to just a handful of colloidal particles. In short, an atomistic approach is currently too computationally intensive to address mesoscale systems of particles in a suspension at large enough length scales to predict bulk properties. At the other extreme, a second approach with no such particle shape limitations would be direct continuum formulation for particle-solvent interactions coupled with grid/elementbased immersed boundary methods (cf. Sect. 5.1). In this case the mechanics of each particle can be handled in the con-tinuum, i.e. particle motion and deformation [162] or with a discrete element approach [134]. Particle-solvent interactions can be accommodated with a variety of approaches as recently reviewed by Noble et al. [163] Beyond atomistic approaches, no such work was found on coarse-grained mesoscale modeling of graphene flake suspensions or carbon nanotubes. The major challenge here would be the approximate long-range inter-particle forces (which may be modulated by the presence of polymers/additives, in addition to their hydrodynamic interactions), as well as obtaining accurate representations of the shapes (triangular meshes or composite spheres). Because of the simpler shapes, solving a microstate orientation equation in the continuum together with the appropriate constitutive relations in a FEM/FDM framework is another option [141,142], but detailed mechanisms of aggregation/agglomeration cannot be obtained in this approach.
Perhaps the most significant challenge pertains to models for the proper colloidal interaction potentials for nonspherical particles. More specifically, the challenge is to develop integrated colloidal potentials that can be efficiently and accurately applied to coarse-grained colloidal model systems. For mildly aspherical particles like spheroids, several colloidal potentials are available [68]. For composite spherical representations (see Fig. 8c), long-range inter-particle interactions can in principle be treated as a sum of pairwise interactions among constituent spheres, but the physical relevance of such an approach is questionable, especially for cases of highly overlapping spheres. On the other hand, short-range (granular) contact potentials can be readily implemented in this context [81,178].
While several atomistic studies addressing this challenge have been undertaken [93], the best and most accurate approach remains to be determined. in't Veld et al. [221] recently compared pair-wise interactions in atomistic simulations of composite nanoparticles with integrated spherical forms. They compared pairwise forces between colloidal particles composed of smaller aggregated LJ atoms to the integrated case due to Everaers [68]. They determined that simply computing the interaction via the Lennard-Jones atoms on the surface of the composite colloidal particle led to incorrect temperature/vapor-pressure behavior, unlike the properly integrated case. Similar results were observed for composite nanoparticles, with the atom-atom interactions truncated at a finite distance. In summary, while DLVO and related potentials developed for spherical systems can be extended with minor modifications for mildly aspherical systems (spheroids, dicolloids, etc.), much work remains for generic shapes. A promising alternative yet to be explored in great detail may be the surface element integration method [15] which addresses in detail the colloidal interactions between spheres and plate-like particles. Fig. 9 Phases of coating and drying of colloidal suspensions: a Convective assembly process for highly ordered films, b fast drying and fully densified films lead to c capillary stresses

Drying and solidification
Addressing the underlying mechanics of processing colloidal suspensions into functional materials (films, composites, fibers) requires not just a firm understanding of rheological behavior of the suspension, as discussed in Sect. 5.1, but also the stability of the suspension under processing and during solidification. Solidification typically starts with volume reduction through drying. Two distinct processing routes are of technological relevance here: the first is aimed at the production of highly-ordered films (monolayers or multiple layers) through what is known as colloidal directed assembly, viz. self assembly influenced by external forces. Although a number of approaches have been explored that deploy electromagnetic forces, the most noteworthy and scalable route is so-called convective-assembly/drying [208,229], which deploys a metering blade to apply the film and subsequent drying to assemble the ordered layer(s) of particles (see Fig. 9). In this case, the initial suspension is highly dilute in order to achieve sufficient particle mobility. Unfortunately, long-range ordering is prone to thermally-induced defects, and so this process is typically run slowly (of order 1 cm/min in typical coating processes). The second distinct drying/processing route is often used in the production of disordered nanocomposites. Dispersions typically being cast into such materials are highly loaded (of order 20-40 % or higher by volume) in order to reduce energy requirements and to force the microstructure to be amorphous. At high particle concentrations, long-range ordering is not the goal. Typically a plethora of high-speed coating and drying processes can be brought to bear. Managing suspension rheology is paramount to successful processing, and subsequent drying and curing is often fraught with defects (residual stress and cracking).
Modeling and simulation tools at the continuum scale have been central to process design of coating and drying for decades [179,196]. Colloidal particle concentration tracking during drying of drops and films has been treated with suspension balance and population balance approaches at the continuum scale [29,156]. However, this work is aimed at understanding the connection of processing to microstructure, which is best achieved with meso-scale models of the sort we have discussed above. At this scale, modeling the effects of a true drying process is still in need of attention.
Foremost are the challenges associated with model formulations that best represent an actual drying process. If the concern is basically the physics underpinning uniform, bulk volume reduction, then this can be achieved with a uniform compression of the simulation box in a meso-scale model, using the usual periodic boundary conditions in all directions. For implicit solvent methods like FLD, these simulations are straightforward without any significant modeling advances, and essentially mimic the case in which solvent is removed from the volume in a uniform way. If particle diffusion is much faster than the volume reduction rate, this simple simulation does represent a slow drying regime, viz. one in which the diffusion time scale τ D = a 2 /D 0 is much less than a/v int , where v int represents the wall speed and a the particle radius. With explicit solvent methods, the solvent particles at the boundaries must be deleted from the system in a way that maintains thermodynamic consistency. These considerations can affect the temperature of the system. Some consideration should be made of latent heat consumption (cooling) during drying, as some workers have addressed at the atomistic [38] and meso-scales [184].
Actual drying or solvent removal drives the predominant solvent and particle transport in one direction relative to the drying interface. Whether the medium is a drop, fiber or film, drying-induced microstructural formation is still a threedimensional problem and ripe for exploration with mesoscale models (see Fig. 9). Special boundary conditions on the dry-ing interface must be contrived for such simulations. Boundaries normal to the interface can be treated as periodic, but at the interface solvent must be removed in a thermodynamically consistent way, depending on the drying regime. We identify two distinct drying regimes: -"Slow drying"-particle diffusion rate much faster than the drying rate -"Fast drying"-drying rate much faster than the particle diffusion rate Within the "fast drying" regime of colloidal systems, as pictured in Fig. 9, there are several phases of the process with distinct physical descriptions: -Constant rate regime-low volume fraction of colloids, which do not impede diffusion path for solvent significantly -Falling rate regime-particle concentration increases at the drying front (free surface), impeding solvent mobility and slowing evaporation -Percolated network-solvent recedes in particle space and drying controlled by pore-diffusion/flow of solvent In a slow-drying regime the rate of drying continually falls, but less precipitously, as no surface layer ("skin") forms. When drying is slow, simply moving a single boundary along its normal direction at the desired drying speed for implicit solvent methods, and/or deleting solvent particles in explicit methods commensurately with the moving boundary, is useful for exploring the microstructural implications. In the case of atomistic explicit solvents such as LJ or atomistically represented water, particles naturally form a liquid-vapor interface, and move from the liquid to the vapor phase; the evaporation rate can be controlled by removing particles from the vapor phase at various rates [38]. If drying is "fast", then the effect on the drying rate of surface-microstructure (skin) formation must be accounted for. Clearly, explicit solvents should capture this effect: evaporation of particles from the surface, or boundary, will induce diffusion of solvent particles to feed the evaporation. That diffusion, if fast enough, leads to heightened convection of colloids towards the surface, induces tension in the liquid, and ultimately leads to a compression of particles near the surface. No known mesoscale approach deploying an implicit solvent (e.g. FLD) can capture this effect due to their quasi-analytical nature. Challenges also remain for coarse-grained explicit solvent approaches related to local depletion in regions where colloids are closed-packed. Simulations of drying processes clearly require additional effort in algorithmic development.
Some relevant work in this area is noteworthy. Interfaces between dense (liquid) and gaseous regions can be modeled with coarse-grained explicit solvents, provided that some energetic interaction among solvent particles exists. If additional energetic interactions between solvent molecules and colloids are present, these approaches can also account for capillary interactions among colloids. However, colloidsolvent interactions must be adjusted appropriately in order to reproduce realistic capillary effects, and the choice of parameters is not straightforward. Most of this work deploys lattice-Boltzmann solvents [112,113,199], but is typically limited to small systems. For models where colloid-solvent energetic interactions are not present (such as MPCD or FLD, as well as some variants of DPD), capillary effects will not arise spontaneously. Instead, additional capillary forces must be applied to particles at the interface (in the case of implicit solvent, no true interface exists, but it can be simulated for the slow drying regime with a moving flat wall, as described above). These capillary forces are based on approximate analytical solutions to the Young-Laplace equation, and have been studied extensively by Kralchevsky and coworkers [47,[120][121][122] and more recently by Vassileva et al. [220]. The resulting expressions have been used in analytical equilibrium models [135] as well as dynamic simulation methods [79]. In the latter case, Fujita and Yamaguchi used an immersed boundary method to couple the resulting DEM to a Navier-Stokes solver, a rigorous but computationally expensive approach. They do however capture the dryinginduced convective effects critical in fast-drying regimes.
As drying proceeds to a complete percolated, stresssupporting network, evaporation continues as the solvent recedes into the film (Fig. 9b). Capillary stresses can further induce consolidation at this point, and in some cases particles can deform and even crack. At a larger scale the film can form defects (mud-cracking). Simulating this problem at the mesoscale requires granular contact potentials (e.g. Hookean/Hertzian with friction) and possibly some accommodation for particle deformation and even grain-boundary diffusion (at high temperature). Theoretical treatments in the continuum have been reported, such as the quasi-analytical approach of Singh and Tirumkudulu [204]. Additionally, network approaches with statistically-based bond-breaking probabilities [149,206] are limited to simplified phenomenological models.
To make matters more complicated, in most practical applications, drying for composite material production is accompanied by the curing of a binder. We found no relevant literature in which this effect was accommodated in mesoscale simulations. With coarse-grained explicit solvent approaches, the chemistry underpinning curing would be difficult to implement, as it is often based on free-radical or condensation reactions, which would require the tracking of species concentration, etc. One could implement some sort of viscosity increase through a thermal effect, which would account for some of the dynamics of curing during the consolidation of particles. This is readily implemented in implicit solvent models (including SD, BD, FLD), so long as a suitable rate model for the viscosity is available. Accounting for curing effects, including their contribution to the final residual stress of a film, remains an outstanding research challenge for DEM mesoscale modeling.

Variable complex flow geometry
Up to this point we have only addressed the applicability and extensibility of mesoscale DEM modeling in simple geometries, with very little discussion of their application to flows not simply driven by Brownian motion or impressed simple shear in simulation boxes with obvious periodicity. The drying problems discussed in Sect. 5.3 present the first need for non-standard boundary conditions (non-periodic), and it is clear that much work is needed, especially with respect to the quasi-analytical, implicit solvent formulations like FLD. Extending the methods highlighted in this paper to accommodate no-slip walls, non-periodic boundaries that account for drying fronts, or flow regimes that are not simple bulk shearing is of great interest for other applications as well.
As with the other advanced capabilities addressed herein, a general approach with no restrictions and complete generality in complex geometries is to solve the Navier-Stokes equations in the actual application of interest. With meshbased methods like finite element, or lattice-particle methods (LB) and even with the explicit solvent methods discussed earlier, such simulations in principle can be carried out, so long as wall interactions with the solvent particles can be designed to produce the desired boundary behavior (e.g. no slip). That is, there are no inherent limitations to imposing no-slip boundaries or other sources of impressed flow (pressure driven flow, stretching flow, etc.). Recent work by Zhao et al. [237] used SRD and DPD successfully in modeling electro-osmotic flow in a microfluidic cell. The same methods have also been recently extended to complex geometries and free surfaces between two phases (like capillary free surfaces), as motivated by related work using lattice Boltzmann techniques [1,112,113]. Although implicit solvent models such as SD and its expedients can be extended to include the effects of planar walls [203], more general boundaries cannot easily be treated, and require significant theoretical and algorithmic modifications.
Regarding generalized Navier-Stokes solvers coupled with DEM particle solvers, the work of Sasic et al. [191] is noteworthy. Building on a finite volume solver and the volume-of-fluid multiphase flow model, similar to the finite element based approach of Baaijens [5], Sasic et al. advance a method to couple a particle solver using the volume-offluid marker function to imprint the volume displacement on the solvent mesh. The method is similar to a multiphase flow approach, but seems to be limited to a small number of particles. In this same generalized framework, Chrispell and Fauci [40] used finite difference methods and markerand-cell free boundary tracking. In general, these grid-based approaches are computationally intensive when coupled with DEM, and thus limited to small systems until work on scalability to massively parallel platforms is undertaken. However, as with all such FEM, FDM, or FVM solvers, these approaches can accommodate all boundary conditions.
What remains an outstanding challenge in the rheological arena is the simulation of mixed shear and extensional deformations, or even pure extensional deformation. No obvious periodic simulation domains exist to impress a purely extensional flow. However, the work of Kraynik and Reinelt [123] advanced a clever method to do just this for the flow of foams, even though it seems to only recently have been discovered for use in atomistic simulations. Their work solves the latticecompatibility condition for planar extension. While straightforward for simple shearing flow, this condition is subtle for extensional flows. At the time of this writing no such work in the meso-scale particle flow arena has taken advantage of this work.
In conclusion on this topic, the most pressing challenge, which may in fact not have an elegant solution, is the extension of implicit solvent methods such as SD and FLD to boundary and physical conditions often encountered in real-life applications. Specifically, because these methods are quasi-analytical, extension of the underpinning theory to accommodate long-range hydrodynamic interactions between particles and walls, and particles and free surfaces is difficult and may only be accessible with mesh-based or particle-based solvents. With regard to explicit solvent methods such as MPCD and DPD, some challenges remain in selecting surface-solvent interactions to enforce the desired boundary conditions, but in principle arbitrary flow geometries can be simulated.

Conclusions
We have presented an overview of particle-based mesoscale simulation techniques for colloidal suspensions. In Sect. 4, we focused on a quantitative comparison of three commonly used techniques: FLD, an implicit solvent method that represents a significant simplification of Stokesian dynamics, as well as DPD and MPCD, two distinct explicit solvent approaches to simulating coarse-grained fluids. We deployed each method to model thermal motion (diffusion) and viscometric flow (shear) of suspensions with a prescribed solvent viscosity. Our approach stands in contrast to some previous works, in which only key non-dimensional parameters were reproduced, but various specific properties and parameters (e.g. viscosity and particle diameter) were not explicitly controlled or mapped to a particular solvent [11,22,36,169,181,217]. While this latter approach would likely have been adequate for the simple diffusion and bulk shearing simulations presented in Sect. 4, the exact mapping of solvent properties is an essential step in simulations of more realistic, complex systems.
With FLD, both the solvent viscosity and particle size are direct inputs to the simulation, thus making for a trivial mapping procedure. With MPCD, the availability of analytical expressions greatly simplifies the selection of MPCD fluid parameters, and the particle size is in principle a direct simulation input for the colloid-solvent coupling schemes considered here; however, in practice, it turns out that both colloid diffusivity and viscosity are sensitive to many algorithmic details of the MPCD method, and no clear choice can be identified to produce the desired results. With DPD, solvent properties cannot be accurately predicted a priori for a given set of parameters, and must be measured from simulations of the pure DPD fluid. Furthermore, no clear choice or physical basis for colloid-solvent coupling exists, and ad hoc potentials must be introduced, with parameters that can only be selected based on a trial and error approach. While this is feasible for the systems considered in this work, it becomes unwieldy for more complex cases (e.g. polydispersity, nonspherical particles, etc.), and any parameters thus selected do not transfer readily to different situations. As such, we conclude that FLD and similar Stokesian dynamics-based methods are in general superior to explicit solvent methods with regard to mapping of key physical parameters. We identify this as a notable shortcoming of both MPCD and DPD in the simulation of colloid suspensions, and an area that requires further method development.
With regard to accuracy, the FLD method is once again superior as measured by the accuracy of the early-and latetime diffusivities at various volume fractions. FLD results are in excellent agreement with theoretical and experimental predictions, whereas notable deviations are apparent in both MPCD and DPD data. With MPCD, most variants and parameters tested yield values that are in reasonable agreement with the expected diffusivities, but no clear picture emerges as to systematic effects of different parameters. With DPD, parameters can in principle be selected to match the target diffusivity data at low volume fractions, but the method fails at higher volume fractions. MPCD and FLD viscosity data are in excellent agreement, but MPCD is limited to moderate shear rates; while DPD performs well with regard to viscosity at low volume fractions, it once again fails at high volume fractions due to the finite size of DPD particles.
Given the superior accuracy and performance of SD-based methods such as FLD over explicit solvent methods, it would seem that the former is a clear choice for colloidal suspension simulations. However, the merits of explicit solvent methods become more apparent when moving away from the canonical system of monodisperse spheres in an infinite domain to more complex, realistic systems, such as those discussed in Sect. 5. Although the difficulties noted for the simple systems persist in these cases, there is significant flexibility to be gained through the use of explicit solvents. In particular, modeling hydrodynamics for non-spherical particles and complex flow geometries is relatively straightforward with explicit solvent methods, whereas only a limited number of simple shapes and geometries can be treated using SD-based implicit solvent methods. Explicit solvents extend more readily to models of solvent evaporation, although simple models of curing can be more readily implemented in the context of FLD. Finally, modeling solvents with non-Newtonian viscosity appears to be more mature for MPCD and DPD, although much work is needed for all methods to achieve realistic complex solvent rheologies.
Overall this work provides a direct comparison of several particle-based simulation methods for colloid suspensions. For simple systems, Stokesian dynamics-like methods, such as FLD, are the clear method of choice; presumably, more rigorous variants of the SD method yield even better results. We have quantified a number of shortcomings for MPCD and DPD explicit solvent methods even when applied to a suspension of monodisperse spherical particles. We feel that these deficiencies can be resolved, and that the potential for applications to realistic systems is promising. While greater flexibility can be attained using direct numerical simulations of solvent hydrodynamics (e.g. finite element solutions of relevant continuum equations), the computational simplicity and ease of parallelization for particle methods cannot be matched. In general, we believe particle-based methods such as those discussed in this work will remain an important tool for the simulation of colloidal suspensions.