A new model for simulating colloidal dynamics

We present a new hybrid lattice-Boltzmann and Langevin molecular dynamics scheme for simulating the dynamics of suspensions of spherical colloidal particles. The solvent is modeled on the level of the lattice-Boltzmann method while the molecular dynamics is done for the solute. The coupling between the two is implemented through a frictional force acting both on the solvent and on the solute, which depends on the relative velocity. A spherical colloidal particle is represented by interaction sites at its surface. We demonstrate that this scheme quantitatively reproduces the translational and rotational diffusion of a neutral spherical particle in a liquid and show preliminary results for a charged spherical particle. We argue that this method is especially advantageous in the case of charged colloids.


Introduction
Understanding the dynamics of colloidal dispersions has been a long-standing problem in condensed matter physics. Despite the continuous need for accurate theoretical predictions concerning particle mobilities or sedimentation rates, the development of the field meets a nearly unpassable barrier. The many-body character of the hydrodynamic interactions (HI) between the colloidal particles (i. e. their correlated motion as a result of solvent-mediated fast momentum transport) poses a serious challenge for analytic theory. On the other hand, the rapid growth of available computer power in the last decades has made it more and more feasible to descend to a more basic description of colloidal systems rather than trying to develop an adequate macroscopic description. Such a step was made recently for the statics of charged colloids, where simulations of the primitive electrolyte model, taking the particle character of the charges and fluctuations in the charge distribution fully into acount, resolved some important strong correlation issues and thus have resulted in a significant refinement of the existing theories of screening [1]. In a similar fashion, many-body effects in the HI of suspensions of only moderate density are quite important. Although their analytical form is, in principle, known in terms of a multipole expansion [2], it is numerically very cumbersome to take these higher-order many-body terms into account in a Brownian or Stokesian dynamics [3] simulation, where only the positions of the colloidal particles enter. For this reason, it is more convenient and for large particle numbers also more efficient, to take these effects into account by simulating the solvent degrees of freedom, and, in particular, the momentum transport through the solvent explicitly. An additional bonus is that retardation effects, if present, are automatically taken into account.
Standard Molecular Dynamics (MD), where just Newton's equations of motion are integrated numerically, is of course able to simulate HI correctly. However, it is very time-consuming to simulate the motion of a huge number of particles explicitly in such detail. The time step is governed by the local oscillations of the solvent particles in their temporary "cages", which is much faster than the motion of the colloidal particles. This information is however not needed for correctly reproducing HI. One just needs an appropriate mechanism for momentum transfer on somewhat larger time scales (however still small compared to the motion of the colloidal particles). This can be done in various ways: Dissipative Particle Dynamics (DPD) [4] is essentially MD, where however the particles are made very soft in order to increase the time step, and where a momentum-conserving Langevin thermostat is added [5]. Another possibility is multi-particle collision dynamics [6], where the solvent is modeled as an ideal gas, and collisions are modeled as simple stochastic updating rules which conserve energy and momentum. Grid-based methods can be either the direct solution of the Navier-Stokes equation by a finite difference method, or the Lattice Boltzmann (LB) [7] approach, where essentially a linearized Boltzmann equation is solved in a fully discretized version (i. e. space, time, and velocities are discretized).
In our opinion, all of these approaches have their advantages and disadvantages.
Nevertheless, it seems that all of them are roughly equivalent with respect to their computational efficiency, and with respect to the level of accuracy of describing the physical phenomena. Therefore, we believe that the choice of method is mainly a matter of convenience. For our model we have chosen the LB route, mainly because this algorithm is particularly simple to implement and parallelize. Moreover, we were inspired by previous highly successful simulations of hard-sphere colloidal suspensions [8,9,10,11,12,13,14,15,16]. A further advantage of a grid-based method is that thermal fluctuations can be both turned on and off, depending on the physical situation under consideration (in a particle method they are always present). This approach is essentially a hybrid method, where the solvent is modeled via LB, while the colloidal particles are run with MD.
The original approach by Ladd [8,9,10,11,12] models the colloidal particles as extended hollow spheres, while stick boundary conditions at the surface are implemented (roughly spoken) via bounce-back collision rules. One should note that for moving spheres this involves some minor fluctuations in the fluid mass, which is included within a sphere: The moving shell incorporates some new fluid at the front, while it releases mass at the rear. The fluctuations are a natural effect of the thermal density fluctuations of the fluid, and of the lattice discretization. Using this method, the translational and the rotational dynamics of colloidal spheres were accurately described.
This model has been recently extended to the case of charged systems [17], where the colloidal particles carry a charge, while the counterions and salt ions are taken into account as LB populations, such that the electrostatics is essentially taken into account on the level of the Poisson-Boltzmann equation. This method has two disadvantages: firstly, the discrete nature of the ions and correlations beyond the Poisson-Boltzmann level are not taken into account; secondly, one cannot avoid a leakage of charge into (and out of) the sphere (as in the case of mass), such that it is hard (if not impossible) to maintain a well-defined Debye layer of charges around it [18]. For these reasons, it is advisable to take the counterions explicitly into account. The disadvantage, however, is that only rather small size ratios (size of colloidal particle vs. size of ion) are easily accessible. Nevertheless, we believe it is useful to study such a system, with respect to both its equilibrium and its nonequilibrium properties. The purpose of the present paper is to describe first steps in the development and validation of the corresponding model.
For the coupling of the small ions to the LB hydrodynamics, one would like to simply model the former as point particles. Fortunately, such a coupling has been recently developed by Ahlrichs and Dünweg [19,20] with the purpose of studying the dynamics of polymer solutions [21]. Each Brownian point particle is assigned a phenomenological friction coefficient ζ, and the coupling is implemented just as a dissipative force F = −ζ( V − u) acting on the particle, where V is the particle velocity, while u is the solvent flow velocity at the particle's position, obtained via linear interpolation from the surrounding lattice sites. Furthermore, the LB variables at these sites are adjusted to ensure momentum balance, and thermal fluctuations are added to both types of degrees of freedom. For more details on implementation and validation of this method, see Refs. [19,20].
It is then convenient to use this coupling not only for the ions, but also for the colloidal particles, whose larger size is taken into account by modeling them as some arrangement of interaction sites. For reasons of efficiency, we take only points at the surface of the colloid. Furthermore, the coupling constant ζ is chosen rather large, which approximates stick boundary conditions. It should be noted that the dissipative nature of the coupling does not model any squeezing-out of solvent, hence the sphere is filled with the same amount of solvent regardless if it is hollow or if it were filled with particles. The solvent within the sphere follows the motion of the surrounding shell (with some minor time lag, see below), and therefore one can view the solvent inside the sphere as just belonging to the colloidal particle. Adding further particles in the sphere's volume would have no effect except coupling the fluid within the sphere even tighter to its motion. Within our desired level of accuracy, this turned out not to be necessary. The long-time motion of the sphere is the same as that of a corresponding hard sphere, as we will show in the present paper.
We analyze basic dynamic properties of such a model colloid, where we restrict ourselves to the case of a single sphere. A detailed analysis is done for the neutral case, while some preliminary data for the case of a charged sphere are presented. We show that the (neutral) model exhibits the essential dynamical features of a spherical colloidal particle in a liquid. In our view, our method provides comparable efficiency to Ladd's approach, while being quite straightforward to implement. Furthermore, it provides substantial flexibility with respect to the properties of the colloidal surface, namely, deformable, permeable, and non-stick surfaces can be easily simulated.
The remainder of this article is organized as follows: In Sec. 2, we describe our simulation model, while Sec. 3 contains the numerical results on translational and rotational diffusion. Finally, Sec. 4 concludes with a brief summary.

Model
Our hybrid simulation method involves two subsystems: the solvent that is modeled via LB with fluctuating stress tensor (i. e. we run a constant-temperature version of the LB method) and a Langevin MD simulation for the particles immersed in the solvent. The LB simulation is performed using the 18-velocity model [9], using the protocol described in [19,20]. The fluid simulation consists of collision and propagation steps, the former being performed with inclusion of the momentum transfer from the solute particles (surface beads, and, for charged systems, ions).
The colloidal particle is represented by a two-dimensional tethered bead-spring network consisting of 100 beads, which is wrapped around a ball of a radius σ cs (for notation, see below), so that the whole construction resembles a raspberry (see Fig. 1). The network connectivity is maintained via finitely extendable nonlinear elastic (FENE) where k is the spring constant, and R 0 the maximum bond extension. Furthermore, the beads repel each other by a modified Lennard-Jones (LJ) potential An additional repulsive LJ bead is introduced at the center of the sphere in order to maintain its shape. In Eq. 2, i, j denote either a central ("c") or a surface ("s") bead.
The unit system is completely defined by the surface bead parameters by setting ǫ ss , σ ss , and the surface bead mass m s to unity. The interaction between the central bead and the surface beads is described by σ cs = 3, which is thus the sphere radius, and ǫ cs = 8. Furthermore, the FENE spring constant for the surface beads is k = 300 and the maximum bond extension is R 0 = 1.25. To simulate a charged colloidal particle, we place the charge at the central bead, and add an appropriate number of counterions (LJ beads with "s" properties) outside the sphere. The electrostatic interaction is taken into account via the Coulomb potential between the various charges, where the standard Ewald summation technique [22] is applied. In Eq. 3, λ B = e 2 / (4πε 0 ε r k B T ) is the Bjerrum length, k B the Boltzmann constant, q i the charge of species i in units of the elementary charge e, and T the temperature.
The LB lattice constant is chosen as one (in our LJ unit system), and the fluid is simulated in a cubic box with periodic boundary conditions. The force between the LB fluid and the surface beads (or ions) is given by Here, ζ is the "bare" [20] friction coefficient, V and u are the velocities of the bead and the fluid (at the position of the bead), respectively, while f is a Gaussian white noise force with zero mean, whose strength is given via the standard fluctuation-dissipation theorem [19,20] to keep the surface beads and ions at the same temperature as the solvent. The central bead is not coupled to the solvent (here ζ = f = 0); as discussed in the Introduction, the behavior of the model would change only marginally if such a coupling were included. In our simulation we used a friction constant ζ = 20, a temperature k B T = 1, a fluid mass density ρ = 0.85, and a kinematic viscosity ν = 3, resulting in a dynamic viscosity η = 2.55. At least 20000 MD steps were performed to equilibrate the initial random bead configuration before the interaction with the LB solvent was turned on. Further details on the method can be found in [19,20].

Results
We test the simulation method against basic relations for an isolated sphere in solvent. First, we look at the center of mass' velocity relaxation. The simplest experiment to perform is a "kick". The sphere is placed in a LB fluid at rest (i. e. without thermal noise), and at time t = 0 all particles of the sphere are assigned an identical velocity V = 1 in x direction. Fig. 2 monitors the time behavior of the sphere's center of mass velocity, normalized by the initial value. According to linear response theory, this relaxation function must be identical to the normalized center-of-mass velocity autocorrelation function for Brownian motion in thermal equilibrium, if the initial kick is weak enough. This is indeed satisfied, as a comparison of the two curves in Fig. 2 shows. For the experiment in thermal equilibrium, we performed 10 runs with different random number generator initializations, in order to reduce the statistical uncertainty. It is well-known that simulations of Brownian motion in a hydrodynamic solvent are always strongly affected by finite size effects. The diffusion constant, and therefore also the relaxation function, depend on the linear system size L due to hydrodynamic interactions with the periodic images. The diffusion constant exhibits a finite-size correction of order R/L [20,23]. Asymptotic behavior can therefore only be expected for R/L ≪ 1, and this is why we performed the experiment in a rather large box of size L = 80. For the same reason, the equivalence between "kick" experiment and Brownian motion will only hold if the comparison is done for the same box sizes.   One clearly sees that the initial decay of the relaxation function is characterized by two relaxation processes, one initial fast decay followed by a somewhat slower relaxation. Qualitatively, this may be explained as follows: A compact sphere of radius R and mass M should exhibit a velocity relaxation which, after transient ballistic motion, is initially characterized by an exponential decay exp (−t/τ ), with a relaxation time τ = M/ζ tot . Here, ζ tot is the total friction coefficient, which we estimate via Stokes' law for stick boundary conditions as ζ tot ≈ 6πηR ≈ 144. However, the effective mass is time-dependent: While initially only the mass of the beads M = 101 contributes, at later times the fluid within the sphere is dragged as well, such that then the mass is roughly estimated as M ≈ 101 + 4πρR 3 /3 ≈ 214. This gives rise to the initial and final relaxation times τ in ≈ 0.7 and τ f in ≈ 1.5. However, the initial relaxation time cannot be observed, since in the extreme short-time regime ballistic effects play a role. Conversely, the decay with τ f in is clearly visible (see Fig. 2).
After t ≈ 1, the famous long-time tail [24,25] (normalized relaxation function V (t)/V (0) = Bt −3/2 ) sets in. The physical mechanism of this slow relaxation is the fact that the momentum is conserved and hence transported away diffusively from the particle [25]. For this reason, it can only be observed in a sufficiently large box. The prefactor of the power law is known in the colloidal limit, B = (1/12)(Nm/ρ)(πν) −3/2 , where N is the number of beads, and m the bead mass [26,27]. It should be noted that the prefactor of the power-law decay of the unnormalized velocity autocorrelation function does not depend on the properties of the sphere at all, but only on the temperature, and the hydrodynamic properties of the solvent [26,27]. The mass dependence of B results only from the normalization, i. e. the value of V 2 according to the equipartition theorem, V 2 = 3k B T /(Nm). From this consideration it is clear that the short-time value of the sphere's mass (Nm) enters (at t = 0 the shell and the inner fluid are not yet coupled). For our simulation parameters, we find B = 0.342. As Fig. 2 shows, the data exhibit the expected behavior very nicely. Figure 3 illustrates the finite size effect in the diffusive properties by plotting the mean square displacement of the colloidal sphere as a function of time, for two different box sizes L = 6.33R and L = 20R, for Brownian motion at k B T = 1. In the long-time regime, the slope is given by 6D S , where D S is the self-diffusion coefficient. The figure clearly shows that D S increases with box size. This can be explained in terms of hydrodynamic interactions with the periodic images, or, equivalently, in terms of suppression of long-wavelength hydrodynamic modes [20,23]. We have therefore performed a systematic finite-size analysis, and measured D S for various box sizes by integrating the velocity autocorrelation function over time [25]. In Fig. 4 we plot D S /D 0 (D 0 denoting the asymptotic Stokes-Einstein value D 0 = k B T /(6πηR) = 6.9 × 10 −3 ) as a function of 1/L; the expected linear behavior [23] is clearly seen.
We now look at the relaxation of the rotational motion. We performed a similar kick experiment as described above, where now an initial angular velocity ω 0 = 1 was provided to the sphere. The data for the normalized decay function ω(t)/ω(0) are presented in Fig. 5. The sphere dynamics shows a characteristic "raw egg" (damped) rotation pattern with a fast initial decay and a subsequent slower one. For the rotational relaxation of a hard sphere with moment of inertia I and rotational friction coefficient 8πηR 3 , theory predicts a decay according to the Debye law ω(t)/ω(0) = exp (−t/τ r ), where the relaxation time is given by τ r = I/(8πηR 3 ) [28]. Similarly to the effective mass, the effective moment of inertia is expected to increase as a function of time, due to the time-delayed dragging of the fluid. Initially we expect a hollow-sphere value of roughly I in ≈ (2/3)MR 2 ≈ 600, where we used the mass of the outer shell M = 100. A more accurate value is obtained by direct summation over the contributing beads, yielding I in = 546. In the late-time regime the moment of inertia is expected to be the compact-sphere value I f in ≈ (2/5)MR 2 ≈ 770, where the total mass M ≈ 214 has been used. This results in relaxation times τ r,in = 0.34 and τ r,f in = 0.44. These values are roughly consistent with a fit to the data in the interval 0.1 < t < 1.0, which yields a somewhat larger relaxation time τ r = 0.68 (see Fig. 5). Given the general inaccuracy of the estimates, and the crossover of the Debye relaxation function both to short-time ballistic behavior, and long-time hydrodynamic behavior, this deviation is not too surprising.
Similar to the translational motion, the exponential decay is then followed by a power-law long-time tail. Theory predicts ω(t)/ω(0) = (πI/ρ)(4πνt) −5/2 [29]. As we have seen before, the long-time tail in the translational case is governed by the short-time mass as a result of normalization. Similarly, the rotational tail must be controlled by the short-time moment of inertia I in , which also governs the mean square fluctuations of the angular velocity via the equipartition theorem ω 2 = 3k B T /I in . For this reason, we can determine I in by fitting a t −5/2 law to the data, which is much more accurate than our rough geometrical estimate. The fit results in I in = 533, which is reasonably consistent. If we insert this value into the relaxation time expression, we obtain τ r = 0.64, in quite good agreement with the data. It hence seems that fluid dragging effects are not yet very important for the initial Debye relaxation.
At longer times, one can again notice a significant finite-size effect: the curves obtained in the smaller simulation box depart from the asymptotic power-law line earlier, i. e. the long-time diffusion is hindered at the small system sizes. In fact, the power-law regime is inaccessible at L = 10 while for L = 100 it extends up to t = 200. This value characterizes the interval after which the particle starts feeling its own periodic images. The rotational diffusion constant D R is given by the Green-Kubo integral [30] We can evaluate this by again making use of linear response theory: The correlation function ω(t) · ω(0) is identical to the relaxation function presented in Fig. 5, multiplied with the initial value ω 2 , which we know from the equipartition theorem (see above). Using this approach, we have calculated D R for different box sizes L. In an infinite fluid, D R has the Stokes-Einstein-Debye value D R = k B T /(8πηR 3 ) = 0.58×10 −3 , towards which the results indeed converge for L → ∞. Again a 1/L behavior (but weaker than for translational diffusion) is observed, as shown in Fig. 4. It should be noted that the accuracy of the data is slightly hampered by the fact that we had to use a somewhat arbitrary criterion for cutting off the integral function, which does not arrive at a constant value but rather ends up with a linear increase originating from a small nonzero constant angular velocity in the long-time limit (as a manifestation of the finite-size effect).
Finally, we have started to study the effect of charge on the self-diffusion of the colloidal particle. We performed simulations of a sphere with central charge Z C = 10 in an LB box of size L = 40. Ten counterions of charge −1 were also added. The Bjerrum length was set to 2. Technically, the simulation ran without any problems just as well as for the neutral system. In Fig. 6 we compare the decay of the velocity autocorrelation function for the neutral and charged spheres. The difference between the two is not detectable within our error bars at this charge. We cannot say yet much about the influence of the charges on the dynamic properties; we do however expect that for strongly charged systems the Debye layer will effectively increase the hydrodynamic radius and slow down the diffusion. More work needs to be done to resolve this issue.

Summary
In summary, we introduced and tested a new model for simulating colloidal dynamics with inclusion of the hydrodynamic effects on the level of the lattice-Boltzmann equation. The suggested "raspberry"-like colloidal object exhibits the essential features of the diffusion of a spherical colloidal particle. For the first time we combined a model containing explicit charges with accurate treatment of the hydrodynamics. We expect this method to have some advantages over the previously applied scheme [17] which builds upon the model of Ladd [8] and has the problem of charge leakage. In our method: (i) the stick boundary conditions are realized not strictly, i. e. the particle surface is softly bound to the fluid; (ii) the method is easily adaptable to the multiple timestep method where the ionic action on the colloid can be averaged out somewhat within the integration step for the solvent. We should note that such a hybrid simulation requires still some computational effort -two million hybrid steps for the colloid consisting of 100 beads in a box of 20 × 20 × 20 LB cells take about 12 hours on a 2 GHz Pentium 4 processor, while it is expected to run a few times faster on more sophisticated hardware. Here, one step is a MD step of the total system of beads, plus one LB update of the whole lattice, where the latter part completely dominates the CPU effort. Some improvement is possible by not updating the LB fluid at every MD step, as it was done here. For dense systems, where the MD part is no longer negligible, one should also optimize the interaction potentials. Finally, we would like to mention that the algorithm can be fully parallelized and we are working on the parallel version.