Multi-Particle Collision Dynamics -- a Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids
MMulti-Particle Collision Dynamics — aParticle-Based Mesoscale Simulation Approach tothe Hydrodynamics of Complex Fluids
G. Gompper , T. Ihle , D.M. Kroll , and R.G. Winkler Theoretical Soft Matter and Biophysics, Institut f¨ur Festk¨orperforschung,Forschungszentrum J¨ulich, 52425 J¨ulich, Germany Department of Physics, North Dakota State University, Fargo, North Dakota,58105-5566In this review, we describe and analyze a mesoscale simulation method for fluidflow, which was introduced by Malevanets and Kapral in 1999, and is now calledmulti-particle collision dynamics (MPC) or stochastic rotation dynamics (SRD). Themethod consists of alternating streaming and collision steps in an ensemble of pointparticles. The multi-particle collisions are performed by grouping particles in colli-sion cells, and mass, momentum, and energy are locally conserved. This simulationtechnique captures both full hydrodynamic interactions and thermal fluctuations.The first part of the review begins with a description of several widely used MPCalgorithms and then discusses important features of the original SRD algorithm andfrequently used variations. Two complementary approaches for deriving the hydro-dynamic equations and evaluating the transport coefficients are reviewed. It is thenshown how MPC algorithms can be generalized to model non-ideal fluids, and binarymixtures with a consolute point. The importance of angular-momentum conserva-tion for systems like phase-separated liquids with different viscosities is discussed.The second part of the review describes a number of recent applications of MPCalgorithms to study colloid and polymer dynamics, the behavior of vesicles and cellsin hydrodynamic flows, and the dynamics of viscoelastic fluids.
PACS number(s): 47.11.-j, 05.40.-a, 02.70.Ns “Soft Matter” is a relatively new field of research that encompasses tradi-tional complex fluids such as amphiphilic mixtures, colloidal suspensions, andpolymer solutions, as well as a wide range of phenomena including chemicallyreactive flows (combustion), the fluid dynamics of self-propelled objects, andthe visco-elastic behavior of networks in cells. One characteristic feature ofall these systems is that phenomena of interest typically occur on mesoscopiclength-scales—ranging from nano- to micrometers—and at energy scales com-parable to the thermal energy k B T . a r X i v : . [ c ond - m a t . s o f t ] A ug G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler
Because of the complexity of these systems, simulations have played a par-ticularly important role in soft matter research. These systems are challengingfor conventional simulation techniques due to the presence of disparate time,length, and energy scales. Biological systems present additional challenges be-cause they are often far from equilibrium and are driven by strong spatiallyand temporally varying forces. The modeling of these systems often requiresthe use of “coarse-grained” or mesoscopic approaches that mimic the behaviorof atomistic systems on the length scales of interest. The goal is to incorporatethe essential features of the microscopic physics in models which are compu-tationally efficient and are easily implemented in complex geometries and onparallel computers, and can be used to predict behavior, test physical the-ories, and provide feedback for the design and analysis of experiments andindustrial applications.In many situations, a simple continuum description based on the Navier-Stokes equation is not sufficient, since molecular-level details—including ther-mal fluctuations—play a central role in determining the dynamic behavior.A key issue is to resolve the interplay between thermal fluctuations, hydro-dynamic interactions, and spatio-temporally varying forces. One well-knownexample of such systems are microemulsions—a dynamic bicontinuous net-work of intertwined mesoscopic patches of oil and water—where thermal fluc-tuations play a central role in creating this phase. Other examples includeflexible polymers in solution, where the coil state and stretching elasticity aredue to the large configurational entropy. On the other hand, atomistic molec-ular dynamics simulations retain too many microscopic degrees of freedom,consequently requiring very small time steps in order to resolve the high fre-quency modes. This makes it impossible to study long timescale behavior suchas self-assembly and other mesoscale phenomena.In order to overcome these difficulties, considerable effort has been devotedto the development of mesoscale simulation methods such as Dissipative Par-ticle Dynamics [1–3], Lattice-Boltzmann [4–6], and Direct Simulation MonteCarlo [7–9]. The common approach of all these methods is to “average out”irrelevant microscopic details in order to achieve high computational efficiencywhile keeping the essential features of the microscopic physics on the lengthscales of interest. Applying these ideas to suspensions leads to a simplified,coarse-grained description of the solvent degrees of freedom, in which embed-ded macromolecules such as polymers are treated by conventional moleculardynamics simulations.All these approaches are essentially alternative ways of solving the Navier-Stokes equation and its generalizations. This is because the hydrodynamicequations are expressions for the local conservation laws of mass, momentum,and energy, complemented by constitutive relations which reflect some aspectsof the microscopic details. Frisch et al. [10] demonstrated that discrete algo-rithms can be constructed which recover the Navier-Stokes equation in thecontinuum limit as long as these conservation laws are obeyed and space isdiscretized in a sufficiently symmetric manner. ulti-Particle Collision Dynamics The first model of this type was a cellular automaton, called the Lattice-Gas-Automaton (LG). The algorithm consists of particles which jump betweennodes of a regular lattice at discrete time intervals. Collisions occur whenmore than one particle jumps to the same node, and collision rules are cho-sen which impose mass and momentum conservation. The Lattice-Boltzmannmethod (LB)—which follows the evolution of the single-particle probabilitydistribution at each node—was a natural generalization of this approach. LBsolves the Boltzmann equation on a lattice with a small set of discrete veloci-ties determined by the lattice structure. The price for obtaining this efficiencyis numerical instabilities in certain parameter ranges. Furthermore, as origi-nally formulated, LB did not contain any thermal fluctuations. It became clearonly very recently (and only for simple liquids) how to restore fluctuations byintroducing additional noise terms to the algorithm [11].Except for conservation laws and symmetry requirements, there are rel-atively few constraints on the structure of mesoscale algorithms. However,the constitutive relations and the transport coefficients depend on the detailsof the algorithm, so that the temperature and density dependencies of thetransport coefficients can be quite different from those of real gases or liq-uids. However, this is not a problem as long as the functional form of theresulting hydrodynamic equations is correct. The mapping to real systems isachieved by tuning the relevant characteristic numbers, such as the Reynoldsand Peclet numbers [12, 13], to those of a given experiment. When it is notpossible to match all characteristic numbers, one concentrates on those whichare of order unity, since this indicates that there is a delicate balance betweentwo effects which need to be reproduced by the simulation. On occasion, thiscan be difficult, since changing one internal parameter, such as the mean freepath, usually affects all transport coefficients in different ways, and it mayhappen that a given mesoscale algorithm is not at all suited for a given ap-plication [14–17].In this review we focus on the development and application of a particle-based mesoscopic simulation technique which was recently introduced by Mal-evanets and Kapral [18, 19]. The algorithm, which consists of discrete stream-ing and collision steps, shares many features with Bird’s Direct SimulationMonte Carlo (DSMC) approach [7]. Collisions occur at fixed discrete time in-tervals, and although space is discretized into cells to define the multi-particlecollision environment, both the spatial coordinates and the velocities of theparticles are continuous variables. Because of this, the algorithm exhibits un-conditional numerical stability and has an H -theorem [18, 20]. In this review,we will use the name multi-particle collision dynamics (MPC) to refer to thisclass of algorithms. In the original and most widely used version of MPC,collisions consist of a stochastic rotation of the relative velocities of the par-ticles in a collision cell. We will refer to this algorithm as stochastic rotationdynamics (SRD) in the following.One important feature of MPC algorithms is that the dynamics is well-defined for an arbitrary time step, ∆t . In contrast to methods such as molecu- G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler lar dynamics simulations (MD) or dissipative particle dynamics (DPD), whichapproximate the continuous-time dynamics of a system, the time step doesnot have to be small. MPC defines a discrete-time dynamics which has beenshown to yield the correct long-time hydrodynamics; one consequence of thediscrete dynamics is that the transport coefficients depend explicitly on ∆t . Infact, this freedom can be used to tune the Schmidt number, Sc [15]; keepingall other parameters fixed, decreasing ∆t leads to an increase in Sc . For smalltime steps, Sc is larger than unity (as in a dense fluid), while for large timesteps, Sc is of order unity, as in a gas.Because of its simplicity, SRD can be considered an “Ising model” for hy-drodynamics, since it is Galilean invariant (when a random grid shift of thecollision cells is performed before each collision step [21]) and incorporates allthe essential dynamical properties in an algorithm which is remarkably easyto analyze. In addition to the conservation of momentum and mass, SRD alsolocally conserves energy, which enables simulations in the microcanonical en-semble. It also fully incorporates both thermal fluctuations and hydrodynamicinteractions. Other more established methods, such as Brownian Dynamics(BD) can also be augmented to include hydrodynamic interactions. However,the additional computational costs are often prohibitive [22, 23]. In addition,hydrodynamic interactions can be easily switched off in MPC algorithms,making it easy to study the importance of hydrodynamic interactions [24,25].It must, however, be emphasized that all local algorithms such as MPC,DPD, and LB model compressible fluids, so that it takes time for the hydro-dynamic interactions to “propagate” over longer distances. As a consequence,these methods become quite inefficient in the Stokes limit, where the Reynoldsnumber approaches zero. Algorithms which incorporate an Oseen tensor donot share this shortcoming.The simplicity of the SRD algorithm has made it possible to derive ana-lytic expressions for the transport coefficients which are valid for both largeand small mean free paths [26–28]. This is usually very difficult to do for othermesoscale particle-based algorithms. Take DPD as an example: the viscositymeasured in Ref. [29] is about 50% smaller than the value predicted theo-retically in the same paper. For SRD, the agreement is generally better than1%.MPC is particularly well suited for (i) studying phenomena where boththermal fluctuations and hydrodynamics are important, (ii) for systems withReynolds and Peclet numbers of order 0 . ulti-Particle Collision Dynamics are not crucial, but fluctuations are, one might be better served using Langevinor Brownian Dynamics.This review consists of two parts. The first part begins with a descriptionof several widely used MPC algorithms in Sec. 2, and then discusses impor-tant features of the original SRD algorithm and a frequently used variation(MPC-AT), which effectively thermostats the system by replacing the relativevelocities of particles in a collision cell with newly generated Gaussian randomnumbers in the collision step. After a qualitative discussion of the static anddynamic properties of MPC fluids in Sec. 3, two alternative approaches forderiving the hydrodynamic equations and evaluating the transport coefficientsare described. First, in Sec. 4, discrete-time projection operator methods arediscussed and the explicit form of the resulting Green-Kubo relations for thetransport coefficients are given and evaluated. Subsequently, in Sec. 5, analternative non-equilibrium approach is described. The two approaches com-plement each other, and the predictions of both methods are shown to be incomplete agreement. It is then shown in Sec. 6 how MPC algorithms can begeneralized to model non-ideal fluids and binary mixtures, Finally, variousapproaches for implementing slip and no-slip boundary conditions—as well asthe coupling of embedded objects to a MPC solvent—are described in Sec. 7.In Sec. 8, the importance of angular-momentum conservation is discussed,in particular in systems of phase-separated fluids with different viscositiesunder flow. An important aspect of mesoscale simulations is the possibilityto directly assert the effect of hydrodynamic interactions by switching themoff, while retaining the same thermal fluctuations and similar friction coeffi-cients; in MPC, this can be done very efficiently by an algorithm describedin Sec. 9. The second part of the review describes a number of recent appli-cations of MPC algorithms to study colloid and polymer dynamics, and thebehavior of vesicles and cells in hydrodynamic flows. Sec. 10 focuses on thenon-equilibrium behavior of colloidal suspensions, the dynamics of dilute so-lutions of linear polymers both in equilibrium and under flow conditions, andthe properties of star polymers—also called ultra-soft colloids—in shear flow.Sec. 11 is devoted to the review of recent simulation results for membranes inflow. After a short introduction to the modeling of membranes with differentlevels of coarse-graining, the behavior of fluid vesicles and red blood cells,both in shear and capillary flow, is discussed. Finally, a simple extension ofMPC for viscoelastic solvents is described in Sec. 12, where the point particlesof MPC for Newtonian fluids are replaced by harmonic dumbbells.A discussion of several complementary applications—such as chemicallyreactive flows and self-propelled objects—can be found in a recent review ofMPC by R. Kapral [30]. G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler
In the following, we use the term multi-particle collision dynamics (MPC) todescribe the generic class of particle-based algorithms for fluid flow which con-sist of successive free-streaming and multi-particle collision steps. The namestochastic rotation dynamics (SRD) is reserved for the most widely used algo-rithm which was introduced by Malevanets and Kapral [18]. The name refersto the fact that the collisions consist of a random rotation of the relative ve-locities δ v i = v i − u of the particles in a collision cell, where u is the meanvelocity of all particles in a cell. There are a number of other MPC algorithmswith different collision rules [31–33]. For example, one class of algorithms usesmodified collision rules which provide a nontrivial “collisional” contributionto the equation of state [33,34]. As a result, these models can be used to modelnon-ideal fluids or multi-component mixtures with a consolute point. In SRD, the solvent is modeled by a large number N of point-like particlesof mass m which move in continuous space with a continuous distribution ofvelocities. The algorithm consists of individual streaming and collision steps.In the streaming step, the coordinates, r i ( t ), of all solvent particles at time t are simultaneously updated according to r i ( t + ∆t ) = r i ( t ) + ∆t v i ( t ) , (1)where v i ( t ) is the velocity of particle i at time t and ∆t is the value of thediscretized time step.In order to define the collisions, particles are sorted into cells, and theyinteract only with members of their own cell. Typically, the system is coarse-grained into cells of a regular, typically cubic, grid with lattice constant a .In practice, lengths are often measured in units of a , which corresponds tosetting a = 1. The average number of particles per cell, M , is typically chosento be between three and 20. The actual number of particles in cell at a giventime, which fluctuates, will be denoted by N c . The collision step consists of arandom rotation R of the relative velocities δ v i = v i − u of all the particlesin the collision cell, v i ( t + ∆t ) = u ( t ) + R · δ v i ( t ) . (2)All particles in the cell are subject to the same rotation, but the rotationsin different cells and at different times are statistically independent. Thereis a great deal of freedom in how the rotation step is implemented, and anystochastic rotation matrix which satisfies semi-detailed balance can be used.Here, we describe the most commonly used algorithm. In two dimensions, R is a rotation by an angle ± α , with probability 1 /
2. In three dimensions, a ulti-Particle Collision Dynamics rotation by a fixed angle α about a randomly chosen axis is typically used.Note that rotations by an angle − α need not be considered, since this amountsto a rotation by an angle α about an axis with the opposite orientation. If wedenote the randomly chosen rotation axis by ˆ R , the explicit collision rule inthree dimensions is v i ( t + ∆t ) = u ( t ) + δ v i, ⊥ ( t ) cos( α )+ ( δ v i, ⊥ ( t ) × ˆ R ) sin( α ) + δ v i, (cid:107) ( t ) , (3)where ⊥ and (cid:107) are the components of the vector which are perpendicularand parallel to the random axis ˆ R , respectively. Malevanets and Kapral [18]have shown that there is an H -theorem for the algorithm, that the equilib-rium distribution of velocities is Maxwellian, and that it yields the correcthydrodynamic equations with an ideal-gas equation of state.In its original form [18,19], the SRD algorithm was not Galilean invariant.This is most pronounced at low temperatures or small time steps, where themean free path, λ = ∆t (cid:112) k B T /m , is smaller than the cell size a . If the particlestravel a distance between collisions which is small compared to the cell size,essentially the same particles collide repeatedly before other particles enterthe cell or some of the participating particles leave the cell. For small λ , largenumbers of particles in a given cell remain correlated over several time steps.This leads to a breakdown of the molecular chaos assumption—i.e., particlesbecome correlated and retain information of previous encounters. Since thesecorrelations are changed by a homogeneous imposed flow field, V , Galileaninvariance is destroyed, and the transport coefficients depend on both themagnitude and direction of V .Ihle and Kroll [20, 21] showed that Galilean invariance can be restoredby performing a random shift of the entire computational grid before everycollision step. The grid shift constantly groups particles into new collisionneighborhoods; the collision environment no longer depends on the magni-tude of an imposed homogeneous flow field, and the resulting hydrodynamicequations are Galilean invariant for arbitrary temperatures and Mach number.This procedure is implemented by shifting all particles by the same randomvector with components uniformly distributed in the interval [ − a/ , a/
2] be-fore the collision step. Particles are then shifted back to their original positionsafter the collision.In addition to restoring Galilean invariance, this grid-shift procedure accel-erates momentum transfer between cells and leads to a collisional contributionto the transport coefficients. If the mean free path λ is larger than a/
2, theviolation of Galilean invariance without grid shift is negligible, and it is notnecessary to use this procedure.
G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler
As noted by Pooley and Yeomans [35] and confirmed in Ref. [28], the macro-scopic stress tensor of SRD is not symmetric in ∂ α v β . The reason for thisis that the multi-particle collisions do not, in general, conserve angular mo-mentum. The problem is particularly pronounced for small mean free paths,where asymmetric collisional contributions to the stress tensor dominate theviscosity (see Sec. 4.1.1). In contrast, for mean free paths larger than the cellsize, where kinetic contributions dominate, the effect is negligible.An anisotropic stress tensor means that there is non-zero dissipation ifthe entire fluid undergoes a rigid-body rotation, which is clearly unphysical.However, as emphasized in Ref. [28], this asymmetry is not a problem formost applications in the incompressible (or small Mach number) limit, sincethe form of the Navier-Stokes equation is not changed. This is in accordancewith results obtained in SRD simulations of vortex shedding behind an obsta-cle [36], and vesicle [37] and polymer dynamics [14]. In particular, it has beenshown that the linearized hydrodynamic modes are completely unaffected intwo dimensions; in three dimensions only the sound damping is slightly mod-ified [28].However, very recently G¨otze et al. [38] identified several situations involv-ing rotating flow fields in which this asymmetry leads to significant deviationsfrom the behavior of a Newtonian fluid. This includes (i) systems in whichboundary conditions are defined by torques rather than prescribed velocities,(ii) mixtures of liquids with a viscosity contrast, and (iii) polymers with alocally high monomer density and a monomer-monomer distance on the orderof or smaller than the lattice constant, a , embedded in a MPC fluid. A moredetailed discussion will be presented in Sec. 8 below.For the SRD algorithm, it is possible to restore angular momentum con-servation by having the collision angle depend on the specific positions of theparticles within a collision cell. Such a modification was first suggested by Ry-der [39] for SRD in two dimensions. She showed that the angular momentumof the particles in a collision cell is conserved if the collision angle α is chosensuch thatsin( α ) = − AB/ ( A + B ) and cos( α ) = ( A − B ) / ( A + B ) (4)where A = N c (cid:88) [ r i × ( v i − u )] | z and B = N c (cid:88) r i · ( v i − u ) . (5)When the collision angles are determined in this way, the viscous stress tensoris symmetric. Note, however, that evaluating Eq. (4) is time-consuming, sincethe collision angle needs to be computed for every collision cell every timestep. This typically increases the CPU time by a factor close to two. ulti-Particle Collision Dynamics A general procedure for implementing angular-momentum conservation inmulti-particle collision algorithms was introduced by Noguchi et al. [32]; it isdiscussed in the following section.
A stochastic rotation of the particle velocities relative to the center-of-massvelocity is not the only possibility for performing multi-particle collisions.In particular, MPC simulations can be performed directly in the canonicalensemble by employing an Anderson thermostat (AT) [31, 32]; the resultingalgorithm will be referred to as MPC-AT-a. In this algorithm, instead of per-forming a rotation of the relative velocities, { δ v i } , in the collision step, newrelative velocities are generated. The components of { δ v rani } are Gaussianrandom numbers with variance (cid:112) k B T /m . The collision rule is [32, 38] v i ( t + ∆t ) = u ( t ) + δ v rani = u ( t ) + v rani − (cid:88) j ∈ cell v ranj /N c , (6)where N c is the number of particles in the collision cell, and the sum runs overall particles in the cell. It is important to note that MPC-AT is both a colli-sion procedure and a thermostat. Simulations are performed in the canonicalensemble, and no additional velocity rescaling is required in non-equilibriumsimulations, where there is viscous heating.Just as SRD, this algorithm conserves momentum at the cell level butnot angular momentum. Angular momentum conservation can be restored[32, 39] by imposing constraints on the new relative velocities. This leads toan angular-momentum conserving modification of MPC-AT [32, 38], denotedMPC-AT+ a . The collision rule in this case is v i ( t + ∆t ) = u ( t ) + v i,ran − (cid:88) cell v i,ran /N c + m Π − (cid:88) j ∈ cell (cid:2) r j,c × ( v j − v ranj ) (cid:3) × r i,c , (7)where Π is the moment of inertia tensor of the particles in the cell, and r i,c = r i − R c is the relative position of particle i in the cell and R c is thecenter of mass of all particles in the cell.When implementing this algorithm, an unbiased multi-particle collision isfirst performed, which typically leads to a small change of angular momen-tum, ∆ L . By solving the linear equation − ∆ L = Π · ω , the angular velocity ω which is needed to cancel the initial change of angular momentum is then de-termined. The last term in Eq. (7) restores this angular momentum deficiency.MPC-AT can be adapted for simulations in the micro-canonical ensemble byimposing an additional constraint on the values of the new random relativevelocities [32]. Because d Gaussian random numbers per particle are required at every it-eration, where d is the spatial dimension, the speed of the random numbergenerator is the limiting factor for MPC-AT. In contrast, the efficiency ofSRD is rather insensitive to the speed of the random number generator sinceonly d − decrease when the number of particlesper cell is increased, while they increase for SRD. A longer relaxation timemeans that a larger number of time steps is required for transport coefficientsto reach their asymptotic value. This could be of importance when fast oscil-latory or transient processes are investigated. As a consequence, when usingSRD, the average number of particles per cell should be in range of 3 − M ) − , where M is the average number ofparticles in a collision cell. The MPC-AT algorithm discussed in Sec. 2.2 provides a very efficient particle-level thermostating of the system. However, it is considerably slower thanthe original SRD algorithm, and there are situations in which the additionalfreedom offered by the choice of SRD collision angle can be useful.Thermostating is required in any non-equilibrium MPC simulation, wherethere is viscous heating. A basic requirement of any thermostat is that itdoes not violate local momentum conservation, smear out local flow profiles,or distort the velocity distribution too much. When there is homogeneousheating, the simplest way to maintain a constant temperature is to just rescalevelocity components by a scale factor S , v newα = Sv α , which adjusts the totalkinetic energy to the desired value. This can be done with just a single globalscale factor, or a local factor which is different in every cell. For a knownmacroscopic flow profile, u , like in shear flow, the relative velocities v − u canbe rescaled. This is known as a profile-unbiased thermostat; however, it hasbeen shown to have deficiencies in molecular dynamics simulations [41].Here we describe an alternative thermostat which exactly conserves mo-mentum in every cell and is easily incorporated into the MPC collision step. Itwas originally developed by Heyes for constant-temperature molecular dynam-ics simulations; however, the original algorithm described in Ref. [42] violates ulti-Particle Collision Dynamics detailed balance. The thermostat consists of the following procedure which isperformed independently in every collision cell as part of the collision step.1. Randomly select a real number ψ ∈ [1 , c ], where c is a small numberbetween 0.05 and 0.3 which determines the strength of the thermostat.2. Accept this number as a scaling factor S = ψ with probability 1 /
2; oth-erwise, take S = 1 /ψ .3. Create another random number ξ ∈ [0 , ξ issmaller than the acceptance probability p A = min(1 , A ), where A = S d ( N c − exp (cid:32) − m k B T N c (cid:88) i =1 ( v i − u ) { S − } (cid:33) . (8) d is the spatial dimension, and N c is the number of particles in the cell.The prefactor in Eq. (8) is an entropic contribution which accounts forthe fact that the phase-space volume changes if the velocities are rescaled.4. If the attempt is accepted, perform a stochastic rotation with the scaledrotation matrix S R . Otherwise, use the rotation matrix R .This thermostat reproduces the Maxwell velocity distribution and does notchange the viscosity of the fluid. It gives excellent equilibration, and the devi-ation of the measured kinetic temperature from T is smaller than 0 . c controls the rate at which the kinetic temperature relaxes to T ,and in agreement with experience from MC-simulations, an acceptance rate inthe range of 50% to 65% leads to the fastest relaxation. For these acceptancerates, the relaxation time is of order 5 −
10 time steps. The correspondingvalue for c depends on the particle number N c ; in two dimensions, it is about0 . N c = 7 and decreases to 0 .
05 for N c = 100. This thermostat has beensuccessfully applied to SRD simulations of sedimenting charged colloids [16]. The previous section outlines several multi-particle algorithms. A detaileddiscussion of the link between the microscopic dynamics described by Eqs. (1)and (2) or (3) and the macroscopic hydrodynamic equations, which describethe behavior at large length and time scales, requires a more careful analysisof the corresponding Liouville operator L . Before describing this approach inmore detail, we provide a more heuristic discussion of the equation of stateand of one of the transport coefficients, the shear viscosity, using more familiarapproaches for analyzing the behavior of dynamical systems. In a homogeneous fluid, the pressure is the normal force exerted by the fluidon one side of a unit area on the fluid on the other side; expressed somewhat differently, it is the momentum transfer per unit area per unit time across animaginary (flat) fixed surface. There are both kinetic and virial contributionsto the pressure. The first arises from the momentum transported across thesurface by particles that cross the surface in the unit time interval; it yields theideal-gas contribution, P id = N k B T /V , to the pressure. For classical particlesinteracting via pair-additive, central forces, the intermolecular “potential”contribution to the pressure can be determined using the method introducedby Irving and Kirkwood [43]. A clear discussion of this approach is given byDavis in Ref. [44], where it is shown to lead to the virial equation of state ofa homogeneous fluid, P = N k B TV + 13 V (cid:88) i (cid:104) r i · F i (cid:105) , (9)in three dimensions, where F i is the force on particle i due to all the otherparticles, and the sum runs over all particles of the system.The kinetic contribution to the pressure, P id = N k B T /V , is clearly presentin all MPC algorithms. For SRD, this is the only contribution. The reason isthat the stochastic rotations, which define the collisions, transport (on aver-age) no net momentum across a fixed dividing surface. More general MPCalgorithms (such as those discussed in Sec. 6) have an additional contributionto the virial equation of state. However, instead of an explicit force F i asin Eq. (9), the contribution from the multi-particle collisions is a force of theform m∆v i /∆t , and the role of the particle position, r i , is played by a variablewhich denotes the cell-partners which participate in the collision [33, 45]. Just as for the pressure, there are both kinetic and collisional contributions tothe transport coefficients. We present here a heuristic discussion of these con-tributions to the shear viscosity, since it illustrates rather clearly the essentialphysics and provides background for subsequent technical discussions.Consider a reference plane (a line in two dimensions) with normal in the y -direction embedded in a homogeneous fluid in equilibrium. The fluid belowthe plane exerts a mean force p y per unit area on the fluid above the plane;by Newton’s third law, the fluid above the plane must exert a mean force − p y on the fluid below the plane. The normal force per unit area is justthe mean pressure, P , so that p yy = P . In a homogeneous simple fluid inwhich there are no velocity gradients, there is no tangential force, so that, forexample, p yx = 0. p αβ is called the pressure tensor , and the last result is justa statement of the well known fact that the pressure tensor in a homogeneoussimple fluid at equilibrium has no off-diagonal elements; the diagonal elementsare all equal to the mean pressure P .Consider a shear flow with a shear rate ˙ γ = ∂u x ( y ) /∂y . In this case, thereis a tangential stress on the reference surface because of the velocity gradient ulti-Particle Collision Dynamics normal to the plane. In the small gradient limit, the dynamic viscosity , η , isdefined as the coefficient of proportionality between the tangential stress, p yx ,and the normal gradient of the imposed velocity gradient, p yx = − η ˙ γ. (10)The kinematic viscosity , ν , is related to η by ν = η/ρ , where ρ = nm is themass density, with n the number density of the fluid and m the particle mass. Kinetic contribution to the shear viscosity:
The kinetic contribution to theshear viscosity comes from transverse momentum transport by the flow offluid particles. This is the dominant contribution to the viscosity of gases.The following analogy may make this origin of viscosity clearer. Consider twoships moving side by side in parallel, but with different speeds. If the sailorson the two ships constantly throw sand bags from their ship onto the other,there will be a transfer of momentum between to two ships so that the slowership accelerates and the faster ship decelerates. This can be interpreted as aneffective friction, or kinetic viscosity, between the ships. There are no directforces between the ships, and the transverse momentum transfer originatessolely from throwing sandbags from one ship to the other.A standard result from kinetic theory is that the kinetic contribution tothe shear viscosity in simple gases is [46] η kin ∼ nm ¯ vλ, (11)where λ is the mean free path and ¯ v is the thermal velocity. Using the factthat λ ∼ ¯ v∆t for SRD and that ¯ v ∼ (cid:112) k B T /m , relation (11) implies that η kin ∼ nk B T ∆t, or equivalently , ν kin ∼ k B T ∆t/m, (12)which is, as more detailed calculations presented later will show, the correctdependence on n , k B T , and ∆t . In fact, the general form for the kineticcontribution to the kinematic viscosity is ν kin = k B T ∆tm f kin ( d, M, α ) , (13)where d is the spatial dimension, M is the mean number of particles percell, and α is the SRD collision angle. Another way of obtaining this resultis to use the analogy with a random walk: The kinematic viscosity is thediffusion coefficient for momentum diffusion. At large mean free path, λ/a (cid:29)
1, momentum is primarily transported by particle translation (as in the shipanalogy). The mean distance a particle streams during one time step, ∆t ,is λ . According to the theory of random walks, the corresponding diffusioncoefficient scales as ν kin ∼ λ /∆t ∼ k B T ∆t/m .Note that in contrast to a “real” gas, for which the viscosity has a squareroot dependence on the temperature, ν kin ∼ T for SRD. This is because themean free path of a particle in SRD does not depend on density; SRD allowsparticles to stream right through each other between collisions. Note, however,that SRD can be easily modified to give whatever temperature dependence is desired. For example, an additional temperature-dependent collision prob-ability can be introduced; this would be of interest, e.g., for a simulation ofrealistic shock-wave profiles. Collisional contribution to the shear viscosity:
At small mean free paths, λ/a (cid:28)
1, particles “stream” only a short distance between collisions, and themulti-particle “collisions” are the primary mechanism for momentum trans-port. These collisions redistribute momenta within cells of linear size a . Thismeans that momentum “hops” an average distance a in one time step, leadingto a momentum diffusion coefficient ν col ∼ a /∆t . The general form of thecollisional contribution to the shear viscosity is therefore ν col = a ∆t f col ( d, M, α ) . (14)This is indeed the scaling observed in numerical simulations at small meanfree path.The kinetic contribution dominates for λ (cid:29) a , while the collisional con-tribution dominates in the opposite limit. Two other transport coefficients ofinterest are the thermal diffusivity, D T , and the single particle diffusion co-efficient, D . Both have the dimension m /sec. As dimensional analysis wouldsuggest, the kinetic and collisional contributions to D T exhibit the same char-acteristic dependencies on λ , a , and ∆t described by Eqs. (13) and (14). Sincethere is no collisional contribution to the diffusion coefficient, D ∼ λ /∆t .Two complementary approaches have been used to derive the transport co-efficients of the SRD fluid. The first is an equilibrium approach which utilizesa discrete projection operator formalism to obtain Green-Kubo (GK) relationswhich express the transport coefficients as sums over the autocorrelation func-tions of reduced fluxes. This approach was first utilized by Malevanets andKapral [19], and later extended by Ihle, Kroll and T¨uzel [20, 27, 28] to in-clude collisional contributions and arbitrary rotation angles. This approach isdescribed in Sec. 4.1.The other approach uses kinetic theory to calculate the transport coef-ficients in stationary non-equilibrium situation such as shear flow. The firstapplication of this approach to SRD was presented in Ref. [21], where the col-lisional contribution to the shear viscosity for large M , where particle numberfluctuations can be ignored, was calculated. This scheme was later extendedby Kikuchi et al. [26] to include fluctuations in the number of particles percell, and then used to obtain expressions for the kinetic contributions to shearviscosity and thermal conductivity [35]. This non-equilibrium approach is de-scribed in Sec. 5. A projection operator formalism for deriving the linearized hydrodynamicequations and Green-Kubo (GK) relations for the transport coefficients of ulti-Particle Collision Dynamics molecular fluids was originally introduced by Zwanzig [47–49] and lateradapted for lattice gases by Dufty and Ernst [50]. With the help of this formal-ism, explicit expressions for both the reversible (Euler) as well as dissipativeterms of the long-time, large-length-scale hydrodynamics equations for thecoarse-grained hydrodynamic variables were derived. In addition, the result-ing GK relations enable explicit calculations of the transport coefficients ofthe fluid. This work is summarized in Sec. 4.1. An analysis of the equilibriumfluctuations of the hydrodynamic modes can then be used to directly measurethe shear and bulk viscosities as well as the thermal diffusivity. This approachis described in Sec. 4.2, where SRD results for the dynamic structure factorare discussed. The Green-Kubo (GK) relations for SRD differ from the well-known continu-ous versions due to the discrete-time dynamics, the underlying lattice struc-ture, and the multi-particle interactions. In the following, we briefly outlinethis approach for determining the transport coefficients. More details can befound in Refs. [20, 27].The starting point of this theory are microscopic definitions of local hy-drodynamic densities A β . These “slow” variables are the local number, mo-mentum, and energy density. At the cell level, they are defined as A β ( ξ ) = N (cid:88) i =1 a β,i d (cid:89) γ =1 Θ (cid:16) a − (cid:12)(cid:12)(cid:12) ξ γ + a − r iγ (cid:12)(cid:12)(cid:12)(cid:17) , (15)with the discrete cell coordinates ξ = a m , where m β = 1 , . . . , L , for eachspatial component. a ,i = 1 is the particle density, { a β,i } = m { v i ( β − } , with β = 2 , ..., d + 1, are the components of the particle momenta, and a d +2 ,i = mv i / i . d is the spatial dimension, and r i and v i are position and velocity of particle i , respectively. A β ( ξ ), for β = 2 , . . . , d + 2, are cell level coarse-grained densities. Forexample, A ( ξ ) is the x -component of the total momentum of all the particlesin cell ξ at the given time. Note that the particle density, A , was not coarse-grained in Ref. [20], i.e., the Θ functions in Eq. (15) were replaced by a δ -function. This was motivated by the fact that during collisions the particlenumber is trivially conserved in areas of arbitrary size, whereas energy andmomentum are only conserved at the cell level.The equilibrium correlation functions for the conserved variables are de-fined by (cid:104) δA β ( r , t ) δA γ ( r (cid:48) , t (cid:48) ) (cid:105) , where (cid:104) δA (cid:105) = A − (cid:104) A (cid:105) , and the brackets denotean average over the equilibrium distribution. In a stationary, translationallyinvariant system, the correlation functions depend only on the differences r − r (cid:48) and t − t (cid:48) , and the Fourier transform of the matrix of correlation functions is G αβ ( k , t ) = 1 V (cid:104) δA ∗ β ( k , δA γ ( k , t ) (cid:105) , (16) where the asterisk denotes the complex conjugate, and the spatial Fouriertransforms of the densities are given by A β ( k ) = (cid:88) j a β,j e i k · ξ j , (17)where ξ j is the coordinate of the cell occupied by particle j . k = 2 π n / ( aL ) isthe wave vector, where n β = 0 , ± , . . . , ± ( L − , L for the spatial components.To simplify notation, we omit the wave-vector dependence of G αβ in thissection.The collision invariants for the conserved densities are (cid:88) j e i k · ξ sj ( t + ∆t ) [ a β,j ( t + ∆t ) − a β,j ( t )] = 0 , (18)where ξ sj is the coordinate of the cell occupied by particle j in the shifted system. Starting from these conservation laws, a projection operator can beconstructed that projects the full SRD dynamics onto the conserved fields [20].The central result is that the discrete Laplace transform of the linearizedhydrodynamic equations can be written as[ s + ikΩ + k Λ ] G ( k , s ) = 1 ∆t G (0) R ( k ) , (19)where R ( k ) = [1 + ∆t ( ikΩ + k Λ ] − is the residue of the hydrodynamicpole [20]. The linearized hydrodynamic equations describe the long-time large-length-scale dynamics of the system, and are valid in the limits of small k and s . The frequency matrix Ω contains the reversible (Euler) terms of thehydrodynamic equations. Λ is the matrix of transport coefficients. The discreteGreen-Kubo relation for the matrix of viscous transport coefficients is [20] Λ αβ (ˆ k ) ≡ ∆tN k B T ∞ (cid:88) t =0 (cid:48) (cid:104) ˆ k λ σ αλ (0) | ˆ k λ (cid:48) σ βλ (cid:48) ( t ) (cid:105) , (20)where the prime on the sum indicates that the t = 0 term has the relativeweight 1/2. σ αβ = P δ αβ − p αβ is the viscous stress tensor. The reduced fluxesin Eq. (20) are given byˆ k λ σ αλ ( t ) = m∆t (cid:88) j (cid:18) − v jα ( t ) ˆk · [ ∆ ξ j ( t ) + ∆v jα ( t ) ∆ ξ sj ( t )] + ∆td ˆ k α v j ( t ) (cid:19) (21)for α = 1 , . . . , d , with ∆ ξ j ( t ) = ξ j ( t + ∆t ) − ξ j ( t ), ∆ ξ sj ( t + ∆t ) = ξ j ( t + ∆t ) − ξ sj ( t + ∆t ), and ∆v xj ( t ) = v xj ( t + ∆t ) − v xj ( t ). ξ j ( t ) is the cell coordinateof particle j at time t , while ξ sj is it’s cell coordinate in the (stochastically)shifted frame. The corresponding expressions for the thermal diffusivity andself-diffusion coefficient can be found in Ref. [20].The straightforward evaluation of the GK relations for the viscous (21) andthermal transport coefficients leads to three—kinetic, collisional, and mixed—contributions. In addition, it was found that for mean free paths λ smallerthan the cell size a , there are finite cell-size corrections which could not be ulti-Particle Collision Dynamics summed in a controlled fashion. The origin of the problem was the explicitappearance of ∆ ξ in the stress correlations. However, it was subsequentlyshown [28,51] that the Green-Kubo relations can be re-summed by introducinga stochastic variable, B i , which is the difference between change in the shiftedcell coordinates of particle i during one streaming step and the actual distancetraveled, ∆t v i . The resulting microscopic stress tensor for the viscous modesis ¯ σ αβ = (cid:88) i (cid:104) mv iα v iβ + m∆t v iα B iβ (cid:105) (22)where B jβ ( t ) = ξ sjβ ( t + ∆t ) − ξ sjβ ( t ) − ∆t v jβ ( t ). It is interesting to comparethis result to the corresponding expression σ αβ = (cid:88) i δ ( r − r i ) mv iα v iβ + 12 (cid:88) j (cid:54) = i r ijα F ijβ ( r ij ) (23)for molecular fluids. The first term in both expressions, the ideal-gas contri-bution, is the same in both cases. The collisional contributions, however, arequit different. The primary reason is that in SRD, the collisional contributioncorresponds to a nonlocal (on the scale of the cell size) force which acts onlyat discrete time intervals. B i has a number of important properties which simplify the calculationof the transport coefficients. In particular, it is shown in Refs. [28, 51] thatstress-stress correlation functions involving one B i in the GK relations forthe transport coefficients are zero, so that, for example, Λ αβ (ˆ k ) = Λ kinαβ (ˆ k ) + Λ colαβ (ˆ k ), with Λ kinαβ (ˆ k ) = ∆tN mk B T ∞ (cid:88) n =0 (cid:48) (cid:104) ˆ k λ σ kinαλ (0) | ˆ k λ (cid:48) σ kinβλ ( n∆t ) (cid:105) (24)and Λ colαβ (ˆ k ) = ∆tN mk B T ∞ (cid:88) n =0 (cid:48) (cid:104) ˆ k λ σ colαλ (0) | ˆ k λ (cid:48) σ colβλ ( n∆t ) (cid:105) ] , (25)with σ kinαβ ( n∆t ) = (cid:88) j mv jα ( n∆t ) v jβ ( n∆t ) (26)and σ colαβ ( n∆t ) = 1 ∆t (cid:88) j mv jα ( n∆t ) B jβ ( n∆t ) , (27)where B jβ ( n∆t ) = ξ sjβ ([ n + 1] ∆t ) − ξ sjβ ( n∆t ) − ∆tv jβ ( n∆t ). Similar relationswere obtained for the thermal diffusivity in Ref. [28]. Analytical calculations of the SRD transport coefficients are greatly simplifiedby the fact that collisional and kinetic contributions to the stress-stress auto-correlation functions decouple. Both the kinetic and collisional contributionshave been calculated explicitly in two and three dimension, and numerous nu-merical tests have shown that the resulting expressions for all the transport co-efficients are in excellent agreement with simulation data. Before summarizingthe results of this work, it is important to emphasize that because of the cellstructure introduced to define coarse-grained collisions, angular momentumis not conserved in a collision [28, 35, 39]. As a consequence, the macroscopicviscous stress tensor is not, in general, a symmetric function of the derivatives ∂ α v β . Although the kinetic contributions to the transport coefficients lead toa symmetric stress tensor, the collisional do not. Before evaluating the trans-port coefficients, we discuss the general form of the macroscopic viscous stresstensor.Assuming only cubic symmetry and allowing for a non-symmetric stresstensor, the most general form of the linearized Navier-Stokes equation is ∂ t v α ( k ) = − ∂ α p + Λ αβ (ˆ k ) v β ( k ) , (28)where Λ αβ ( ˆk ) ≡ ν (cid:18) δ α,β + d − d ˆ k α ˆ k β (cid:19) (29)+ ν (cid:16) δ α,β − ˆ k α ˆ k β (cid:17) + γ ˆ k α ˆ k β + κ ˆ k α δ α,β . In a normal simple liquid, κ = 0 (because of invariance with respect to in-finitesimal rotations) and ν = 0 (because the stress tensor is symmetric in ∂ α v β ), so that the kinematic shear viscosity is ν = ν . In this case, Eq. (29)reduces to the well-known form [20] Λ αβ ( ˆk ) = ν (cid:18) δ α,β + d − d ˆ k α ˆ k β (cid:19) + γ ˆ k α ˆ k β , (30)where γ is the bulk viscosity. Kinetic contributions:
Kinetic contributions to the transport coefficients dom-inate when the mean free path is larger than the cell size, i.e. , λ > a . As canbe seen from Eqs. (24) and (26), an analytic calculations of these contribu-tions requires the evaluation of time correlation functions of products of theparticle velocities. This is straightforward if one makes the basic assumptionof molecular chaos that successive collisions between particles are not corre-lated. In this case, the resulting time-series in Eq. (24) is geometrical, andcan be summed analytically. The resulting expression for the shear viscosityin two dimensions is ν kin = k B T ∆t m (cid:18) M ( M − e − M ) sin ( α ) − (cid:19) . (31) ulti-Particle Collision Dynamics Fig. 1. (a) Normalized kinetic contribution to the viscosity, ν kin / ( ∆tk B T ), in threedimensions as a function of the collision angle α . Data were obtained by time aver-aging the Green-Kubo relation over 75,000 iterations using λ/a = 2 .
309 for M = 5( (cid:4) ) and M = 20 ( • ). The lines are the theoretical prediction, Eq. (32). Parameters: L/a = 32, ∆t = 1. From Ref. [53].(b) Normalized collisional contribution to the viscosity, ν col ∆t/a , in three dimen-sions as a function of the collision angle α . The solid line is the theoretical predic-tion, Eq. (39). Data were obtained by time averaging the Green-Kubo relation over300,000 iterations. Parameters: L/a = 16, λ/a = 0 . M = 3, and ∆t = 1. FromRef. [54]. Fluctuations in the number of particles per cell are included in (31). Thisresult agrees with the non-equilibrium calculations of Pooley et al. [35, 52],measurements in shear flow [26], and the numerical evaluation of the GKrelation in equilibrium simulations (see Fig. 1).The corresponding result in three dimensions for collision rule (3) is ν kin = k B T ∆t m (cid:18) M ( M − e − M )[2 − cos( α ) − cos(2 α )] − (cid:19) . (32)The kinetic contribution to the stress tensor is symmetric, so that ν kin = 0and the kinetic contribution to the shear viscosity is ν kin ≡ ν kin . Collisional contributions:
Explicit expressions for the collisional contributionsto the viscous transport coefficients can be obtained by considering variouschoices for ˆk and α and β in Eqs. (25), (27) and (29). Taking ˆk in the y -direction and α = β = 1 yields ν col + ν col = 1 ∆tN k B T ∞ (cid:88) t =0 (cid:48) (cid:88) i,j (cid:104) v ix (0) B iy (0) v ix ( t ) B iy ( t ) (cid:105) . (33)Other choices lead to relations between the collisional contributions to theviscous transport coefficients, namely[1 + ( d − /d ] ν col + γ col + κ col = ν col + ν col . (34)and[( d − /d ] ν col − ν col + γ col = 0 . (35) These results imply that κ col = 0, and γ col − ν col /d = ν col − ν col . It followsthat the collision contribution to the macroscopic viscous stress tensor isˆ σ colαβ /ρ = ν col ( ∂ β v α + ∂ α v β ) + ν col ( ∂ β v α − ∂ α v β ) + ( ν col − ν col ) δ αβ ∂ λ v λ = ( ν col + ν col ) ∂ β v α + ( ν col − ν col ) Q αβ , (36)where Q αβ ≡ δ αβ ∂ λ v λ − ∂ α v β . Since Q αβ has zero divergence, ∂ β Q αβ = 0, theterm containing Q in Eq. (36) will not appear in the linearized hydrodynamicequation for the momentum density, so that ρ ∂ v ∂t = −∇ p + ρ ( ν kin + ν col ) ∆ v + d − d ν kin ∇ ( ∇ · v ) , (37)where ν col = ν col + ν col . In writing Eq. (37) we have used the fact that thekinetic contribution to the microscopic stress tensor, ¯ σ kin , is symmetric, and γ kin = 0 [27]. The viscous contribution to the sound attenuation coefficient is ν col + 2( d − ν kin /d instead of the standard result, 2( d − ν/d + γ , for simpleisotropic fluids. The collisional contribution to the effective shear viscosity is ν col ≡ ν col + ν col . It is interesting to note that the kinetic theory approachdiscussed in Ref. [35] is able to show explicitly that ν col = ν col , so that ν col = 2 ν col .It is straightforward to evaluate the various contributions to the right handside of (33). In particular, note that since velocity correlation functions areonly required at equal times and for a time lag of one time step, molecularchaos can be assumed [51]. Using the relation [28] (cid:104) B iα ( n∆t ) B jβ ( m∆t ) (cid:105) = a δ αβ (1 + δ ij ) [2 δ n,m − δ n,m +1 − δ n,m − ] , (38)and averaging over the number of particles in a cell assuming that the numberof particles in any cell is Poisson distributed at each time step, with an averagenumber M of particles per cell, one then finds ν col = ν col + ν col = a d∆t (cid:18) M − e − M M (cid:19) [1 − cos( α )] , (39)for the SRD collision rules in both two and three dimensions. Eq. (39) agreeswith the result of Refs. [26] and [35] obtained using a completely differentnon-equilibrium approach in shear flow. Simulation results for the collisionalcontribution to the viscosity are in excellent agreement with this result (seeFig. 1). Thermal diffusivity and self-diffusion coefficient:
As with the viscosity, thereare both kinetic and collisional contributions to the thermal diffusivity, D T .A detailed analysis of both contributions is given in Ref. [28], and the resultsare summarized in Table 1. The self-diffusion coefficient, D , of particle i isdefined by D = lim t →∞ dt (cid:104) [ r i ( t ) − r i (0)] (cid:105) = ∆td ∞ (cid:88) n =0 (cid:48) (cid:104) v i ( n∆t ) · v i (0) (cid:105) , (40) ulti-Particle Collision Dynamics where the second expressions is the corresponding discrete GK relation. Theself-diffusion coefficient is unique in that the collisions do not explicitly con-tribute to D . With the assumption of molecular chaos, the kinetic contribu-tions are easily summed [27] to obtain the result given in Table 1. The kinetic contributions to the transport coefficients presented in Table 1have all been derived under the assumption of molecular chaos, i.e., thatparticle velocities are not correlated. Simulation results for the shear viscosityand thermal diffusivity have generally been found to be in good agreementwith these results. However, it is known that there are correlation effects for λ/a smaller than unity [15, 55]. They arise from correlated collisions betweenparticles that are in the same collision cell for more than one time step.For the viscosity and thermal conductivity, these corrections are generallynegligible, since they are only significant in the small λ/a regime, where thecollisional contribution to the transport coefficients dominates. In this regard,it is important to note that there are no correlation corrections to ν col and D colT [28]. For the self-diffusion coefficient—for which there is no collisionalcontribution—correlation corrections dramatically increase the value of thistransport coefficient for λ (cid:28) a , see Refs. [15,55]. These correlation corrections,which arise from particles which collide with the same particles in consecutivetime steps, are distinct from the correlations effects which are responsible forthe long-time tails. This distinction is important, since long-time tails are alsovisible at large mean free paths, where these corrections are negligible. Spontaneous thermal fluctuations of the density, ρ ( r , t ), the momentum den-sity, g ( r , t ), and the energy density, (cid:15) ( r , t ), are dynamically coupled, and ananalysis of their dynamic correlations in the limit of small wave numbers andfrequencies can be used to measure a fluid’s transport coefficients. In par-ticular, because it is easily measured in dynamic light scattering, x-ray, andneutron scattering experiments, the Fourier transform of the density-densitycorrelation function—the dynamics structure factor—is one of the most widelyused vehicles for probing the dynamic and transport properties of liquids [56].A detailed analysis of equilibrium dynamic correlation functions—the dy-namic structure factor as well as the vorticity and entropy-density correlationfunctions—using the SRD algorithm is presented in Ref. [57]. The results—which are in good agreement with earlier numerical measurements and the-oretical predictions—provided further evidence that the analytic expressionsor the transport coefficients are accurate and that we have an excellent un-derstanding of the SRD algorithm at the kinetic level. d Kinetic ( × k B T ∆t/ m ) Collisional ( × a /∆t ) ν M ( M − e − M ) sin ( α ) − M − e − M ) dM [1 − cos( α )]3 M ( M − e − M )[2 − cos( α ) − cos(2 α )] − D T d − cos( α ) − dM ˆ − d − csc ( α/ ˜ (1 − /M )3( d +2) M [1 − cos( α )]3 D dM (1 − cos( α ))( M − e − M ) − Table 1.
Theoretical expressions for the kinematic shear viscosity ν , the thermaldiffusivity, D T , and the self-diffusion coefficient, D , in both two ( d = 2) and three( d = 3) dimensions. M is the average number of particles per cell, α is the collisionangle, k B is Boltzmann’s constant, T is the temperature, ∆t is the time step, m isthe particle mass, and a is the cell size. Except for self-diffusion constant, for whichthere is no collisional contribution, both the kinetic and collisional contributions arelisted. The expressions for shear viscosity and self-diffusion coefficient include theeffect of fluctuations in the number of particles per cell; however, for brevity, therelations for thermal diffusivity are correct only up to O (1 /M ) and O (1 /M ) for thekinetic and collisional contributions, respectively. For the complete expressions, seeRefs. [28, 53, 54]. Here, we briefly summarize the results for the dynamic structure factor.The dynamic structure exhibits three peaks, a central “Rayleigh” peak causedby the thermal diffusion, and two symmetrically placed “Brillouin peaks”caused by sound. The width of the central peak is determined by the thermaldiffusivity, D T , while that of the two Brillouin peaks is related to the soundattenuation coefficient, Γ . For the SRD algorithm [57], Γ = D T (cid:18) c p c v − (cid:19) + 2 (cid:18) d − d (cid:19) ν kin + ν col . (41)Note that in two-dimensions, the sound attenuation coefficient for a SRD fluidhas the same functional dependence on D T and ν = ν kin + ν col as an isotropicfluid with an ideal-gas equation of state (for which γ = 0).Simulation results for the structure factor in two-dimensions with λ/a =1 . α = 120 ◦ , and λ/a = 0 . α = 60 ◦ are shown in Figs. 2a and 2b, respectively. The solid lines are the theoreticalprediction for the dynamic structure factor (see Eq. (36) of Ref. [57]) using c = (cid:112) k B T /m and values for the transport coefficients obtained using theexpressions in Table 1, assuming that the bulk viscosity γ = 0. As can beseen, the agreement is excellent. ulti-Particle Collision Dynamics Fig. 2.
Normalized dynamic structure, S cρρ ( kω ) /χ ρρ ( k ), for k = 2 π (1 , /L and (a) λ/a = 1 . α = 120 ◦ , and (b) λ/a = 0 . α = 60 ◦ . The solid lines are thetheoretical prediction for the dynamic structure factor (see Eq. (36) of Ref. [57])using values for the transport coefficients obtained using the expressions in Table 1.The dotted lines show the predicted positions of the Brillouin peaks, ω = ± ck , with c = p k B T /m . Parameters:
L/a = 32, M = 15, and ∆t = 1 .
0. From Ref. [57].
MPC transport coefficients have also be evaluated by calculating the linearresponse of the system to imposed gradients. This approach was introducedby Kikuchi et al. [26] for the shear viscosity and then extended and refinedin Ref. [35] to determine the thermal diffusivity and bulk viscosity. Here, wesummarize the derivation of the shear viscosity.
Linear response theory provides an alternative, and complementary, approachfor evaluating the shear viscosity. This non-equilibrium approach is relatedto equilibrium calculations described in the previous section through thefluctuation-dissipation theorem. Both methods yield identical results. For themore complicated analysis of the hydrodynamic equations, the stress tensor,and the longitudinal transport coefficients such as the thermal conductivity,the reader is referred to Ref. [35].Following Kikuchi et al. [26], we consider a two-dimensional liquid withan imposed shear ˙ γ = ∂u x ( y ) /∂y . On average, the velocity profile is givenby v = ( ˙ γy, η is the proportionality constantbetween the velocity gradient ˙ γ and the frictional force acting on a planeperpendicular to y ; i.e. σ xy = η ˙ γ, (42)where σ xy is the off-diagonal element of the viscous stress tensor. During thestreaming step, particles will cross this plane only if | v y ∆t | is greater than the distance to the plane. Assuming that the fluid particles are homogeneouslydistributed, the momentum flux is obtained by integrating over the coordi-nates and velocities of all particles that cross the plane from above and belowduring the time step ∆t . The result is [26] σ xy = ρ (cid:18) ˙ γ∆t (cid:104) v y (cid:105) − (cid:104) v x v y (cid:105) (cid:19) , (43)where the mass density ρ = mM/a d , and the averages are taken over thesteady-state distribution P ( v x − ˙ γy, v y ). It is important to note that this is not the Maxwell-Boltzmann distribution, since we are in a non-equilibriumsteady state where the shear has induced correlations between v x and v y . Asa consequence, (cid:104) v x v y (cid:105) is nonzero. To determine the behavior of (cid:104) v x v y (cid:105) , theeffect of streaming and collisions are calculated separately. During streaming,particles which arrive at y with positive velocity v y have started from y − v y ∆t ; these particles bring a velocity component v x which is smaller than thatof particles originally located at y . On the other hand, particles starting out at y > y with negative v y bring a larger v x . The velocity distribution is thereforesheared by the streaming, so that P after ( v x , v y ) = P before ( v x + ˙ γv y ∆t, v y ).Averaging v x v y over this distribution gives [26] (cid:104) v x v y (cid:105) after = (cid:104) v x v y (cid:105) − ˙ γ∆t (cid:104) v y (cid:105) , (44)where the superscript denotes the quantity after streaming. The streamingstep therefore reduces correlations by − ˙ γ∆t (cid:104) v y (cid:105) , making v x and v y increas-ingly anti-correlated.The collision step redistributes momentum between particles and tends toreduce correlations. Making the assumption of molecular chaos, i.e., that isthat the velocities of different particles are uncorrelated, and averaging overthe two possible rotation directions, one finds, (cid:104) v x v y (cid:105) after = (cid:20) − N c − N c [1 − cos(2 α )] (cid:21) (cid:104) v x v y (cid:105) before (45)The number of particles in a cell, N c is not constant, and density fluctuationshave to be included. The probability to find n uncorrelated particles in agiven cell is given by the Poisson distribution, w ( n ) = exp( − M ) M n /n !; theprobability of a given particle being in a cell together with n − nw ( n ) /M . Taking an average over this distribution gives (cid:104) v x v y (cid:105) after = f (cid:104) v x v y (cid:105) before , (46)with f = (cid:20) − M − − M ) M [1 − cos(2 α )] (cid:21) . (47)The difference between this result and just replacing N c by M in Eq. (45) issmall, and only important for M ≤
3. One sees that (cid:104) v x v y (cid:105) is first modified bystreaming and then multiplied by a factor f in the subsequent collision step. ulti-Particle Collision Dynamics In the steady state, it therefore oscillates between two values. Using Eqs. (44),(46), and (47), we obtain the self-consistency condition ( (cid:104) v x v y (cid:105)− ˙ γ∆t (cid:104) v y (cid:105) ) f = (cid:104) v x v y (cid:105) . Solving for (cid:104) v x v y (cid:105) , assuming equipartition of energy, (cid:104) v y (cid:105) = k B T /m ,and substituting into (43), we have σ xy = ˙ γ M ∆tk B Tm (cid:18)
12 + f − f (cid:19) , (48)Inserting this result into the definition of the viscosity, (42), yields the sameexpression for the kinetic viscosity in two-dimensional as obtained by theequilibrium Green-Kubo approach discussed in Sec. 4.1.1. The collisional contribution to the shear viscosity is proportional to a /∆t ; asdiscussed in Sec. 3.2, it results from the momentum transfer between particlesin a cell of size a during the collision step. Consider again a collision cell oflinear dimension a with a shear flow u x ( y ) = ˙ γy . Since the collisions occur ina shifted grid, they cause a transfer of momentum between neighboring cellsof the original unshifted reference frame [21, 27]. Consider now the momen-tum transfer due to collisions across the line y = h , the coordinate of a cellboundary in the unshifted frame. If we assume a homogeneous distributionof particles in the collision cell, the mean velocities in the upper ( y > h ) andlower partitions are u = 1 M M (cid:88) i =1 v i and u = 1 M M (cid:88) i = M +1 v i , (49)respectively, where M = M ( a − h ) /a and M = M h/a . Collisions trans-fer momentum between the two parts of the cell. The x -component of themomentum transfer is ∆p x ( h ) ≡ M (cid:88) i =1 [ v afterix − v beforeix ] . (50)The use of the rotation rule (2) together with an average over the sign of thestochastic rotation angle yields ∆p x ( h ) = [cos( α ) − M ( u x − u x ) . (51)Since M u = M u + M u , ∆p x ( h ) = [1 − cos( α )] M ( u x − u x ) ha (cid:18) − ha (cid:19) . (52)Averaging over the position h of the dividing line, which corresponds to av-eraging over the random shift, we find (cid:104) ∆p x (cid:105) = 1 a (cid:90) a ∆p x ( h ) dh = 16 [1 − cos( α )] M ( u x − u x ) . (53)Since the dynamic viscosity η is defined as the ratio of the tangential stress, P yx , to ∂u x /∂y , we have η = (cid:104) ∆p x (cid:105) / ( a ∆t ) ∂u x /∂y = (cid:104) ∆p x (cid:105) / ( a ∆t )( u x − u x ) / ( a/ , (54)so that the kinematic viscosity, ν = η/ρ , in two-dimensions for SRD is ν col = a ∆t [1 − cos( α )] (55)in the limit of small mean free path. Since we have neglected the fluctuationsin the particle number, this expression corresponds to the limit M → ∞ . Eventhough this derivation is somewhat heuristic, it gives a remarkably accurateexpression; in particular, it contains the correct dependence on the cell size, a , and the time step, ∆t , in the limit of small free path, ν col = a ∆t f col ( d, M, α ) , (56)as expected from simple random walk arguments. Kikuchi et al. [26] includedparticle number fluctuations and obtained identical results for the collisionalcontribution to the viscosity as was obtained in the Green-Kubo approach(see Table 1). For MPC-AT, the viscosities have been calculated in Ref. [32] using the meth-ods described in Secs. 5.1 and 5.2. The total viscosity of MPC-AT is given bythe sum of two terms, the collisional and kinetic contributions. For MPC-AT-a, it was found for both two and three dimensions that [32] ν kin = k B T ∆tm (cid:18) MM − − M − (cid:19) and ν col = a ∆t (cid:18) M − − M M (cid:19) . (57)The exponential terms e − M are due to the fluctuation of the particle num-ber per cell and become important for M ≤
3. As was the case for SRD,the kinetic viscosity has no anti-symmetric component; the collisional con-tribution, however, does. Again, as discussed in Sec. 4.1.1 for SRD, one finds ν col = ν col = ν col /
2. This relation is true for all − a versions of MPC discussedin Refs. [32, 58, 59]. Simulation results were found to be in good agreementwith theory.For MPC-AT+a it was found for sufficiently large M that [38, 59] ulti-Particle Collision Dynamics ν kin = k B T ∆tm (cid:18) MM − ( d + 2) / − (cid:19) ,ν col = a ∆t (cid:18) M − / M (cid:19) . (58)MPC-AT-a and MPC-AT+a both have the same kinetic contribution to theviscosity in two dimensions; however, imposing angular-momentum conserva-tion makes the collisional contribution to the stress tensor symmetric, so thatthe asymmetric contribution, ν , discussed in Sec. 4.1.1 vanishes. The result-ing collisional contribution to the viscosity is then reduced by a factor closeto two. The original SRD algorithm models a single-component fluid with an ideal-gas equation of state. The fluid is therefore very compressible, and the speedof sound, c s , is low. In order to have negligible compressibility effects, as inreal liquids, the Mach number has to be kept small, which means that thereare limits on the flow velocity in the simulation. The SRD algorithm can bemodified to model both excluded volume effects, allowing for a more realisticmodeling of dense gases and liquids, as well as repulsive hard-core interac-tions between components in mixtures, which allow for a thermodynamicallyconsistent modeling of phase separating mixtures. As in SRD, the algorithm consists of individual streaming and collision steps.In order to define the collisions, a second grid with sides of length 2 a is intro-duced, which (in d = 2) groups four adjacent cells into one “supercell”. Thecell structure is sketched in Fig. 3 (left panel). To initiate a collision, pairsof cells in every supercell are chosen at random. Three different choices arepossible: a) horizontal (with σ = ˆ x ), b) vertical ( σ = ˆ y ), and c) diagonalcollisions (with σ = (ˆ x + ˆ y ) / √ σ = (ˆ x − ˆ y ) / √ u n = (1 /M n ) (cid:80) M n i =1 v i , of cell n , the projec-tion of the difference of the mean velocities of the selected cell pairs on σ j , ∆u = σ j · ( u − u ), is then used to determine the probability of collision. If ∆u <
0, no collision will be performed. For positive ∆u , a collision will occurwith an acceptance probability, p A , which depends on ∆u and the numberof particles in the two cells, N and N . The choice of p A determines boththe equation of state and the values of the transport coefficients. While thereis considerable freedom in choosing p A , the requirement of thermodynamicconsistency imposes certain restrictions [33, 34, 55]. One possible choice is Fig. 3. (left panel) Schematic of collision rules. Momentum is exchanged in threeways: a) horizontally along σ , b) vertically along σ , and c) diagonally along σ and σ . w and w d denote the probabilities of choosing collisions a), b), and c)respectively.(right panel) Static structure factor S (¯ k, t = 0) as a function of ∆t for M = 3.The open circles ( ◦ ) show results obtained by taking the numerical derivative ofthe pressure. The bullets ( • ) are data obtained from direct measurements of thedensity fluctuations. The solid line is the theoretical prediction obtained using thefirst term in Eq. (61) and Eq. (63). ¯ k is the smallest wave vector, ¯ k = (2 π/L )(1 , L/a = 32, A = 1 /
60, and k B T = 1 .
0. From Ref. [33]. p A ( M , M , ∆u ) = Θ ( ∆u ) tanh( Λ ) with Λ = A ∆u N N , (59)where Θ is the unit step function and A is a parameter which is used to tunethe equation of state. The choice Λ ∼ N N leads to a non-ideal contributionto the pressure which is quadratic in the particle density.The collision rule chosen in Ref. [33] maximizes the momentum trans-fer parallel to the connecting vector σ j and does not change the transversemomentum. It exchanges the parallel component of the mean velocities ofthe two cells, which is equivalent to a “reflection” of the relative velocities, v (cid:107) i ( t + ∆t ) − u (cid:107) = − ( v (cid:107) i ( t ) − u (cid:107) ), where u (cid:107) is the parallel component of themean velocity of the particles of both cells. This rule conserves momentumand energy in the cell pairs.Because of x − y symmetry, the probabilities for choosing cell pairs in the x - and y - directions (with unit vectors σ and σ in Fig. 3) are equal, andwill be denoted by w . The probability for choosing diagonal pairs ( σ and σ in Fig. 3) is given by w d = 1 − w . w and w d must be chosen so that thehydrodynamic equations are isotropic and do not depend on the orientation ofthe underlying grid. An equivalent criterion is to guarantee that the relaxationof the velocity distribution is isotropic. These conditions require w = 1 / w d = 1 /
2. This particular choice also ensures that the kinetic part of theviscous stress tensor is isotropic [45]. ulti-Particle Collision Dynamics The transport coefficients can be determined using the same Green-Kubo for-malism as was used for the original SRD algorithm [21, 51]. Alternatively, thenon-equilibrium approach describe in Sec. 5 can be used. Assuming molecularchaos and ignoring fluctuations in the number of particles per cell, the kineticcontribution to the viscosity is found to be ν kin = k B Tm ∆t (cid:18) p col − (cid:19) with p col = A (cid:114) k B Tmπ M / , (60)which is in good agreement with simulation data. p col is essentially the collisionrate, and can be obtained by averaging the acceptance probability, Eq. (59).The collisional contribution to the viscosity is ν col = p col ( a / ∆t ) [60].The self-diffusion constant, D , is evaluated by summing over the velocity-autocorrelation function (see, e.g. Ref. [21]); which yields D = ν kin . The collision rules conserve the kinetic energy, so that the internal energyshould be the same as that of an ideal gas. Thermodynamic consistency there-fore requires that the non-ideal contribution to the pressure is linear in T . Thisis possible if the coefficient A in Eq. (59) is sufficiently small.The mechanical definition of pressure—the average longitudinal momen-tum transfer across a fixed interface per unit time and unit surface area—canbe used to determine the equation of state. Only the momentum transferdue to collisions needs to be considered, since that coming from streamingconstitutes the ideal part of the pressure. Performing this calculation for afixed interface and averaging over the position of the interface, one finds thenon-ideal part of the pressure, P n = (cid:18) √ (cid:19) A M k B Ta∆t + O ( A T ) . (61) P n is quadratic in the particle density, ρ = M/a , as it would be expectedfrom a virial expansion. The prefactor A must be chosen small enough thathigher order terms in this expansion are negligible. Prefactors A leading toacceptance rates of about 15% are sufficiently small to guarantee that thepressure is linear in T .The total pressure is the average of the diagonal part of the microscopicstress tensor, P = P id + P n = 1 ∆tL x L y (cid:42)(cid:88) j { ∆tv jx − ∆v jx z sjlx / } (cid:43) . (62) The first term gives the ideal part of the pressure, P id , as discussed in Ref. [21].The average of the second term is the non-ideal part of the pressure, P n . z sjl is avector which indexes collision partners. The first subscript denotes the particlenumber and the second, l , is the index of the collision vectors σ l in Fig. 3 (leftpanel). The components of z sjl are either 0, 1, or − P n obtained using Eq. (62) are in good agreement with the analyticalexpression, Eq. (61). In addition, measurements of the static structure factor S ( k → , t = 0) agree with the thermodynamic prediction S ( k → , t = 0) = ρk B T ∂ρ/∂P | T (63)when result (61) is used [see Fig. 3 (right panel)]. The adiabatic speed of soundobtained from simulations of the dynamic structure factor is also in goodagreement with the predictions following from Eq. (61). These results providestrong evidence for the thermodynamic consistency of the model. Consistencychecks are particularly important because the non-ideal algorithm does notconserve phase-space volume. This is because the collision probability dependson the difference of collision-cell velocities, so that two different states can bemapped onto the same state by a collision. While the dynamics presumablystill obeys detailed—or at least semi-detailed—balance, this is very hard toprove, since it would require knowledge not only of the transition probabilities,but also of the probabilities of the individual equilibrium states. Nonetheless,no inconsistencies due to the absence of time-reversal invariance or a possibleviolation of detailed balance have been observed.The structure of S ( k ) for this model is also very similar to that of a simpledense fluid. In particular, for fixed M , both the depth of the minimum atsmall k and the height of the first peak increase with decreasing ∆t , untilthere is an order-disorder transition. The four-fold symmetry of the resultingordered state—in which clusters of particles are concentrated at sites with theperiodicity close, but not necessarily equal, to that of the underlying grid—is clearly dictated by the structure of the collision cells. Nevertheless, theseordered structures are similar to the low-temperature phase of particles witha strong repulsion at intermediate distances, but a soft repulsion at shortdistances. The scaling behavior of both the self-diffusion constant and thepressure persists until the order/disorder transition. In a binary mixture of A and B particles, phase separation can occur whenthere is an effective repulsion between A-B pairs. In the current model, this isachieved by introducing velocity-dependent multi-particle collisions betweenA and B particles. There are N A and N B particles of type A and B, respec-tively. In two dimensions, the system is coarse-grained into ( L/a ) cells of asquare lattice of linear dimension L and lattice constant a . The generalizationto three dimensions is straightforward. ulti-Particle Collision Dynamics Collisions are defined in the same way as in the non-ideal model discussedin the previous section. Now, however, two types of collisions are possible foreach pair of cells: particles of type A in the first cell can undergo a collisionwith particles of type B in the second cell; vice versa, particles of type B inthe first cell can undergo a collision with particles of type A in the secondcell. There are no A-A or B-B collisions, so that there is an effective repulsionbetween A-B pairs. The rules and probabilities for these collisions are chosen inthe same way as in the non-ideal single-component fluid described in Refs. [33,55]. For example, consider the collision of A particles in the first cell with theB particles in the second. The mean particle velocity of A particles in thefirst cell is u A = (1 /N c,A ) (cid:80) N c,A i =1 v i , where the sum runs over all A particles, N c,A , in the first cell. Similarly, u B = (1 /N c,B ) (cid:80) N c,B i =1 v i is the mean velocityof B particles in the second cell. The projection of the difference of the meanvelocities of the selected cell-pairs on σ j , ∆u AB = σ j · ( u A − u B ), is thenused to determine the probability of collision. If ∆u AB <
0, no collision willbe performed. For positive ∆u AB , a collision will occur with an acceptanceprobability p A ( N c,A , N c,B , ∆u AB ) = A ∆u AB Θ ( ∆u AB ) N c,A N c,B , (64)where Θ is the unit step function and A is a parameter which allows us totune the equation of state; in order to ensure thermodynamic consistency, itmust be sufficiently small that that p A < v (cid:107) i ( t + ∆t ) − u (cid:107) AB = − ( v (cid:107) i ( t ) − u (cid:107) AB ), is exchanged,where u (cid:107) AB = ( N c,A u (cid:107) A + N c,B u (cid:107) B ) / ( N c,A + N c,B ) is the parallel componentof the mean velocity of the colliding particles. The perpendicular componentremains unchanged. It is easy to verify that these rules conserve momentumand energy in the cell pairs. The collision of B particles in the first cell withA particles in the second is handled in a similar fashion.Because there are no A-A and B-B collisions, additional SRD collisions atthe cell level are incorporated in order to mix particle momenta. The order ofA-B and SRD collision is random, i.e., the SRD collision is performed first witha probability 1 /
2. If necessary, the viscosity can be tuned by not performingSRD collisions every time step. The results presented here were obtained usinga SRD collision angle of α = 90 ◦ .The transport coefficients can be calculated in the same way as for theone-component non-ideal system. The resulting kinetic contribution to theviscosity is ν kin = ∆tk B T (cid:18) A (cid:114) πk B T [ M A M B ( M A + M B )] − / − (cid:19) , (65)where M A = (cid:104) N c,A (cid:105) , M B = (cid:104) N c,B (cid:105) . In deep quenches, the concentrationof the minority component is very small, and the non-ideal contribution tothe viscosity approaches zero. In this case, the SRD collisions provide thedominant contribution to the viscosity. An analytic expression for the equation of state of this model can be derivedby calculating the momentum transfer across a fixed surface, in much thesame way as was done for the non-ideal model in Ref. [33]. Since there areonly non-ideal collisions between A-B particles, the resulting contribution tothe pressure is P n = (cid:18) w + w d √ (cid:19) AM A M B k B Ta∆t = Γ ρ A ρ B , (66)where ρ A and ρ B are the densities of A and B and Γ ≡ ( w + w d / √ a A/∆t .In simulations, the total pressure can be measured by taking the ensembleaverage of the diagonal components of the microscopic stress tensor. In thisway, the pressure can be measured locally, at the cell level. In particular, thepressure in a region consisting of N cell cells is P n = 1 ∆ta N cell (cid:42) N c (cid:88) c =1 (cid:88) j ∈ c [ ∆tv jx − ∆v jx z sjlx / (cid:43) , (67)where the second sum runs over the particles in cell c . The first term inEq. (67) is the ideal-gas contribution to the pressure; the second comes fromthe momentum transfer between cells involved in the collision indexed by z sjl [45].Expression (66) can be used to determine the entropy density, s . Theideal-gas contribution to s has the form [61] s ideal = ρ ϕ ( T ) − k B [ ρ A ln ρ A + ρ B ln ρ B ] , (68)where ρ = ρ A + ρ B . Since ϕ ( T ) is independent of ρ A and ρ B , this term doesnot play a role in the current discussion. The non-ideal contribution to theentropy density, s n , can be obtained from Eq. (66) using the thermodynamicrelation P n /T = − s n + ρ A ∂s n /∂ρ A + ρ B ∂s n /∂ρ B . (69)The result is s n = Γ ρ A ρ B , so that the total configurational contribution tothe entropy density is s = − k B { ρ A ln ρ A + ρ B ln ρ B + Γ ρ A ρ B } . (70)Since there is no configurational contribution to the internal energy in thismodel, the mean-field phase diagram can be determined by maximizing theentropy at fixed density ρ . The resulting demixing phase diagram as a functionof ρ AB = ( ρ A − ρ B ) /ρ is given by the solid line in Fig. 4 (left panel). Thecritical point is located at ρ AB = 0, ρΓ ∗ = 2. For ρΓ <
2, the order parameter ρ AB = 0; for ρΓ >
2, there is phase separation into coexisting A- and B-richphases. As can be seen, the agreement between the mean-field predictions andsimulation results is very good except close to the critical point, where the ulti-Particle Collision Dynamics Fig. 4. (left panel) Binary phase diagram. There is phase separation for ρΓ >
2. Sim-ulation results for ρ AB obtained from concentration histograms are shown as bullets.The dashed line is a plot of the leading singular behavior, ρ AB = p ρΓ − / ρ AB = 0 to ρΓ = 3 .
62 (arrow). The dark (blue)and light (white) spheres are A and B particles, respectively. Parameters:
L/a = 64, M A = M B = 5, k B T = 0 . ∆t = 1, and a = 1. From Ref. [45].(right panel) Dimensionless radial fluctuations, (cid:104)| u k |(cid:105) , as a function of the modenumber k for A = 0 .
45 ( • ) and A = 0 .
60 ( ◦ ) with k B T = 0 . r = 11 . a and r = 15 . a , respectively. The solid lines arefits to Eq. (71). The inset shows a typical droplet configuration for ρ AB = − . ρΓ = 3 .
62 ( A = 0 .
60 and k B T = 0 . L/a = 64, M A = 2, M B = 8, ∆t = 1, and a = 1. From Ref. [45]. histogram method of determining the coexisting densities is unreliable andcritical fluctuations influence the shape of the coexistence curve. A typical configuration for ρ AB = 0, ρΓ = 3 .
62 is shown in the inset to Fig. 4(left panel), and a snapshot of a fluctuating droplet at ρ AB = − . ρΓ = 3 . σ . Usingthe parameterization r ( φ ) = r (cid:2) (cid:80) ∞ k = −∞ u k exp( ikφ ) (cid:3) and choosing u tofix the area of the droplet, it can be shown that [54] (cid:104)| u k | (cid:105) = k B T πr σ (cid:18) k − (cid:19) . (71)Figure 4 (right panel) contains a plot of (cid:104)| u k | (cid:105) as a function of mode number k for ρΓ = 3 .
62 and ρΓ = 2 .
72. Fits to the data yield σ (cid:39) . k B T for ρΓ =3 .
62 and σ (cid:39) . k B T for ρΓ = 2 .
72. Mechanical equilibrium requires that the pressure difference across the interface of a droplet satisfies the Laplaceequation ∆p = p in − p out = ( d − σ/r (72)in d spatial dimensions. Measurements of ∆p [using Eq. (67)] as a functionof the droplet radius for A = 0 .
60 at k B T = 0 . Fig. 5. (left panel) Droplet configuration in a mixture with N A = 8192, N B = 32768,and N d = 1500 dimers after 10 time steps. The initial configuration is a dropletwith a homogeneous distribution of dimers. The dark (blue) and light (white) coloredspheres indicate A and B particles, respectively. For clarity, A particles in the bulkare smaller and B particles in the bulk are not shown. Parameters: L/a = 64, M A = 2, M B = 8, A = 1 . k B T = 0 . ∆t = 1, and a = 1.(right panel) Typical configuration showing the bicontinuous phase for N A = N B =20480 and N d = 3000. Parameters: L/a = 64, M A = 5, M B = 5, A = 1 . k B T =0 . ∆t = 1, and a = 1. From Ref. [45]. The model therefore displays the correct thermodynamic behavior andinterfacial fluctuations. It can also be extended to model amphiphilic mixturesby introducing dimers consisting of tethered A and B particles. If the A and Bcomponents of the dimers participate in the same collisions as the solvent, theybehave like amphiphilic molecules in binary oil-water mixtures. The resultingmodel displays a rich phase behavior as a function of ρΓ and the number ofdimers, N d . Both the formation of droplets and micelles, as shown in Fig. 5(left panel), and a bicontinuous phase, as illustrated in Fig. 5 (right panel),have been observed [45]. The coarse-grained nature of the algorithm thereforeenables the study of large time scales with a feasible computational effort. ulti-Particle Collision Dynamics There have been other generalizations of SRD to model binary mixtures byHashimoto et al. [62] and Inoue et al. [63], in which a color charge, c i = ± α in theSRD rotation step is then chosen such that the color-weighted momentum ina cell, m = (cid:80) N c i =1 c i ( v i − u ), is rotated to point in the direction of the gradientof the color field ¯ c = (cid:80) N c i =1 c i . This rule also leads to phase separation. Sev-eral tests of the model have been performed; Laplace’s equation was verifiednumerically, and simulation studies of spinodal decomposition and the defor-mation of a falling droplet were performed [62]. Later applications include astudy of the transport of slightly deformed immiscible droplets in a bifurcatingchannel [64]. Subsequently, the model was generalized through the additionof dumbbell-shaped surfactants to model micellization [65] and the behaviorof ternary amphiphilic mixtures in both two and three dimensions [66, 67].Note that since the color current after the collision is always parallel to thecolor gradient, thermal fluctuations of the order parameter are neglected inthis approach. A very simple procedure for coupling embedded objects such as colloids orpolymers to a MPC solvent has been proposed in Ref. [68]. In this approach,every colloid particle or monomer in the polymer chain is taken to be a point-particle which participates in the SRD collision. If monomer i has mass m m and velocity w i , the center of mass velocity of the particles in the collisioncell is u = m (cid:80) N c i =1 v i + m m (cid:80) N m i =1 w i N c m + N m m m , (73)where N m is the number of monomers in the collision cell. A stochastic col-lision of the relative velocities of both the solvent particles and embeddedmonomers is then performed in the collision step. This results in an exchangeof momentum between the solvent and embedded monomers. The same proce-dure can of course be employed for other MPC algorithms, such as MPC-AT.The new monomer momenta are then used as initial conditions for a molecular-dynamics update of the polymer degrees of freedom during the subsequentstreaming time step, ∆t . Alternatively, the momentum exchange, ∆p , can beincluded as an additional force ∆p/∆t in the molecular-dynamics integration.If there are no other interactions between monomers—as might be the case for embedded colloids—these degrees of freedom stream freely during this timeinterval.When using this approach, the average mass of solvent particles per cell, mN c , should be of the order of the monomer or colloid mass m m (assumingone embedded particle per cell) [15]. This corresponds to a neutrally buoyantobject which responds quickly to the fluid flow but is not kicked around tooviolently. It is also important to note that the average number of monomersper cell, (cid:104) N m (cid:105) , should be smaller than unity in order to properly resolvehydrodynamic interactions between the monomers. On the other hand, theaverage bond length in a semi-flexible polymer or rod-like colloid should alsonot be much larger than the cell size a , in order to capture the anisotropicfriction of rod-like molecules due to hydrodynamic interactions [69] (whichleads to a twice as large perpendicular than parallel friction coefficient forlong stiff rods [70]), and to avoid an unnecessarily large ratio of the numberof solvent to solute particles. For a polymer, the average bond length shouldtherefore be of the order of a .In order to use SRD to model suspended colloids with a radius of order1 µ m in water, this approach would require approximately 60 solvent particlesper cell in order to match the Peclet number [16]. This is much larger thanthe optimum number (see discussion in Sec. 2.2.1), and the relaxation to theBoltzmann distribution is very slow. Because of its simplicity and efficiency,this monomer-solvent coupling has been used in many polymer [14,71–74] andcolloid simulations [15, 16, 75, 76]. In order to accurately resolve the local flow field around a colloid, more ac-curate methods have been proposed which exclude fluid-particles from theinterior of the colloid and mimic slip [19, 77] or no-slip [78] boundary condi-tions. The latter procedure is similar to what is known in molecular dynamicsas a “thermal wall” boundary condition: fluid particles which hit the colloidparticle are given a new, random velocity drawn from the following probabil-ity distributions for the normal velocity component, v N , and the tangentialcomponent, v T , p N ( v N ) = ( mv N /k B T ) exp (cid:0) − mv N / k B T (cid:1) , with v N > ,p T ( v T ) = (cid:112) m/ πk B T exp (cid:0) − mv T / k B T (cid:1) . (74)These probability distributions are constructed so that the probability distri-bution for particles near the wall remains Maxwellian. The probability distri-bution, p T , for the tangential components of the velocity is Maxwellian, andboth positive and negative values are permitted. The normal component mustbe positive, since after scattering at the surface, the particle must move awayfrom the wall. The form of p N is a reflection of the fact that more particleswith large | v N | hit the wall per unit time than with small | v N | [78]. ulti-Particle Collision Dynamics This procedure models a no-slip boundary condition at the surface of thecolloid, and also thermostats the fluid at the boundaries. For many non-equilibrium flow conditions, this may not be sufficient, and it may also benecessary to thermostat the bulk fluid also (compare Sec. 2.3). It should alsobe noted that Eqs. (74) will be a good approximation only if the radius ofthe embedded objects is much larger than the mean free path λ . For smallerparticles, corrections are needed.If a particle hits the surface at time t in the interval between n∆t and( n + 1) ∆t , the correct way to proceed would be to give the particle its newvelocity and then have it stream the remaining time ( n + 1) ∆t − t . However,such detailed resolution is not necessary. It has been found [16] that goodresults are also obtained using the following simple stochastic procedure. Ifa particle is found to have penetrated the colloid during the streaming step,one simply moves it to the boundary and then stream a distance v new ∆t (cid:15) ,where (cid:15) is a uniformly distributed random number in the interval [0,1].Another subtlety is worth mentioning. If two colloid particles are veryclose, it can happen that a solvent particle could hit the second colloid afterscattering off the first, all in the interval ∆t . Naively, one might be tempted tosimply forbid this from happening or ignore it. However, this would lead to astrong depletion-like attractive force between the colloids [16]. This effect canbe greatly reduced by allowing multiple collisions in which one solvent particleis repeatedly scattered off the two colloids. In every collision, momentum istransferred to one of the colloids, which pushes the colloids further apart.In practice, even allowing for up to ten multiple collisions cannot completelycancel the depletion interaction—one needs an additional repulsive force toeliminate this unphysical attraction. The same effect can occur when a colloidparticle is near a wall.Careful tests of this thermal coupling have been performed by Padding etal. [17, 79], who were able to reproduce the correct rotational diffusion of acolloid. It should be noted that because the coupling between the solvent par-ticles and the surface occurs only through the movement of the fluid particles,the coupling is quite weak for small mean free paths. Another procedure for coupling an embedded object to the solvent has beenpursued by Kapral et al. [19, 30, 80]. They introduce a central repulsive forcebetween the solvent particles and the colloid. This force has to be quite strongin order to prohibit a large number of solvent particles from penetrating thecolloid. When implementing this procedure, a small time step δt is thereforerequired in order to resolve these forces correctly, and a large number of molec-ular dynamics time steps are needed during the SRD streaming step. In itsoriginal form, central forces were used, so that only slip boundary conditions could be modeled. In principle, non-central forces could be used to imposeno-slip conditions.This approach is quite natural and very easy to implement; it does, how-ever, require the use of small time steps and therefore may not be the optimalprocedure for many applications. One of the first approaches employed to impose a non-slip boundary con-dition at an external wall or at a moving object in a MPC solvent was touse “ghost” or “wall” particles [36, 81]. In other mesoscale methods such asLattice-Boltzmann, no-slip conditions are modeled using the bounce-back rule:the velocity of particle is inverted from v to − v when it intersects a wall. Forplanar walls which coincide with the boundaries of the collision cells, the sameprocedure can be used in MPC. However, the walls will generally not coincidewith, or even be parallel to, the cell walls. Furthermore, for small mean freepaths, where a shift of the cell lattice is required to guarantee Galilean invari-ance, partially occupied boundary cells are unavoidable, even in the simplestflow geometries. Fig. 6.
Velocity field of a fluid near a square cylinder in a Poiseuille flow at Reynoldsnumber Re = v max L/ν = 30. The channel width is eight times larger than thecylinder size L . A pair of stationary vortices are seen behind the obstacle, as expectedfor Re ≤
60. From Ref. [81].
The simple bounce-back rule fails to guarantee no-slip boundary conditionsin the case of partially filled cells. The following generalization of the bounce-back rule has therefore been suggested. For all cells that are cut by walls,fill the “wall” part of the cell with a sufficient number of virtual particles inorder to make the total number of particles equal to M , the average number ofparticles per cell. The velocities of the wall particles are drawn from a Maxwell-Boltzmann distribution with zero mean velocity and the same temperature as ulti-Particle Collision Dynamics the fluid. The collision step is then carried out using the mean velocity ofall particles in the cell. Note that since Gaussian random numbers are used,and the sum of Gaussian random numbers is also Gaussian-distributed, thevelocities of the individual wall particles need not to be determined explicitly.Instead, the average velocity u can be written as u = ( (cid:80) ni =1 v i + a ) /M , where a is vector whose components are Gaussian random numbers with zero meanand variance ( M − n ) k B T . Results for Poiseuille flow obtained using thisprocedure, both with and without cell shifting, were found to be in excellentagreement with the correct parabolic profile [36]. Similarly, numerical resultsfor the recirculation length, the drag coefficient, and the Strouhal number forflows around a circular and square cylinder in two dimensions were shownto be in good agreement with experimental results and computational fluiddynamics data for a range of Reynolds numbers between Re = 10 and Re =130 (see Fig. 6) [36, 81]. As an example of a situation in which it is important to use an algorithmwhich conserves angular momentum, consider a drop of a highly viscous fluidinside a lower-viscosity fluid in circular Couette flow. In order to avoid thecomplications of phase-separating two-component fluids, the high viscosityfluid is confined to a radius r < R by an impenetrable boundary with re-flecting boundary conditions (i.e., the momentum parallel to the boundary isconserved in collisions). No-slip boundary conditions between the inner andouter fluids are guaranteed because collision cells reach across the boundary.When a torque is applied to the outer circular wall (with no-slip, bounce-backboundary conditions) of radius R > R , a solid-body rotation of both fluids isexpected. The results of simulations with both MPC-AT − a and MPC-AT+ a are shown in Fig. 7. While MPC-AT+ a reproduces the expected behavior,MPC-AT − a produces different angular velocities in the two fluids, with a low(high) angular velocity in the fluid of high (low) viscosity [38].The origin of this behavior is that the viscous stress tensor in general hassymmetric and antisymmetric contributions (see Sec. 4.1.1), σ αβ = λ ( ∂ γ v γ ) δ αβ + ¯ η ( ∂ β v α + ∂ α v β ) + ˇ η ( ∂ β v α − ∂ α v β ) , (75)where λ is the second viscosity coefficient and ¯ η ≡ ρν and ˇ η ≡ ρν are thesymmetric and anti-symmetric components of the viscosity, respectively. Thelast term in Eq. (75) is linear in the vorticity ∇ × v , and does not conserveangular momentum. This term therefore vanishes ( i.e. , ˇ η = 0) when angularmomentum is conserved.The anti-symmetric part of the stress tensor implies an additional torque,which becomes relevant when the boundary condition is given by forces. Incylindrical coordinates ( r, θ, z ), the azimuthal stress is given by [38] Fig. 7.
Azimuthal velocity of binary fluids in a rotating cylinder with Ω =0 . k B T /m a ) / . The viscous fluids with particle mass m and m are locatedat r < R and R < r < R , respectively, with R = 5 a and R = 10 a . Symbolsrepresent the simulation results of MPC-AT − a with m /m = 2 (+) or m /m = 5( × ), and MPC-AT+ a for m /m = 5 ( ◦ ). Solid lines represent the analytical resultsfor MPC-AT − a at m /m = 5. Error bars are smaller than the size of the symbols.From Ref. [38]. σ rθ = (¯ η + ˇ η ) r∂ ( v θ /r ) ∂r + 2ˇ η v θ r . (76)The first term is the stress of the angular-momentum-conserving fluid, whichdepends on the derivative of the angular velocity Ω = v θ /r . The second termis an additional stress caused by the lack of angular momentum conservation;it is proportional to Ω .In the case of the phase-separated fluids in circular Couette flow, thisimplies that if both fluids rotate at the same angular velocity, the inner andouter stresses do not coincide. Thus, the angular velocity of the inner fluid Ω is smaller than the outer one, with v θ ( r ) = Ω r for r < R and v θ ( r ) = Ar + B/r, with A = Ω R − Ω R R − R , B = ( Ω − Ω ) R R R − R (77)for R < r < R . Ω is then obtained from the stress balance at r = R , i.e. , 2ˇ η Ω = (8 / η ( Ω − Ω ) + 2ˇ η Ω . This calculation reproduces thenumerical results very well, see Fig. 7. Thus, it is essential to employ an+ a version of MPC in simulations of multi-phase flows of binary fluids withdifferent viscosities.There are other situations in which the lack angular momentum conserva-tion can cause significant deviations. In Ref. [38], a star polymer with smallmonomer spacing was placed in the middle of a rotating Couette cell. As in theprevious case, it was observed that the polymer fluid rotated with a smaller ulti-Particle Collision Dynamics angular velocity than the outer fluid. When the angular momentum conser-vation was switched on, everything rotated at the same angular velocity, asexpected. The importance of hydrodynamic interactions (HI) in complex fluids is gener-ally accepted. A standard procedure for determining the influence of HI is toinvestigate the same system with and without HI. In order to compare results,however, the two simulations must differ as little as possible—apart from theinclusion of HI. A well-known example of this approach is Stokesian dynamicssimulations (SD), where the original Brownian dynamics (BD) method canbe extended by including hydrodynamic interactions in the mobility matrixby employing the Oseen tensor [12, 70].A method for switching off HI in MPC has been proposed in Refs. [24,26].The basic idea is to randomly interchange velocities of all solvent particlesafter each collision step, so that momentum (and energy) are not conserved locally . Hydrodynamic correlations are therefore destroyed, while leaving fric-tion coefficients and fluid self-diffusion coefficients largely unaffected. Sincethis approach requires the same numerical effort as the original MPC algo-rithm, a more efficient method has been suggested recently in Ref. [25]. If thevelocities of the solvent particles are not correlated, it is no longer necessaryto follow their trajectories. In a random solvent, the solvent-solute interac-tion in the collision step can thus be replaced by the interaction with a heatbath. This strategy is related to that proposed in Ref. [36] to model no-slipboundary conditions of solvent particles at a planar wall, compare Sec. 7.4.Since the positions of the solvent particles within a cell are not required inthe collision step, no explicit particles have to be considered. Instead, eachmonomer is coupled with an effective solvent momentum P which is directlychosen from a Maxwell-Boltzmann distribution of variance mM k B T and amean given by the average momentum of the fluid field—which is zero at rest,or ( mM ˙ γr iy , ,
0) in the case of an imposed shear flow. The total center-of-massvelocity, which is used in the collision step, is then given by [25] v cm,i = m m v i + P mM + m m , (78)where m m is the mass of the solute particle. The solute trajectory is thendetermined using MD, and the interaction with the solvent is performed everycollision time ∆t .The random MPC solvent therefore has similar properties to the MPCsolvent, except that there are no HI. The relevant parameters in both methodsare the average number of particles per cell, M , the rotation angle α , and thecollision time ∆t which can be chosen to be the same. For small values of the density ( M <
The relevance of hydrodynamic interactions for the dynamics of complexfluids—such as dilute or semidilute polymer solutions, colloid suspensions,and microemulsions—is well known [12, 70]. From the simulation point ofview, however, these systems are difficult to study because of the large gap inlength- and time-scales between solute and solvent dynamics. One possibilityfor investigating complex fluids is the straightforward application of molec-ular dynamics simulations (MD), in which the fluid is course-grained andrepresented by Lennard-Jones particles. Such simulations provide valuable in-sight into polymer dynamics [83–87]. Similarly, mesoscale algorithms such asLattice-Boltzmann and dissipative particle dynamics have been widely usedfor modeling of colloidal and polymeric systems [88–92].Solute molecules, e.g., polymers, are typically composed of a large numberof individual particles, whose interactions are described by a force-field. Asdiscussed in Sec. 7, the particle-based character of the MPC solvent allowsfor an easy and controlled coupling between the solvent and solute particles.Hybrid simulations combining MPC and molecular dynamics simulations aretherefore easy to implement. Results of such hybrid simulations are discussedin the following.
Many applications in chemical engineering, geology, and biology involve sys-tems of particles immersed in a liquid or gas flow. Examples include sedimen-tation processes, liquid-solid fluidized beds, and flocculation in suspensions.Long-ranged solvent-mediated hydrodynamic interactions have a profoundeffect on the non-equilibrium properties of colloidal suspensions, and themany-body hydrodynamic backflow effect makes it difficult to answer evenrelatively simple questions such what happens when a collection of parti-cles sediments through a viscous fluid. Batchelor [93] calculated the lowest-order volume fraction correction to the average sedimentation velocity, v s = v s (1 − . φ ), of hard spheres of hydrodynamic radius R H where v s is thesedimentation velocity of a single sphere. Due to the complicated interplay ulti-Particle Collision Dynamics between short-ranged contact forces and long-ranged HI, it is hard to extendthis result to the high volume fraction suspensions of interest for ceramicsand soil mechanics. An additional complication is that the Brownian motionof solute particles in water cannot be neglected if they are smaller than 1 µm in diameter.The dimensionless Peclet number P e = v s R H /D , where D is the self-diffusion coefficient of the suspended particles, measures the relative strengthof HI and thermal motion. Most studies of sedimentation have focused onthe limit of infinite Peclet number, where Brownian forces are negligible. Forexample, Ladd [94] employed a Lattice-Boltzmann method (LB), and Hoeflerand Schwarzer [95] used a marker-and-cell Navier-Stokes solver to simulatesuch non-Brownian suspensions. The main difficulty with such algorithms isthe solid-fluid coupling which can be very tricky: in LB simulations, special“boundary nodes” were inserted on the colloid surface, while in Ref. [95], thecoupling was mediated by inertia-less markers which are connected to thecolloid by stiff springs and swim in the fluid, effectively dragging the colloid,but also exerting a force on the fluid. Several methods for coupling embeddedparticles to an MPC solvent were discussed in Sec. 7.Using the force-based solvent-colloid coupling described in Sec. 7.3, Paddingand Louis [96] investigated the importance of HI during sedimentation at smallPeclet numbers. Surprisingly, they found that the sedimentation velocity doesnot change if the Peclet number is varied between 0.1 and 15 for a range of vol-ume fractions. For small volume fractions, the numerical results agree with theBatchelor law; for intermediate φ they are consistent with the semi-empiricalRichardson-Zaki law, v s = v s (1 − φ ) n , n = 6 .
55. Even better agreement wasfound with theoretical predictions by Hayakawa and Ichiki [97, 98], who tookhigher order HI into account. Purely hydrodynamic arguments are thereforestill valid in an average sense at low
P e , i.e. , for strong Brownian motion andrelatively weak HI. This also means that pure Brownian simulations withoutHI, which lead to v s = v s (1 − φ ), strongly underestimate the effect of backflow.On the other hand, it is known that the velocity autocorrelation func-tion of a colloidal particle embedded in a fluctuating liquid at equilibriumexhibits a hydrodynamic long-time tail, (cid:104) v ( t ) v (0) (cid:105) ∼ t − d/ , where d is thespatial dimension [99]. These tails have been measured earlier for point-likeSRD particles in two [21, 27] and three [15] spatial dimensions, and found tobe in quantitative agreement with analytic predictions, with no adjustableparameters. It is therefore not surprising that good agreement was also ob-tained for embedded colloids [96]. MPC therefore correctly describes two ofthe most important effects in colloidal suspensions, thermal fluctuations andhydrodynamic interactions.In a series of papers, Hecht et al. [16, 76, 100] used hybrid SRD-moleculardynamics simulations to investigate a technologically important colloidalsystem— Al O -particles of diameter 0 . µm (which is often used in ceram-ics) suspended in water—with additional colloid-colloid interactions. Thesecolloids usually carry a charge which, by forming an electric double layer with ions in water, results in a screened electrostatic repulsion. The interaction canbe approximated by the Derjaguin-Landau-Verwey-Overbeek (DLVO) the-ory [101, 102]. The resulting potential contains a repulsive Debye-H¨uckel con-tributions, V EL ∼ exp( − κ [ r − d ]) /r , where d is the particle diameter, κ isthe inverse screening length, and r is the distance of the particle centers. Thesecond part of the DLVO-potential is a short-range van der Waals attraction, V vdW = − A H (cid:20) d r − d + d r + 2ln (cid:18) r − d r (cid:19)(cid:21) , (79)which turns out to be important at the high volume fractions ( φ > A H is the Hamaker constant which in-volves the polarizability of the particles and the solvent. DLVO theory makesthe assumption of linear polarizability and is valid only at larger distances. Ittherefore does not include the so-called primary potential minimum at particlecontact, which is observed experimentally and is about 30 k B T deep. Becauseof this potential minimum, colloids which come in contact rarely become freeagain. In order to ensure numerical stability for reasonable values of the timestep, this minimum was modeled by an additional parabolic potential withdepth of order 6 k B T . The particle Reynolds number of the real system isvery small, of order 10 − to 10 − . Since it would be too time-consuming tomodel this Reynolds number, the simulations were performed at Re ≈ . V Hertz ∼ ( d − r ) / if r < d . (80)SRD correctly describes long-range HI, but it can only resolve hydrody-namic interactions on scales larger than both the mean free path λ and the cellsize a . In a typical simulation with about 1000 colloid particles, a relativelysmall colloid diameter of about four lattice units was chosen for computationalefficiency. This means that HI are not fully resolved at interparticle distancescomparable to the colloid diameter, and lubrication forces have to be insertedby hand. Only the most divergent mode, the so-called squeezing mode, wasused, F lub ∼ v rel /r rel , where r rel and v rel are the relative distance and ve-locity of two colloids, respectively. This system of interacting Al O -particleswas simulated in order to study the dependence of the suspension’s viscosityand structure on shear rate, pH, ionic strength and volume fraction. The re-sulting stability diagram of the suspension as a function of ionic strength andpH value is shown in Fig. 8 (plotted at zero shear) [76]. The pH controls thesurface charge density which, in turn, affects the electrostatic interactions be-tween the colloids. Increasing the ionic strength, experimentally achieved by“adding salt”, decreases the screening length 1 /κ , so that the attractive forces ulti-Particle Collision Dynamics become more important; the particles start forming clusters. Three differentstates are observed: (i) a clustered regime, where particles aggregate when vander Waals attractions dominate, (ii) a suspended regime where particles aredistributed homogeneously and can move freely—corresponding to a stablesuspension favored when electrostatic repulsion prevents clustering but is notstrong enough to induce order. At very strong Coulomb repulsion the repul-sive regime (iii) occurs, where the mobility of the particles is restricted, andparticles arrange in local order which maximizes nearest neighbor distances.The location of the phase boundaries in Fig. 8 depends on the shear rate. Inthe clustered phase, shear leads to a breakup of clusters, and for the shear rate˙ γ = 1000 s − , there are many small clusters which behave like single particles.In the regime where the particles are slightly clustered, or suspended, shearthinning is observed. Shear thinning is more pronounced in the slightly clus-tered state, because shear tends to reduce cluster size. Reasonable agreementwith experiments was achieved, and discrepancies were attributed to poly-dispersity and the manner in which lubrication forces were approximated, aswell as uncertainties how the pH and ionic strength enter the model forceparameters. Fig. 8.
Phase diagram of a colloidal suspension (plotted at zero shear and volumefraction φ = 35%) in the ionic-strength-pH plane depicting three regions: a clusteredregion, a suspended regime, and a repulsive structure. From Ref. [76].6 G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler In the simulations of Hecht et al. [16], the simple collisional coupling proce-dure described in Sec. 7.1 was used. This means that the colloids were treatedas point particles, and solvent particles could flow right through them. Hydro-dynamic interactions were therefore only resolved in an average sense, whichis acceptable for studies of the general properties of an ensemble of manycolloids. The heat from viscous heating was removed using the stochasticthermostat described in Sec. 2.3.Various methods for modeling no-slip boundary conditions at colloidsurfaces—such as the thermal wall coupling described in Sec. 7.2—were sys-tematically investigated in Ref. [79]. No-slip boundary conditions are impor-tant, since colloids are typically not completely spherical or smooth, and thesolvent molecules also transfer angular momentum to the colloid. Using theSRD algorithm without angular momentum conservation, it was found thatthe rotational friction coefficient was larger than predicted by Enskog-theorywhen the ghost-particle coupling was used [82]. On the other hand, in a de-tailed study of the translational and rotational velocity autocorrelation func-tion of a sphere coupled to the solvent by the thermal-wall boundary condition,quantitative agreement with Enskog theory was observed at short times, andwith mode-coupling theory at long times. However, it was also noticed thatfor small particles, the Enskog and hydrodynamic contributions to the frictioncoefficients were not clearly separated. Specifically, mapping the system to adensity matched colloid in water, it appeared that the Enskog and the hydro-dynamic contributions are equal at a particle radius of 6 nm for translationand 35 . nm for rotation; even for a particle radius of 100 nm , the Enskogcontribution to the friction is still of order 30% and cannot be ignored.In order to clarify the detailed character of the hydrodynamic interactionsbetween colloids in SRD, Lee and Kapral [103] numerically evaluated the fixed-particle friction tensor for two nano-spheres embedded in an SRD solvent.They found that for intercolloidal spacings less than 1 . d , where d is thecolloid diameter, the measured friction coefficients start to deviate from theexpected theoretical curve. The reader is referred to the review by Kapral [30]for more details. The dynamical behavior of macromolecules in solution is strongly affected oreven dominated by hydrodynamic interactions [70, 104, 105]. From a theoret-ical point of view, scaling relations predicted by the Zimm model for, e.g.,the dependencies of dynamical quantities on the length of the polymer are, ingeneral, accepted and confirmed [106]. Recent advances in experimental single-molecule techniques provide insight into the dynamics of individual polymers,and raise the need for a quantitative theoretical description in order to deter-mine molecular parameters such as diffusion coefficients and relaxation times.Mesoscale hydrodynamic simulations can be used to verify the validity of the- ulti-Particle Collision Dynamics oretical models. Even more, such simulations are especially valuable whenanalytical methods fail, as for more complicated molecules such as polymerbrushes, stars, ultrasoft colloids, or semidilute solutions, where hydrodynamicinteractions are screened to a certain degree. Here, mesoscale simulations stillprovide a full characterization of the polymer dynamics.We will focus on the dynamics of polymer chains in dilute solution. In orderto compare simulation results with theory—in particular the Zimm approach[70, 107]—and scaling predictions, we address the dynamics of Gaussian aswell as self-avoiding polymers. Polymer molecules are composed of a large number of equal repeat units calledmonomers. To account for the generic features of polymers, such as their con-formational freedom, no detailed modeling of the basic units is necessary. Acoarse-grained description often suffices, where several monomers are com-prised in an effective particle. Adopting such an approach, a polymer chainis introduced into the MPC solvent by adding N m point particles, each ofmass m m , which are connected linearly by bonds. Two different models areconsidered, a Gaussian polymer and a polymer with excluded-volume (EV)interactions. Correspondingly, the following potentials are applied:(i) Gaussian chain : The monomers, with the positions r i ( i = 1 , . . . , N m ), areconnected by the harmonic potential U G = 3 k B T b N m − (cid:88) i =1 ( r i +1 − r i ) , (81)with zero mean bond length, and b the root-mean-square bond length. Here,the various monomers freely penetrated each other. This simplification allowsfor an analytical treatment of the chain dynamics as in the Zimm model[70, 107].(ii) Excluded-volume chain : The monomers are connected by the harmonicpotential U B = κ N m − (cid:88) i =1 ( | r i +1 − r i | − b ) , (82)with mean bond length b . The force constant κ is chosen such that the fluc-tuations of the bond lengths are on the order of a percent of the mean bondlength. In addition, non-bonded monomers interact via the repulsive, trun-cated Lennard-Jones potential U LJ = (cid:40) (cid:15) (cid:104)(cid:0) σr (cid:1) − (cid:0) σr (cid:1) (cid:105) + (cid:15), r < / σ , otherwise . (83) The excluded volume leads to swelling of the polymer structure compared to aGaussian chain, which is difficult to fully account for in analytical calculations[73].The dynamics of the chain monomers is determined by Newtons’ equationsof motion between the collisions with the solvent. These equations are inte-grated using the velocity Verlet algorithm with the time step ∆t p . The latteris typically smaller than the collision time ∆t . The monomer-solvent interac-tion is taken into account by inclusion of the monomer of mass m m = ρm in the collision step [68, 73], compare Sec. 7.1. Alternatively, a Lennard-Jonespotential can be used to account for the monomer-MPC particle interaction,where a MPC particle is of zero interaction range [19, 108].We scale length and time according to ˆ x = x/a and ˆ t = t (cid:112) k B T /ma ,which corresponds to the choice k B T = 1, m = 1, and a = 1. The mean freepath of a fluid particle ∆t (cid:112) k B T /m is then given by λ = ∆ ˆ t . In addition, weset b = a , σ = a , and (cid:15)/k B T = 1.The equilibrium properties of a polymer are not affected by hydrodynamicinteractions. Indeed, the results for various equilibrium quantities—such asthe radius of gyration—of MPC simulation are in excellent agreement withthe results of molecular dynamics of Monte Carlo simulations without explicitsolvent [73].Simulations of Gaussian chains, i.e., polymers with the bond potential(81), can be compared with analytical calculations based on the Zimm ap-proach [70,107]. Note, however, that the simulations are not performed in theZimm model. The Zimm approach relies on the preaveraging approximationof hydrodynamic interactions, whereas the simulations take into account theconfigurational dependence of the hydrodynamic interactions, and thereforehydrodynamic fluctuations. Hence, the comparison can serve as a test of thevalidity of the approximations employed in the Zimm approach.The Zimm model rests upon the Langevin equation for over-damped mo-tion of the monomers, i.e., it applies for times larger than the Brownian timescale τ B (cid:29) m m /ζ , where ζ is Stokes’ friction coefficient [12]. On such timescales, velocity correlation functions have decayed to zero and the monomermomenta are in equilibrium with the solvent. Moreover, hydrodynamic in-teractions between the various parts of the polymer are assumed to propa-gate instantaneously. This is not the case in our simulations. First of all, themonomer inertia term is taken into account, which implies non-zero velocityautocorrelation functions. Secondly, the hydrodynamic interactions build upgradually. The center-of-mass velocity autocorrelation function displayed inFig. 9 reflects these aspects. The correlation function exhibits a long-time tail,which decays as (cid:104) v cm ( t ) v cm (0) (cid:105) ∼ t − / on larger time scales. The algebraicdecay is associated with a coupling between the motion of the polymer andthe hydrodynamic modes of the fluid [99, 109, 110]. A scaling of time withthe diffusion coefficient D shows that the correlation function is a universalfunction of Dt . This is in agreement with results of DPD simulations of dilutepolymer systems [92]. ulti-Particle Collision Dynamics Fig. 9.
Center-of-mass velocity autocorrelation functions for Gaussian polymers oflength N m = 20, N m = 40, and λ = 0 . Dt . The solid line isproportional to ( Dt ) − / . From Ref. [73]. The polymer center-of-mass diffusion coefficient follows either via theGreen-Kubo relation from the velocity autocorrelation function or by the Ein-stein relation from the center-of-mass mean square displacement. Accordingto the Kirkwood formula [104, 105, 111] D ( K ) = D N m + k B T πη R H , (84)where the hydrodynamic radius R H is defined as1 R H = 1 N m (cid:42) N m (cid:88) i =1 N m (cid:88) j =1 (cid:48) | r i − r j | (cid:43) (85)and the prime indicates that the term with j = i has to be left out in thesummation. The diffusion coefficient is composed of the local friction contri-bution D /N m , where D is the diffusion coefficient of a single monomer inthe same solvent, and the hydrodynamic contribution.Simulation results for the hydrodynamic contribution, D H = D − D /N m ,to the diffusion coefficient are plotted in Fig. 10 as a function of the hydrody-namic radius (85). In the limit N m (cid:29)
1, the diffusion coefficient D is domi-nated by the hydrodynamic contribution D H , since D H ∼ N − / m . For shorterchains, D /N m cannot be neglected, and therefore has to be subtracted inorder to extract the scaling behavior of D H . The hydrodynamic part of thediffusion coefficient D H exhibits the dependence predicted by the Kirkwoodformula and the Zimm theory, i.e., D H ∼ /R H . The finite-size corrections to D show a dependence D = D ∞ − const ./L on the size L of a periodic system, Fig. 10.
Dependence of the hydrodynamic part of the diffusion coefficient, D H = D − D /N m , on the hydrodynamic radius for Gaussian chains of lengths N m = 5,10, 20, 40, 80, and 160 (right to left). The mean free path is λ = 0 .
1. From Ref. [73]. in agreement with previous studies [68, 91, 112]. Simulations for various sys-tem sizes for polymers of lengths N m = 10, 20, and 40 allow an extrapolationto infinite system size, which yields D / (cid:112) k B T a /m ≈ . × − , in goodagreement with the diffusion coefficient of a monomer in the same solvent.The values of D ∞ are about 30% larger than the finite-system-size values pre-sented in Fig. 10. Similarly the diffusion coefficient for a polymer chain withexcluded volume interactions displays the dependence D H ∼ /R H [73].The Kirkwood formula neglects hydrodynamic fluctuations and is thusidentical with the preaveraging result of the Zimm approach. When only thehydrodynamic part is considered, the Zimm model yields the diffusion coeffi-cient D Z = 0 . k B Tbη √ N m . (86)MPC simulations for polymers of length N m = 40 yield D Z / (cid:112) k B T a /m =0 . D H / (cid:112) k B T a /m = 0 . D ( K ) , in agreement with previous studies pre-sented in Refs. [70,111,113]. Note that the experimental values are also smallerby about 15% than those predicted by the Zimm approach [70, 114, 115].To further characterize the internal dynamics of the molecular chain, amode analysis in terms of the eigenfunctions of the discrete Rouse model [70,116] has been performed. The mode amplitudes χ p are calculated accordingto ulti-Particle Collision Dynamics Fig. 11.
Correlation functions of the Rouse-mode amplitudes for the modes p = 1 − N m = 20 (right) and N m = 40 (left).From Ref. [73]. χ p = (cid:114) N m N m (cid:88) i =1 r i cos (cid:20) pπN m (cid:18) i − (cid:19)(cid:21) , p = 1 . . . N m . (87)Due to hydrodynamic interactions, Rouse modes are no longer eigenfunctionsof the chain molecule. However, within the Zimm theory, they are reasonableapproximations and the autocorrelation functions of the mode amplitudesdecay exponentially, i.e., (cid:10) χ p ( t ) χ p (0) (cid:11) = (cid:10) χ p (cid:11) exp ( − t/τ p ) . (88)For the Rouse model, the relaxation times τ p depend on chain length and modenumber according to τ p ∼ / sin ( pπ/ N m ), whereas for the Zimm model thedependence τ p ∼ ( p/N m ) / / sin ( pπ/ N m ) (89)is obtained. The extra contribution (cid:112) p/N m follows from the eigenfunctionrepresentation of the preaveraged hydrodynamic tensor, under the assumptionthat its off-diagonal elements do not significantly contribute to the relaxationbehavior.In Fig. 11, the autocorrelation functions for the mode amplitudes areshown for the mean free path λ = 0 .
1. Within the accuracy of the simula-tions, the correlation functions decay exponentially and exhibit the scalingbehavior predicted by the Zimm model. Hence, for the small mean free path,hydrodynamic interactions are taken into account correctly. This is no longer
Fig. 12.
Dependence of the longest relaxation time τ on the radius of gyration forGaussian chains of the lengths given in Fig. 10. From Ref. [73]. true for the large mean free path, λ = 2. In this case, a scaling behavior be-tween that predicted by the Rouse and Zimm models is observed. This impliesthat hydrodynamic interactions are present, but are not fully developed or aresmall compared to the local friction of the monomers. We obtain pure Rousebehavior for a system without solvent by simply rotating the velocities of theindividual monomers [73].The dependence of the longest relaxation time on the radius of gyrationis displayed in Fig. 12 for λ = 0 .
1. The scaling behavior τ ∼ R g is in verygood agreement with the predictions of the Zimm theory. We even find al-most quantitative agreement; the relaxation time of the p = 1 mode of oursimulations is approximately 30 % larger than the Zimm value [70].The scaling behavior of equilibrium properties of single polymers withexcluded-volume interactions has been studied extensively [70, 117–120]. Ithas been found that even very short chains already follow the scaling be-havior expected for much longer chains. In particular, the radius of gyrationincreases like R G ∼ N νm with the number of monomers, and the static struc-ture factor S ( q ) exhibits a scaling regime for 2 π/R G (cid:28) q (cid:28) π/σ , with a q − /ν decay as a function of the scattering vector q and the exponent ν ≈ . b = σ = a , (cid:15)/k B T = 1, the exponent ν ≈ .
62 is obtained from the chain-length depen-dence of the radius of gyration, the mean square end-to-end distance, as wellas the q − dependence of the static structure factor [73].An analysis of the intramolecular dynamics in terms of the Rouse modesyields non-exponentially decaying autocorrelation functions of the mode am- ulti-Particle Collision Dynamics Fig. 13.
Correlation functions of the Rouse-mode amplitudes for various modes as afunction of the scaled time tp α for polymers with excluded volume interactions. Thechain lengths are N m = 20 (left) and N m = 40 (right). The calculated correlationswhere fitted by A p exp( − t/τ p ) and have been divided by A p . The scaling exponentsof the mode numbers are α = 1 .
93 ( N m = 20) and α = 1 .
85 ( N m = 40), respectively.From Ref. [73]. plitudes. At very short times, a fast decay is found, which turns into a slowerexponential decay which is well fitted by A p exp( − t/τ p ), see Fig. 13. Withinthe accuracy of these calculations, the correlation functions exhibit universalbehavior. Zimm theory predicts the dependence τ p ∼ p − ν for the relax-ation times on the mode number for polymers with excluded-volume interac-tions [70]. With ν = 0 .
62, the exponent α for the polymer of length N m = 40is found to be in excellent agreement with the theoretical prediction. Theexponent for the polymers with N m = 20 is slightly larger.Zimm theory predicts that the dynamic structure factor, which is definedby S ( q , t ) = 1 N m N m (cid:88) i =1 N m (cid:88) j =1 (cid:104) exp ( i q [ r i ( t ) − r j (0)]) (cid:105) , (90)scales as [70] S ( q , t ) = S ( q , f ( q α t ) (91)with α = 3 for qR G (cid:29)
1, independent of the solvent conditions ( Θ or good sol-vent). To extract the scaling relation for the intramolecular dynamics, whichcorresponds to the prediction (91), we resort to the following considerations.As is well known, the dynamic structure factor for a Gaussian distribution of Fig. 14.
Normalized dynamic structure factor S ( q , t ) / ( S ( q ,
0) exp ` − Dq t ´ ) of poly-mers with excluded volume interactions for N m = 40 and various q -values in therange 0 . < qa < q t / . From Ref. [73]. the differences r i ( t ) − r j (0) and a linear equation of motion is given by [70,121] S ( q , t ) = S ( q , e − Dq t N m N m (cid:88) i =1 N m (cid:88) j =1 exp (cid:0) − q (cid:10) ( r (cid:48) i ( t ) − r (cid:48) j (0)) (cid:11) / (cid:1) , (92)where Dq t accounts for the center-of-mass dynamics and r (cid:48) i denotes theposition of monomer i in the center-of-mass reference frame. Therefore,in order to obtain the dynamics in the center-of-mass reference frame, weplot S ( q , t ) / ( S ( q ,
0) exp (cid:0) − Dq t (cid:1) ). The simulation results for the polymer oflength N m = 40, shown in Fig. 14, confirm the predicted scaling behavior.Thus, MPC-MD hybrid simulations are very well suited to study the dynam-ics of even short polymers in dilute solution.As mentioned above, the structure of a polymer depends on the nature ofthe solvent. In good solvent excluded volume interactions lead to expandedconformations and under bad solvent conditions the polymer forms a densecoil. In a number of simulations the influence of hydrodynamic interactions onthe transition from an extended to a collapsed state has been studied, whenthe solvent quality is abruptly changed. Both, molecular dynamics simula-tions with an explicit solvent [122] as well as MPC simulations [108, 123, 124]yield a significant different dynamics in the presence of hydrodynamic inter-actions. Specifically, the collapse is faster, with a much weaker dependence ofthe characteristic time on polymer length [108, 124], and the folding path isaltered.Similarly, a strong influence of hydrodynamic interactions has been foundon the polymer translocation dynamics through a small hole in a wall [125] or ulti-Particle Collision Dynamics in polymer packing in a virus capsid [126, 127]. Cooperative backflow effectslead to a rather sharp distribution of translocation times with a peak at rela-tively short times. The fluid flow field, which is created as a monomer movesthrough the hole, guides following monomers to move in the same direction. Simulations of a MPC fluid confined between surfaces and exposed to a con-stant external force yield the expected parabolic velocity profile for appropri-ate boundary conditions [31, 36, 128]. The ability of MPC to account for theflow behavior of mesoscale objects, such as polymers, under non-equilibriumconditions has been demonstrate for a number of systems. Rod-like colloidsin shear flow exhibit flow induced alignment [72]. The various diagonal com-ponents of the radius of gyration tensor exhibit a qualitative and quantitativedifferent behavior. Due to the orientation, the component in the flow direc-tion increases with increasing Peclet number larger than united and saturatesat large shear rates because of finite size effects. The transverse componentsdecrease with shear rate, where the component in the gradient direction isreduced to a greater extent. The rod rotational velocity in the shear planeshows two distinct regimes. For Peclet numbers much smaller than unity, therotational velocity increases linearly with the shear rate, because the systemis isotropic. At Peclet numbers much larger than unity, the shear inducedanisotropies lead to a slower increase of the rotational velocity with the shearrate [72].The simulations of a tethered polymer in a Poiseuille flow [74] yield a seriesof morphological transitions from sphere to deformed sphere to trumpet tostem and flower to rod, similar to theoretically predicted structures [129–131]. The crossovers between the various regimes occur at flow rates close tothe theoretical estimates for a similar system. Moreover, the simulations inRef. [74] show that backflow effects lead to an effective increase in viscosity,which is attributed to the fluctuations of the free polymer end rather than itsshape.The conformational, structural, and transport properties of free flexiblepolymers in microchannel flow have been studied in Refs. [128,132] by hybridMPC-molecular dynamics simulations. These simulations confirm the crossstreamline migration of the molecules as previously observed in Refs. [133–138]. In addition, various other polymer properties are addressed in Ref. [132].All these hybrid simulations confirm that MPC is an excellent method tostudy the non-equilibrium behavior of polymers in flow fields. In the next sec-tion, we will provide a more detailed example for a more complicated object,namely an ultrasoft colloid in shear flow.
Fig. 15.
Fluid flow lines in the flow-gradient plane of the star polymer’s center ofmass reference frame for f = 10 (top) and f = 50 (bottom) arms, both with L f = 20monomers per arm and an applied shear field with Wi = ˙ γτ = 22. From Ref. [25]. Star polymers present a special macromolecular architecture, in which sev-eral linear polymers of identical length are linked together by one of theirends at a common center. This structure is particularly interesting because itallows for an almost continuous change of properties from that of a flexiblelinear polymer to a spherical colloidal particle with very soft interactions. Starpolymers are therefore also often called ultrasoft colloids. The properties ofstar polymers in and close to equilibrium have been studied intensively, boththeoretically [139,140] and experimentally [141]. A star polymer is a ultrasoft ulti-Particle Collision Dynamics Fig. 16.
Normalized component G xx of the average gyration tensor as a functionof the Weissenberg number Wi, for star polymers of f = 50 arms and the armlengths L f = 10, 20, and 40 monomers. Power-law behaviors with quadratic andlinear dependencies on Wi are indicated by lines. The inset shows all three diagonalcomponents of the gyration tensor (for L f = 20). From Ref. [142]. colloid, where the core extension is very small compared to the length of anarm. By anchoring polymers on the surface of a hard colloid, the softnesscan continuously be changed from ultrasoft to hard by increasing the ratiobetween the core and shell radius at the expense of the thickness of the softpolymer corona. Moreover, star polymers have certain features in commonwith vesicles and droplets. Although their shell can be softer than that of theother objects, the dense packing of the monomers will lead to a cooperativedynamical behavior resembling that of vesicles or droplets [142].Vesicles and droplets encompass fluid which is not exchanged with thesurrounding. In contrast, for star-like molecules fluid is free to penetrate intothe molecule and internal fluid is exchanged with the surrounding in the courseof time. This intimate coupling of the star-polymer dynamics and the fluidflow leads to a strong modification of the flow behavior at and next to theultrasoft colloid particular in non-equilibrium systems.In the following, we will discuss a few aspects in the behavior of starpolymers in shear flow as a function of arm number f , arm length L f , andshear rate ˙ γ . The polymer model is the same as described in Sec. 10.2.1,where the chain connectivity is determined by the bond potential (82) andthe excluded-volume interaction is described by the Lennard-Jones potential(83). A star polymer of functionality f is modeled as f linear polymer chainsof L f monomers each, with one of their ends linked to a central particle.Linear polymer molecules are a special case of star polymers with function-ality f = 2. Lees-Edwards boundary conditions [42] are employed in order toimpose a linear velocity profile ( v x , v y , v z ) = ( ˙ γr y , ,
0) in the fluid in the ab- sence of a polymer. For small shear rates, the conformations of star polymersremain essentially unchanged compared to the equilibrium state. Only whenthe shear rate exceeds a characteristic value, a structural anisotropy as wellas an alignment is induced by the flow (cf. Fig. 15). The shear rate depen-dent quantities are typically presented in terms of the Weissenberg numberWi = ˙ γτ rather than the shear rate itself, where τ is the longest characteristicrelaxation time of the considered system. For the star polymers, the best datacollapse for stars of various arm lengths is found when the relaxation time τ = ηb L f /k B T is used [142]. Remarkably, there is essentially no dependenceon the functionality. Within the range of investigated star sizes, this relax-ation time has to be considered as consistent with the prediction for the blobmodel of Ref. [143], where τ ∼ L . f f . for the Flory exponent ν = 0 . stagnation points of vanishing fluid velocity [25].The fluid flow in the vicinity of the star polymer is distinctively differentfrom that of a sphere but resembles the flow around an ellipsoid [144]. Incontrast to the latter, the fluid penetrates into the area covered by the starpolymer. While the fluid in the core of the star rotates together with thepolymer, the fluid in the corona follows the external flow to a certain extent.A convenient quantity to characterize the structural properties and align-ment of polymers in flow is the average gyration tensor, which is defined as G αβ ( ˙ γ ) = 1 N m N m (cid:88) i =1 (cid:104) r (cid:48) i,α r (cid:48) i,β (cid:105) , (93)where N m = f L f + 1 is the total number of monomers, r (cid:48) i,α is the position ofmonomer i relative to the polymer center of mass, and α ∈ { x, y, z } . The av-erage gyration tensor is directly accessible in scattering experiments. Its diag-onal components G αα ( ˙ γ ) are the squared radii of gyration of the star polymeralong the axes of the reference frame. In the absence of flow, scaling consid-erations predict [139] G xx (0) = G yy (0) = G zz (0) = R g (0) / ∼ L νf f − ν .The diagonal components G αα of the average gyration tensor are shownin Fig. 16 as a function of the Weissenberg number for various functionalitiesand arm lengths. We find that the extension of a star increases with increasingshear rate in the shear direction ( x ), decreases in the gradient direction ( y ),and is almost independent of Wi in the vorticity direction ( z ). The deviationfrom spherical symmetry exhibits a Wi power-law dependence for small shear ulti-Particle Collision Dynamics Fig. 17.
Orientational resistance m G as a function of the Weissenberg number Wifor star polymers with functionalities f = 2 (circles), f = 15 (triangles), and f = 50(squares), and different arm lengths indicated in the figure. Lines correspond to thepower law m G = ˜ m G ( f )Wi . . The inset shows that the amplitude also follows apower-law behavior with ˜ m G ( f ) ∼ f . . From Ref. [142]. rates for all functionalities. A similar behavior has been found for rod-likecolloids [72] (due to the increasing alignment with the flow direction) and forlinear polymers [70]. However, for stars of not too small functionality, a newscaling regime appears, where the deformation seems to scales linearly withthe Weissenberg number. For large shear rates, finite-size effects appear dueto the finite monomer number. These effects emerge when the arms are nearlystretched, and therefore occur at higher Weissenberg numbers for larger armlengths.The average flow alignment of a (star) polymer can be characterized bythe orientation angle χ G , which is the angle between the eigenvector of thegyration tensor with the largest eigenvalue and the flow direction. It followsstraightforwardly [87] from the simulation data viatan(2 χ G ) = 2 G xy / ( G xx − G yy ) ≡ m G / Wi , (94)where the right-hand-side of the equation defines the orientation resistanceparameter m G [145]. It has been shown for several systems including rod-likecolloids and linear polymers without self-avoidance [70] that close to equilib-rium G xy ∼ ˙ γ and ( G xx − G yy ) ∼ ˙ γ , so that m G is independent of Wi. Ourresults for the orientation resistance are presented in Fig. 17 for various func-tionalities f and arm lengths L f . Data for different L f collapse onto universalcurves, which approach a plateau for small shear rates, as expected. For largershear rates, Wi (cid:29)
1, a power-law behavior [142] m G (Wi) ∼ f α Wi µ , (95) is obtained with respect to the Weissenberg number and the functionality,where α = 0 . ± .
02 and µ = 0 . ± .
05. For self-avoiding linear polymers,a somewhat smaller exponent µ = 0 . ± .
03 was obtained in Refs. [87,146], whereas theoretical calculations predict ∼ Wi / in the limit of largeWeissenberg numbers [147].The data for the average orientation and deformation of a star polymerdescribed so far seem to indicate that the properties vary smoothly andmonotonically from linear polymers to star polymers of high functionality.However, this picture changes when the dynamical behavior is considered.It is well known by now that linear polymers show a tumbling motion inflow, with alternating collapsed and stretched configurations during each cy-cle [146, 148–150]. For large Weissenberg numbers, this leads to very largefluctuations of the largest intramolecular distance of a linear polymer withtime, as demonstrated experimentally in Refs. [146, 150], and reproduced inthe MPC simulations [142], see Fig. 18. A similar behavior is found for f = 3.However, for f >
5, a quantitatively different behavior is observed as displayedin Fig. 19. Now, the fluctuations of the largest intramolecular distance aremuch smaller and decrease with increasing Weissenberg number as shown inFig. 18, and the dynamics resembles much more the continuous tank-treadingmotion of fluid droplets and capsules. The shape and orientation of such starsdepends very little on time, while the whole object is rotating. On the otherhand, a single, selected arm resembles qualitatively the behavior of a linearpolymer—it also collapses and stretches during the tank-treading motion. Thesuccessive snapshots of Fig. 20 illustrate the tank-treating motion. Followingthe top left red polymer in the top left image, we see that the extended poly-mer collapses in the course of time, moves to the right and stretches again. Inparallel, other polymers exhibit a similar behavior on the bottom side. More-over, the images show that the orientation of a star hardly changes in thecourse of time.The rotational dynamics of a star polymer can be characterized quantita-tively by calculating the rotation frequency ω α = z (cid:88) β = x (cid:104) Θ − αβ L β (cid:105) (96)of a star, where Θ αβ = N m (cid:88) i =1 [ r (cid:48) i δ αβ − r (cid:48) i,α r (cid:48) i,β ] (97)is the instantaneous moment-of-inertia tensor and L β is the instantaneous an-gular momentum. Since the rotation frequency for all kinds of soft objects—such as rods, linear polymers, droplets and capsules—depends linearly on ˙ γ for small shear rates, the reduced rotation frequency ω/ ˙ γ is shown in Fig. 21as a function of the Weissenberg number. The data approach ω/ ˙ γ = 1 / ulti-Particle Collision Dynamics Fig. 18.
Temporal evolution of the largest intramolecular distance e = max ij | r i − r j | of a linear polymer and a star polymer with 50 arms, for the Weissenberg numbersWi = 0 .
64 and Wi = 64. In both cases, L f = 40. The time t is measured in units ofthe collision time δt . e s corresponds to the fully stretched arms. From Ref. [142]. Fig. 19.
Widths of the distribution functions of the largest intramolecular distances, σ e = ( (cid:104) e (cid:105) − (cid:104) e (cid:105) ) / (cid:104) e (cid:105) , of a linear polymer and star polymers with up to 50 armsas a function of the Weissenberg number. From Ref. [142]. for small Wi, as expected [87, 151]. For larger shear rates, the reduced fre-quency decreases due to the deformation and alignment of the polymers inthe flow field. With increasing arm number, the decrease of ω/ ˙ γ at a givenWeissenberg number becomes smaller, since the deviation from the sphericalshape decreases. Remarkably, the frequency curves for all stars with f > ω/ ˙ γ is plottedas a function of a rescaled Weissenberg number, see Fig. 21. For high shear Fig. 20.
Successive snapshots of a star polymer of functionality f = 25 and armlength L f = 20, which illustrate the tank-treating motion. Fig. 21.
Scaled rotation frequencies as function of a rescaled Weissenberg numberfor various functionalities. Dashed and full lines correspond to ω/ ˙ γ = 1 / ω/ ˙ γ ∼ / Wi for large Wi, respectively. The inset shows the dependence ofthe rescaling factor ϕ on the functionality. From Ref. [142]. rates, ω/ ˙ γ decays as Wi − , which implies that the rotation frequency becomes independent of ˙ γ . A similar behavior has been observed for capsules at highshear rates [152].The presented results show that star polymers in shear flow show a veryrich structural and dynamical behavior. With increasing functionality, starsin flow change from linear-polymer-like to capsule-like behavior. These macro-molecules are therefore interesting candidates to tune the viscoelastic proper-ties of complex fluids. ulti-Particle Collision Dynamics The flow behavior of fluid droplets, capsules, vesicles and cells is of enormousimportance in science and technology. For example, the coalescence and break-up of fluid droplets is essential for emulsion formation and stability. Capsulesand vesicles are discussed and used as drug carriers. Red blood cells (RBC)flow in the blood stream, and may coagulate or be torn apart under unfa-vorable flow conditions. Red blood cells also have to squeeze through narrowcapillaries to deliver their oxygen cargo. White blood cells in capillary flowadhere to, roll along and detach again from the walls of blood vessels undernormal physiological conditions; in inflamed tissue, leukocyte rolling leads tofirm adhesion and induces an immunological response.These and many other applications have induced an intensive theoreticaland simulation activity to understand and predict the behavior of such soft,deformable objects in flow. In fact, there are some general, qualitative proper-ties in simple shear flow, which are shared by droplets, capsules, vesicles andcells. When the internal viscosity is low and when they are highly deformable,a tank-treading motion is observed (in the case of droplets for not too highshear rates), where the shape and orientation are stationary, but particles lo-calized at the interface or attached to the membrane orbit around the centerof mass with a rotation axis in the vorticity direction. On the other hand,for high internal viscosity or small deformability, the whole object performs atumbling motion, very much like a colloidal rod in shear flow. However, if wetake a more careful look, then the behavior of droplets, capsules, vesicles andcells is quite different. For example, droplets can break up easily at highershear rates, because their shape is determined by the interfacial tension; fluidvesicles can deform much more easily then capsules, since their membrane hasno shear elasticity; etc. We focus here on the behavior of fluid vesicles and redblood cells.
The modeling of lipid bilayer membranes depends very much on the lengthscale of interest. The structure of the bilayer itself or the embedding of mem-brane proteins in a bilayer are best studied with atomistic models of bothlipid and water molecules. Molecular dynamics simulations of such modelsare restricted to about 10 lipid molecules. For larger system sizes, coarse-grained models are required [153–155]. Here, the hydrocarbon chains of lipid molecules are described by short polymer chains of Lennard-Jones particles,which have a repulsive interaction with the lipid head groups as well as withthe water molecules, which are also modeled as single Lennard-Jones spheres.Very similar models, with Lennard-Jones interactions replaced by linear “soft”dissipative-particle dynamics (DPD) potentials, have also been employed in-tensively [156–160]. For the investigation of shapes and thermal fluctuationsof single- or multi-component membranes, the hydrodynamics of the solvent isirrelevant. In this case, it can be advantageous to use a solvent-free membranemodel , in which the hydrophobic effect of the water molecules is replaced by aneffective attraction among the hydrocarbon chains [161–164]. This approach isadvantageous in the case of membranes in dilute solution, because it reducesthe number of molecules—and thus the degrees of freedom to be simulated—by orders of magnitude. However, it should be noticed that the basic lengthscale of atomistic and coarse-grained or solvent-free models is still on the sameorder of magnitude.In order to simulate larger systems, such as giant unilamellar vesicles(GUV) or red blood cells (RBC), which have a radius on the order of sev-eral micrometers, a different approach is required. It has been shown that inthis limit the properties of lipid bilayer membranes are described very well bymodeling the membrane as a two-dimensional manifold embedded in three-dimensional space, with the shape and fluctuations controlled by the curvatureelasticity [165], H = (cid:90) dS κH , (98)where H = ( c + c ) / c and c , and the integral is over the whole membrane area. To makethe curvature elasticity amenable to computer simulations, it has to be dis-cretized. This can be done either by using triangulated surfaces [166, 167], orby employing particles with properly designed interactions which favor theformation of self-assembled, nearly planar sheets [168, 169]. In the latter case,both scalar particles with isotropic multi-particle interactions (and a curvatureenergy obtained from a moving least-squares method) [169] as well as particleswith an internal spin variable and anisotropic, multi-body forces [168] havebeen employed and investigated. In a dynamically-triangulated surface model [166,167,170–172] of vesicles andcells, the membrane is described by N mb vertices which are connected bytethers to form a triangular network of spherical topology, see Fig. 22. Thevertices have excluded volume and mass m mb . Two vertices connected bya bond have an attractive interaction, which keeps their distance below amaximum separation (cid:96) . A short-range repulsive interaction among all vertices ulti-Particle Collision Dynamics makes the network self-avoiding and prevents very short bond lengths. Thecurvature energy can be discretized in different ways [166, 173]. In particular,the discretization [173, 174] U cv = κ (cid:88) i σ i (cid:88) j ( i ) σ i,j r i,j r i,j (99)has been found to give reliable results in comparison with the continuumexpression (98). Here, the sum over j ( i ) is over the neighbors of a vertex i which are connected by tethers. The bond vector between the vertices i and j is r i,j , and r i,j = | r i,j | . The length of a bond in the dual lattice is σ i,j = r i,j [cot( θ ) + cot( θ )] /
2, where the angles θ and θ are opposite tobond ij in the two triangles sharing this bond. Finally, σ i = 0 . (cid:80) j ( i ) σ i,j r i,j is the area of the dual cell of vertex i . Fig. 22.
Triangulated-network model of a fluctuating membrane. All vertices haveshort-range repulsive interactions symbolized by hard spheres. Bonds represent at-tractive interactions with imply a maximum separation (cid:96) of connected vertices.From Ref. [175]. To model the fluidity of the membrane, tethers can be flipped betweenthe two possible diagonals of two adjacent triangles. A number ψN b of bond-flip attempts is performed with the Metropolis Monte Carlo method [173] attime intervals ∆t BF , where N b = 3( N mb −
2) is the number of bonds in thenetwork, and 0 < ψ <
Since the solubility of lipids in water is very low, the number of lipid moleculesin a membrane is essentially constant over typical experimental time scales.Also, the osmotic pressure generated by a small number of ions or macro-molecules in solution, which cannot penetrate the lipid bilayer, keeps the in-ternal volume essentially constant. The shape of fluid vesicles [176] is thereforedetermined by the competition of the curvature elasticity of the membrane,and the constraints of constant volume V and constant surface area S . In thesimplest case of vanishing spontaneous curvature, the curvature elasticity isgiven by Eq. (98). In this case, the vesicle shape in the absence of thermalfluctuations depends on a single dimensionless parameter, the reduced volume V ∗ = V /V , where V = (4 π/ R and R = ( S/ π ) / are the volume andradius of a sphere of the same surface area S , respectively. The calculated vesi-cle shapes are shown in Fig. 23. There are three phases. For reduced volumesnot too far from the sphere, elongated prolate shapes are stable. In a smallrange of reduced volumes of V ∗ ∈ [0 . , . Fig. 23.
Shapes of fluid vesicles as a function of the reduced volume V ∗ . D and D sto denote the discontinuous prolate-oblate and oblate-stomatocyte transitions,respectively. All shapes display rotational symmetry with respect to the verticalaxis. From Ref. [176]. Red blood cells have a biconcave disc shape, which can hardly be distinguishedfrom the discocyte shape of fluid vesicles with reduced volume V ∗ (cid:39) . ulti-Particle Collision Dynamics since a spectrin network is attached to the plasma membrane [181], whichhelps to retain the integrity of the cell in strong shear gradients or capillaryflow. Due to the spectrin network, the red blood cell membrane has a non-zeroshear modulus µ .The bending rigidity κ of RBCs has been measured by micropipette aspira-tion [182] and atomic force microscopy [183] to be approximately κ = 50 k B T .The shear modulus of the composite membrane, which is induced by thespectrin network, has been determined by several techniques; it is found to be µ = 2 × − N/m from optical tweezers manipulation [184], while the value µ = 6 × − N/m is obtained from micropipette aspiration [182]. Thus, thedimensionless ratio µR /κ (cid:39) / k el ( r i − r j ) . This tether network generates a shear modulus µ = √ k el . Solvent-free models, triangulated surfaces and other discretized curvaturemodels have the disadvantage that they do not contain a solvent, and thereforedo not describe the hydrodynamic behavior correctly. However, this apparentdisadvantage can be turned into an advantage by combining these models witha mesoscopic hydrodynamics technique. This approach has been employed fordynamically triangulated surfaces [37,180] and for meshless membrane modelsin combination with MPC [188], as well as for fixed membrane triangulationsin combination with both MPC [187] and the Lattice Boltzmann method(LBM) [189].
The solvent particles of the MPC fluid interact with the membrane in twoways to obtain an impermeable membrane with no-slip boundary conditions.First, the membrane vertices are included in the MPC collision procedure,as suggested for polymers in Ref. [68], compare Sec. 7.1. Second, the solventparticles are scattered with a bounce-back rule from the membrane surface.Here, solvent particles inside (1 ≤ i ≤ N in ) and outside ( N in < i ≤ N s )of the vesicle have to be distinguished. The membrane triangles are assumedto have a finite but very small thickness δ = 2 l bs . The scattering processis then performed at discrete time steps ∆t BS , so that scattering does notoccur exactly on the membrane surface, but the solvent particles can penetrateslightly into the membrane film [180]. A similar procedure has bee suggested inRef. [26] for spherical colloidal particles embedded in a MPC solvent. Particleswhich enter the membrane film, i.e. , which have a distance to the triangulatedsurface smaller than l bs , or interior particles which reach the exterior volumeand vice versa, are scattered at the membrane triangle with the closest centerof mass. Explicitly [180], v ( new ) s ( t ) = v s ( t ) − m mb m s + 3 m mb ( v s ( t ) − v tri ( t )) (100) v ( new ) tri ( t ) = v tri ( t ) + 2 m s m s + 3 m mb ( v s ( t ) − v tri ( t )) , (101)when ( v s ( t ) − v tri ( t )) · n tri <
0, where v s ( t ) and v tri ( t ) are the velocitiesof the solvent particle and of the center of mass of the membrane triangle,respectively, and n tri is the normal vector of the triangle, which is orientedtowards the outside (inside) for external (internal) particles.The bond flips provide a very convenient way to vary the membrane vis-cosity η mb , which increases with decreasing probability ψ for the selectionof a bond for a bond-flip attempt. The membrane viscosity has been deter-mined quantitatively from a simulation of a flat membrane in two-dimensionalPoiseuille flow. The triangulated membrane is put in a rectangular box of size L x × L y with periodic boundary conditions in the x -direction. The edge ver-tices at the lower and upper boundary ( y = ± L y /
2) are fixed at their positions.A gravitational force ( m mb g,
0) is applied to all membrane vertices to inducea flow. Rescaling of relative velocities is employed as a thermostat. Then,the membrane viscosity is calculated from η mb = ρ mb gL y / v max , where ρ mb is average mass density of the membrane particles, and v max the maximumvelocity of the parabolic flow profile. The membrane viscosity η mb , which isobtained in this way, is shown in Fig. 24. As ψ decreases, it takes longer andlonger for a membrane particle to escape from the cage of its neighbors, and η mb increases. This is very similar to the behavior of a hard-sphere fluid withincreasing density. Finally, for ψ = 0, the membrane becomes solid. ulti-Particle Collision Dynamics Fig. 24.
Dependence of the membrane viscosity η mb on the probability ψ for theselection of a bond for a bond-flip attempt, for a membrane with N mb = 1860vertices. From Ref. [180]. The dynamical behavior of fluid vesicles in simple shear flow has been stud-ied experimentally [190–193], theoretically [194–201], numerically with theboundary-integral technique [202, 203] or the phase-field method [203, 204],and with mesoscale solvents [37, 180, 205]. The vesicle shape is now deter-mined by the competition of the curvature elasticity of the membrane, theconstraints of constant volume V and constant surface area S , and the exter-nal hydrodynamic forces.Shear flow is characterized (in the absence of vesicles or cells) by the flowfield v = ˙ γy e x , where e x is a unit vector, compare Sec. 10.4. The controlparameter of shear flow is the shear rate ˙ γ , which has the dimension of aninverse time. Thus, a dimensionless, scaled shear rate ˙ γ ∗ = ˙ γτ can be defined,where τ is the longest relaxation time of a vesicle. It is given by τ = η R /κ ,where η is the solvent viscosity, R the average radius, and κ the bendingrigidity [206]. For ˙ γ ∗ <
1, the internal vesicle dynamics is fast compared tothe external perturbation, so that the vesicle shape is hardly affected by theflow field, whereas for ˙ γ ∗ >
1, the flow forces dominate and the vesicle is in anon-equilibrium steady state.One of the difficulties in theoretical studies of the hydrodynamic effectson vesicle dynamics is the no-slip boundary condition for the embedding fluidon the vesicle surface, which changes its shape dynamically under the effectof flow and curvature forces. In early studies, a fluid vesicle was thereforemodeled as an ellipsoid with fixed shape [194]. This simplified model is stillvery useful as a reference for the interpretation of simulation results.
The theory of Keller and Skalak [194] describes the hydrodynamic behavior ofvesicles of fixed ellipsoidal shape in shear flow, with the viscosities η in and η ofthe internal and external fluids, respectively. Despite of the approximationsneeded to derive the equation of motion for the inclination angle θ , whichmeasures the deviation of the symmetry axis of the ellipsoid with the flowdirection, this theory describes vesicles in flow surprisingly well. It has beengeneralized later [197] to describe the effects of a membrane viscosity η mb .The main result of the theory of Keller and Skalak is the equation ofmotion for the inclination angle [194], ddt θ = 12 ˙ γ {− B cos(2 θ ) } , (102)where B is a constant, which depends on the geometrical parameters of theellipsoid, on the viscosity contrast η ∗ in = η in /η , and the scaled membraneviscosity η ∗ mb = η mb / ( η R ) [180, 194, 197], B = f (cid:26) f + f − f ( η ∗ in −
1) + f f η ∗ mb (cid:27) , (103)where f , f , f , and f are geometry-dependent parameters. In the sphericallimit, B → ∞ . Eq. (102) implies the following behavior: • For
B >
1, there is a stationary solution, with cos(2 θ ) = 1 /B . This cor-responds to a tank-treading motion, in which the orientation of the vesicleaxis is time independent, but the membrane itself rotates around the vor-ticity axis. • For
B <
1, no stationary solution exists, and the vesicle shows a tumbling motion, very similar to a solid rod-like colloidal particle in shear flow. • The shear rate ˙ γ only determines the time scale, but does not affect thetank-treading or tumbling behavior. Therefore, a transition between thesetwo types of motion can only be induced by a variation of the vesicle shapeor the viscosities.However, the vesicle shape in shear flow is often not as constant as as-sumed by Keller and Skalak. In these situations, it is very helpful to comparesimulation results with a generalized Keller-Skalak theory, in which shape de-formation and thermal fluctuations are taken into account. Therefore, a phe-nomenological model has been suggested in Ref. [180], in which in addition tothe inclination angle θ a second parameter is introduced to characterize thevesicle shape and deformation, the asphericity [207] α = ( λ − λ ) + ( λ − λ ) + ( λ − λ ) R g , (104)where λ ≤ λ ≤ λ are the eigenvalues of the moment-of-inertia tensor andthe squared radius of gyration R g = λ + λ + λ . This implies α = 0 for ulti-Particle Collision Dynamics spheres (with λ = λ = λ ), α = 1 for long rods (with λ = λ (cid:28) λ ), and α = 1 / λ (cid:28) λ = λ ). The generalized Keller-Skalakmodel is then defined by the stochastic equations ζ α ddt α = − ∂F/∂α + A ˙ γ sin(2 θ ) + ζ α g α ( t ) (105) ddt θ = 12 ˙ γ {− B ( α ) cos(2 θ ) } + g θ ( t ) , (106)with Gaussian white noises g α and g θ , which are determined by (cid:104) g α ( t ) (cid:105) = (cid:104) g θ ( t ) (cid:105) = (cid:104) g α ( t ) g θ ( t (cid:48) ) (cid:105) = 0 , (cid:104) g α ( t ) g α ( t (cid:48) ) (cid:105) = 2 D α δ ( t − t (cid:48) ) , (107) (cid:104) g θ ( t ) g θ ( t (cid:48) ) (cid:105) = 2 D θ δ ( t − t (cid:48) ) , friction coefficients ζ α and ζ θ , and diffusion constants D α = k B T /ζ α and D θ = k B T /ζ θ . Note that ζ θ does not appear in Eq. (106); it drops out becausethe shear force is also caused by friction.The form of the stochastic equations (105) and (106) is motivated by thefollowing considerations. The first term in Eq. (105), ∂F/∂α , is the thermo-dynamic force due to bending energy and volume constraints; it is calculatedfrom the free energy F ( α ). The second term of Eq. (105) is the deformationforce due to the shear flow. Since the hydrodynamic forces elongate the vesi-cle for 0 < θ < π/ − π/ < θ < θ ) to leading order. The ampli-tude A is assumed to be independent of the asphericity α . ζ α and A can beestimated [205] from the results of a perturbation theory [199] in the quasi-spherical limit. Eq. (106) is adapted from Keller-Skalak theory. While B is aconstant in Keller-Skalak theory, it is now a function of the (time-dependent)asphericity α in Eq. (106). The theory of Keller and Skalak [194] predicts for fluid vesicles a transitionfrom tank-treading to tumbling with increasing viscosity contrast η in /η . Thishas been confirmed in recent simulations based on a phase-field model [203].The membrane viscosity η mb is also an important factor for the vesicle dy-namics in shear flow. For example, the membrane of red blood cells becomesmore viscous on aging [197, 208] or in diabetes mellitus [209]. Experimentsindicate that the energy dissipation in the membrane is larger than that in-side a red blood cell [196, 197]. Furthermore, it has been shown recently thatvesicles can not only be made from lipid bilayers, but also from bilayers ofblock copolymers [210]. The membrane viscosity of these “polymersomes” isseveral orders of magnitude larger than for liposomes, and can be changedover a wide range by varying the polymer chain length [211]. Fig. 25.
Snapshot of a discocyte vesicle in shear flow with reduced shear rate˙ γ ∗ = 0 .
92, reduced volume V ∗ = 0 .
59, membrane viscosity η ∗ mb = 0 and viscositycontrast η in /η = 1. The arrows represent the velocity field in the xz -plane. FromRef. [180]. A variation of the membrane viscosity can be implemented easily in dy-namically triangulated surface models of membranes, as explained in Sec. 11.2.2.An example of a discocyte in tank-treading motion, which is obtained by sucha membrane model [180], is shown in Fig. 25. Simulation results for the inclina-tion angle as a function of the reduced membrane viscosity η ∗ mb = η mb / ( η R )are shown in Fig. 26. This demonstrates the tank-treading to tumbling transi-tion of fluid vesicles with increasing membrane viscosity. The threshold shearrate decreases with decreasing reduced volume V ∗ , since with increasing de-viation from the spherical shape, the energy dissipation within the membraneincreases. Interestingly, the discocyte shape is less affected by the membraneviscosity than the prolate shape for V ∗ = 0 .
59, since it is more compact—incontrast to a vesicle with viscosity contrast η in /η >
1, where the prolateshape is affected less [180].Figure 26 also shows a comparison of the simulation data with results ofK-S theory for fixed shape, both without and with thermal fluctuations. Notethat there are no adjustable parameters. The agreement of the results of the-ory and simulations is excellent in the case of vanishing membrane viscosity, η mb = 0. For small reduced volumes, V ∗ (cid:39) .
6, the tank-treading to tum-bling transition is smeared out by thermal fluctuations, with an intermittenttumbling motion occurring in the crossover region. This behavior is capturedvery well by the generalized K-S model with thermal fluctuations. For largerreduced volumes and non-zero membrane viscosity, significant deviations oftheory and simulations become visible. The inclination angle θ is found to de-crease much more slowly with increasing membrane viscosity than expected ulti-Particle Collision Dynamics Fig. 26.
Average inclination angle (cid:104) θ (cid:105) as a function of reduced membrane viscosity η ∗ mb , for the shear rate ˙ γ ∗ = 0 .
92 and various reduced volumes V ∗ . Results arepresented for prolate (circles) and discocyte (squares) vesicles with V ∗ = 0 .
59, aswell as prolate vesicles with V ∗ = 0 .
66 (triangles), 0 .
78 (diamonds), 0 .
91 (crosses),and V ∗ = 0 .
96 (pluses). The solid and dashed lines are calculated by K-S theory,Eqs. (102) and (103), for prolate ( V ∗ = 0 .
59, 0 .
66, 0 .
78, 0 .
91, and 0 .
96) and oblateellipsoids ( V ∗ = 0 . V ∗ = 0 . V ∗ = 0 .
78, and the rotationalPeclet number ˙ γ/D θ = 600 (where D θ is the rotational diffusion constant). FromRef. [180]. theoretically. This is most likely due to thermal membrane undulations, whichare not taken into account in K-S theory. Fig. 27.
Temporal evolution of vesicle deformation α D and inclination angle θ , for V ∗ = 0 .
78 and η ∗ mb = 2 .
9. Here, α D = ( L − L ) / ( L + L ), where L and L are themaximum lengths in the direction of the eigenvectors of the gyration tensor in thevorticity plane. The solid (red) and dashed (blue) lines represent simulation datafor ˙ γ ∗ = 3 .
68 and 0 .
92 ( κ/k B T = 10 and 40, with ˙ γη R /k B T = 36 . γ ∗ = 1 .
8, 3 .
0, and 10 (from top to bottom). From Ref. [205].
Recently, of a new type of vesicle dynamics in shear flow has been observedexperimentally [192], which is characterized by oscillations of the inclinationangle θ with θ ( t ) ∈ [ − θ , θ ] and θ < π/
2. The vesicles were found totransit from tumbling to this oscillatory motion with increasing shear rate˙ γ . Simultaneously with the experiment, a “vacillating-breathing” mode forquasi-spherical fluid vesicles was predicted theoretically, based on a spherical-harmonics expansion of the equations of motion to leading order (withoutthermal fluctuations) [199]. This mode exhibits similar dynamical behavioras seen experimentally; however, it “coexists” with the tumbling mode, and ulti-Particle Collision Dynamics Fig. 28.
Dynamical phase diagram as a function of viscosity contrast η ∗ in = η in /η ,for η ∗ mb = 0 and various reduced volumes V ∗ , calculated from Eqs. (105), (106)without thermal noise. The tank-treading phase is located on the left-hand-sideof the dashed lines. The solid lines represent the tumbling-to-swinging transitions.From Ref. [205]. its orbit depends on the initial deformation, i.e., it is not a limit cycle. Fur-thermore, the shear rate appears only as the basic time scale, and thereforecannot induce any shape transitions. Hence it does not explain the tumbling-to-oscillatory transition seen in the experiments [192].Simulation data for the oscillatory mode—which has also been denoted“trembling” [192] or “swinging” [205] mode—are shown in Fig. 27. The sim-ulation results demonstrate that the transition can indeed be induced by in-creasing shear rate, and that it is robust to thermal fluctuations. Figure 27also shows that the simulation data are well captured by the generalized K-S model, Eqs. (105) and (106), which takes into account higher-order con-tributions in the curvature energy of a vesicle. The theoretical model cantherefore be used to predict the full dynamic phase diagram of prolate vesi-cles as a function of shear rate and membrane viscosity or viscosity contrast,compare Fig. 28. The swinging phase appears at the boundary between thetank-treading and the tumbling phase for sufficiently large shear rates. Thephase diagram explains under which conditions the swinging phase can bereached from the tumbling phase with increasing shear rate—as observed ex-perimentally [192].The generalized K-S model is designed to capture the vesicle flow behaviorfor non-spherical shapes sufficiently far from a sphere. For quasi-sphericalvesicles, a derivation of the equations of motion by a systematic expansionin the undulation amplitudes gives quantitatively more reliable results. Anexpansion to third order results in a phase diagram [200, 201], which agreesvery well with Fig. 28. Fig. 29.
Dynamical phase diagram of a vesicle in shear flow for reduced volume V ∗ = 0 .
59. Symbols show simulated parameter values, and indicate tank-treadingdiscocyte and tank-treading prolate ( ◦ ), tank-treading prolate and unstable disco-cyte ( (cid:52) ), tank-treading discocyte and tumbling (transient) prolate ( (cid:3) ), tumblingwith shape oscillation ( (cid:5) ), unstable stomatocyte (+), stable stomatocyte ( × ), andnear transition ( (cid:4) ). The dashed lines are guides to the eye. From Ref. [180]. Shear flow does not only induce different dynamical modes of prolate andoblate fluid vesicles, it can also induce phase transformations. The simplestcase is a oblate fluid vesicle with η mb = 0 and viscosity contrast η in /η = 1.When the reduced shear rate reaches ˙ γ ∗ (cid:39)
1, the discocyte vesicles arestretched by the flow forces into a prolate shape [37, 180, 202]. A similar tran-sition is found for stomatocyte vesicles, except that in this case a larger shearrate ˙ γ ∗ (cid:39) V ∗ (cid:39) .
6. Thus, a prolate vesicle in this regime starts tumbling;as the inclination angle becomes negative, shear forces push to shrink thelong axis of the vesicle; when this force is strong enough to overcome thefree-energy barrier between the prolate and the oblate phase, a shape trans- ulti-Particle Collision Dynamics formation can be induced, compare Fig. 30. The vesicle then remains in thestable tank-treading state.At higher shear rates, on intermittent behavior has been observed, in whichthe vesicle motion changes in irregular intervals between tumbling and tank-treading [180]. Fig. 30.
Time dependence of asphericity α and inclination angle θ , for ˙ γ ∗ = 1 . η ∗ mb = 1 .
62, and V ∗ = 0 .
59. The dashed lines are obtained from Eqs. (105) and(106), with ζ α = 100, A = 12, and B ( α ) = 1 . − . α . From Ref. [212]. At finite temperature, stochastic fluctuations of the membrane due to ther-mal motion affect the dynamics of vesicles. Since the calculation of thermalfluctuations under flow conditions requires long times and large membranesizes (in order to have a sufficient range of undulation wave vectors), simu-lations have been performed for a two-dimensional system in the stationarytank-treading state [213]. For comparison, in the limit of small deformationsfrom a circle, Langevin-type equations of motion have been derived, which arehighly nonlinear due to the constraint of constant perimeter length [213].The effect of the shear flow is to induce a tension in the membrane, whichreduces the amplitude of thermal membrane undulations. This tension can be extracted directly from simulation data for the undulation spectrum. Thereduction of the undulation amplitudes also implies that the fluctuations ofthe inclination angle θ get reduced with increasing shear rate. The theoryfor quasi-circular shapes predicts a universal behavior as a function of thescaled shear rate ˙ γ ∗ ∆ / κ/ ( R k B T ), where ˙ γ ∗ = ˙ γη R /κ is the reduced shearrate in two dimensions, and ∆ is the dimensionless excess membrane length.Theory and simulation results for the inclination angle as a function of thereduced shear rate are shown in Fig. 31. There are no adjustable parameters.The agreement is excellent as long as the deviations from the circular shapeare not too large [213]. Fig. 31.
Fluctuations of the inclination angle (cid:104) ∆θ (cid:105) / of a two-dimensional fluc-tuating vesicle in shear flow, as a function of scaled shear rate ˙ γ ∗ ∆ / κ/ ( R k B T ),where ∆ is the dimensionless excess membrane length. Symbols indicate simula-tion data for different internal vesicle areas A for fixed membrane length, with A ∗ ≡ A/πR = 0 .
95 (squares), A ∗ = 0 .
90 (triangles), A ∗ = 0 .
85 (stars), and A ∗ = 0 . ulti-Particle Collision Dynamics The deformation of single RBCs and single fluid vesicles in capillary flows werestudied theoretically by lubrication theories [214–216] and boundary-integralmethods [217–219]. In most of these studies, axisymmetric shapes which arecoaxial with the center of the capillary were assumed and cylindrical coordi-nates were employed. In order to investigate non-axisymmetric shapes as wellas flow-induced shape transformations, a fully three-dimensional simulationapproach is required.We focus here on the behavior of single red blood cells in capillary flow[187], as described by a triangulated surface model for the membrane (compareSec. 11.2.4) immersed in a MPC solvent (see Sec. 11.2.4). The radius of thecapillary, R cap , is taken to be slightly larger than the mean vesicle or RBCradius, R = (cid:112) S/ π , where S is the membrane area. Snapshots of vesicleand RBC shapes in flow are shown in Fig. 32 for a reduced volume of V ∗ =0 .
59, where the vesicle shape at rest is a discocyte. For sufficiently small flowvelocities, the discocyte shape is retained. However, the discocyte is found not in a coaxial orientation; instead the shortest eigenvalue of the gyration tensoris oriented perpendicular to the cylinder axis [187]. Since two opposite sides ofthe rim of the discocyte are closer to wall where the flow velocity is small, therotational symmetry is slightly disturbed and the top view looks somewhattriangular, see Fig. 32A. With increasing flow velocity, a shape transition to anaxisymmetric shape occurs. In the case of fluid vesicles this is a prolate shape,while in the case of RBCs a parachute shape is found. Such parachute shapesof red blood cells have previously been observed experimentally [209, 220].The fundamental difference between the flow behaviors of fluid vesiclesand red blood cells at high flow velocities is due to the shear elasticity of thespectrin network. Its main effect for µR /κ (cid:38) µR . In comparison, the shear stress in the parachute shape is muchsmaller.Some diseases, such as diabetes mellitus and sickle cell anemia, change themechanical properties of RBCs; a reduction of RBC deformability was foundto be responsible for an enhanced flow resistance of blood [221]. Therefore, itis very important to understand the relation of RBC elasticity and flow prop-erties in capillaries. The flow velocity at the discocyte-to-prolate transition offluid vesicles and at the discocyte-to-parachute transition is therefore shownin Fig. 33 as a function of the bending rigidity and the shear modulus, respec-tively. In both cases, an approximately linear dependence is obtained [187], v cm τR cap = 0 . µR k B T + 4 . κk B T . (108)
Fig. 32.
Snapshots of vesicles in capillary flow, with bending rigidity κ/k B T =20 and capillary radius R cap = 1 . R . (A) Fluid vesicle with discoidal shape atthe mean fluid velocity v m τ /R cap = 41, both in side and top views. (B) Elasticvesicle (RBC model) with parachute shape at v m τ /R cap = 218 (with shear modulus µR /k B T = 110). The blue arrows represent the velocity field of the solvent. (C)Elastic vesicle with slipper-like shape at v m τ /R cap = 80 (with µR /k B T = 110).The inside and outside of the membrane are depicted in red and green, respectively.The upper front quarter of the vesicle in (B) and the front half of the vesicle in (C)are removed to allow for a look into the interior; the black circles indicate the lineswhere the membrane has been cut in this procedure. Thick black lines indicate thewalls the cylindrical capillary. From Ref. [187]. This result suggests that parachute shapes of RBCs should appear for flowvelocities larger than v cm = 800 R cap /τ (cid:39) . µR /k B T . ulti-Particle Collision Dynamics Fig. 33.
Critical flow velocity v m of the discocyte-to-parachute transition of elasticvesicles and of the discocyte-to-prolate transition of fluid vesicles, as a function ofthe bending rigidity κ . From Ref. [187]. The flow of many red blood cells in wider capillaries has also been investi-gated by several simulation techniques. Discrete fluid-particle simulations—an extension of dissipative-particle dynamics (DPD)—in combination withbulk-elastic discocyte cells (in contrast to the membrane elasticity of real redblood cells) have been employed to investigate the dynamical clustering of redblood cells in capillary vessels [223,224]. An immersed finite-element model—acombination of the immersed boundary method for the solvent hydrodynam-ics [225] and a finite-element method to describe the membrane elasticity—hasbeen developed to study red blood cell aggregation [226]. Finally, it has beendemonstrated that the Lattice-Boltzmann method for the solvent in combina-tion with a triangulated mesh model with curvature and shear elasticity forthe membrane can be used efficiently to simulate RBC suspensions in widercapillaries [189].
One of the unique properties of soft matter is its viscoelastic behavior [13].Due to the long structural relaxation times, the internal degrees of freedomcannot relax sufficiently fast in an oscillatory shear flow already at moderatefrequencies, so that there is an elastic restoring force which pushes the systemback to its previous state. Well studied examples of viscoelastic fluids arepolymer solutions and polymer melts [13, 70].
The viscoelastic behavior of polymer solutions leads to many unusual flowphenomena, such as viscoelastic phase separation [227]. There is also a secondlevel of complexity in soft matter systems, in which a colloidal component isdispersed in a solvent, which is itself a complex fluid. Examples are spherical orrod-like colloids dispersed in polymer solutions. Shear flow can induce particleaggregation and alignment in these systems [228].It is therefore desirable to generalize the MPC technique to model vis-coelastic fluids, while retaining as much as possible of the computational sim-plicity of standard MPC for Newtonian fluids. This can be done by replacingthe point particles of standard MPC by harmonic dumbbells with spring con-stant K [229].As for point particles, the MPC algorithm consists of two steps, streamingand collisions. In the streaming step, within a time interval ∆t , the motion ofall dumbbells is governed by Newton’s equations of motion. The center-of-masscoordinate of each dumbbell follows a simple ballistic trajectory. The evolutionof the relative coordinates of dumbbell i , which consists of two monomers atpositions r i ( t ) and r i ( t ) with velocities v i ( t ) and v i ( t ), respectively, isdetermined by the harmonic interaction potential, so that r i ( t + ∆t ) − r i ( t + ∆t ) = A i ( t ) cos( ω ∆t ) + B i ( t ) sin( ω ∆t ) ; (109) v i ( t + ∆t ) − v i ( t + ∆t ) = − ω A i ( t ) sin( ω ∆t )+ ω B i ( t ) cos( ω ∆t ) , (110)with angular frequency ω = (cid:112) K/m . The amplitudes A i ( t ) and B i ( t ) aredetermined by the initial positions and velocities at time t . The collision stepis performed for the two point particles constituting a dumbbell in exactlythe same way as for MPC point-particle fluids. This implies, in particular,that the various collision rules of MPC, such as SRD, AT-a or AT+a, canall be employed also for simulations of viscoelastic solvents, depending on therequirements of the system under consideration. Since the streaming step isonly a little more time consuming and the collision step is identical, simula-tions of the viscoelastic MPC fluid can be performed with essentially the sameefficiency as for the standard point-particle fluid.The behavior of harmonic dumbbells in dilute solution has been studied indetail analytically [230]. These results can be used to predict the zero-shearviscosity η and the storage and loss moduli, G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) in oscillatoryshear with frequency ω , of the MPC dumbbell fluid. This requires the solventviscosity and diffusion constant of monomers in the solvent. Since the vis-coelastic MPC fluid consists of dumbbells only, the natural assumption is toemploy the viscosity η MP C and diffusion constant D of an MPC point-particlefluid of the same density. The zero-shear viscosity is then found to be [229] η = η MP C + ρ k B Tω H , (111)where ulti-Particle Collision Dynamics Fig. 34. (a) Storage G (cid:48) and (b) loss moduli G (cid:48)(cid:48) , as function of oscillation frequency ω on a double-logarithmic scale, for systems of dumbbells with various spring constantsranging from K = 0 . K = 1 .
0. Simulations are performed in two dimensions withthe SRD collision rule. The wall separation and the collision time are L y = 10 and ∆t = 0 .
02, respectively. From Ref. [229]. ω H = 4 Kζ = 4 DKk B T . (112)Similarly, the storage and loss modulus, and the average dumbbell extension,are predicted to be [229] G (cid:48) = ρk B T ω/ω H ) ω/ω H ) , (113) G (cid:48)(cid:48) = η MP C ω + ρk B T ω/ω H ω/ω H ) , (114)and (cid:104) r (cid:105)(cid:104) r (cid:105) eq = 1 + 23 ( ˙ γ/ω H ) . (115)Simulation data are shown in Fig. 34, together with the theoretical predica-tions (113) and (114). The comparison shows a very good agreement. Thisincludes not only the linear and quadratic frequency dependence of G (cid:48)(cid:48) and G (cid:48) for small ω , respectively, but also the leveling off when ω reaches ω H . Incase of G (cid:48)(cid:48) , there is quantitative agreement without any adjustable param-eters, whereas G (cid:48) is somewhat overestimated by Eq. (113) for small springconstants K . The good agreement of theory and simulations implies that thecharacteristic frequency decreases linearly with decreasing spring constant K and mean free path λ (since D ∼ λ ). A comparison of other quantities, suchas the zero-shear viscosity, shows a similar quantitative agreement [229]. In the short time since Malevanets and Kapral introduced multi-particle col-lision (MPC) dynamics as a particle-based mesoscale simulation techniquefor studying the hydrodynamics of complex fluids, there has been enormousprogress. It has been shown that kinetic theory can be generalized to cal-culate transport coefficients, several collision algorithms have been proposedand employed, and the method has been generalized to describe multi-phaseflows and viscoelastic fluids. The primary applications to date—which includestudies of the equilibrium dynamics and flow properties of colloids, polymers,and vesicles in solution—have dealt with mesoscopic particles embedded ina single-component Newtonian solvent. An important advantage of this algo-rithm is that it is very straightforward to model the dynamics for the embed-ded particles using a hybrid MPC-molecular dynamics simulations approach.The results of these studies are in excellent quantitative agreement with boththeoretical predictions and results obtained using other simulation techniques.How will the method develop in the future? This is of course difficult topredict. However, it seems clear that there will be two main directions, a fur-ther development of the method itself, and its application to new problems inSoft Matter hydrodynamics. On the methodological front, there are severalvery recent developments, like angular-momentum conservation, multi-phaseflows and viscoelastic fluids, which have to be explored in more detail. Itwill also be interesting to combine them to study, e.g. , multi-phase flows of ulti-Particle Collision Dynamics viscoelastic fluids. On the application side, the trend will undoubtedly be to-wards more complex systems, in which thermal fluctuations are important. Insuch systems, the method can play out its strengths, because the interactionsof colloids, polymers, and membranes with the mesoscale solvent can all betreated on the same basis. Acknowledgments
Financial support from the Donors of the American Chemical Society PetroleumResearch Fund, the National Science Foundation under grant No. DMR-0513393, the German Research Foundation (DFG) within the SFB TR6“Physics of Colloidal Dispersion in External Fields”, and the Priority Program“Nano- and Microfluidics” are gratefully acknowledged. We thank Elshad Al-lahyarov, Luigi Cannavacciuolo, Jens Elgeti, Reimar Finken, Ingo G¨otze, JensHarting, Martin Hecht, Hans Herrmann, Antonio Lamura, Kiaresch Mussaw-isade, Hiroshi Noguchi, Guoai Pan, Marisol Ripoll, Udo Seifert, Yu-Guo Tao,and Erkan T¨uzel for many stimulating discussions and enjoyable collabora-tions.
References
1. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. , 155 (1992).2. P. Espanol, Phys. Rev. E , 1734 (1995).3. P. Espanol and P. B. Warren, Europhys. Lett. , 191 (1995).4. G. R. McNamara and G. Zanetti, Phys. Rev. Lett. , 2332 (1988).5. X. Shan and H. Chen, Phys. Rev. E , 1815 (1993).6. X. He and L.-S. Luo, Phys. Rev. E , 6811 (1997).7. G. A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Oxford University Press, Oxford, 1994).8. F. J. Alexander and A. L. Garcia, Comp. in Phys. , 588 (1997).9. A. L. Garcia, Numerical Methods for Physics (Prentice Hall, 2000).10. U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett. , 1505 (1986).11. R. Adhikari, K. Stratford, M. E. Cates, and A. J. Wagner, Europhys. Lett. ,473 (2005).12. J. K. G. Dhont, An Introduction to Dynamics of Colloids (Elsevier, Amster-dam, 1996).13. R. G. Larson,
The Structure and Rheology of Complex Fluids (Oxford Univer-sity Press, Oxford, 1999).14. M. Ripoll, K. Mussawisade, R. G. Winkler, and G. Gompper, Europhys. Lett. , 106 (2004).15. M. Ripoll, K. Mussawisade, R. G. Winkler, and G. Gompper, Phys. Rev. E , 016701 (2005).16. M. Hecht, J. Harting, T. Ihle, and H. J. Herrmann, Phys. Rev. E , 011408(2005).17. J. T. Padding and A. A. Louis, Phys. Rev. E , 031402 (2006).18. A. Malevanets and R. Kapral, J. Chem. Phys. , 8605 (1999).6 G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler19. A. Malevanets and R. Kapral, J. Chem. Phys. , 7260 (2000).20. T. Ihle and D. M. Kroll, Phys. Rev. E , 066705 (2003).21. T. Ihle and D. M. Kroll, Phys. Rev. E , 020201(R) (2001).22. A. Mohan and P. S. Doyle, Macromolecules , 4301 (2007).23. J. M. Kim and P. S. Doyle, J. Chem. Phys. , 074906 (2006).24. N. Kikuchi, A. Gent, and J. M. Yeomans, Eur. Phys. J. E , 63 (2002).25. M. Ripoll, R. G. Winkler, and G. Gompper, Eur. Phys. J. E , 349 (2007).26. N. Kikuchi, C. M. Pooley, J. F. Ryder, and J. M. Yeomans, J. Chem. Phys. , 6388 (2003).27. T. Ihle and D. M. Kroll, Phys. Rev. E , 066706 (2003).28. T. Ihle, E. T¨uzel, and D. M. Kroll, Phys. Rev. E , 046707 (2005).29. J. A. Backer, C. P. Lowe, H. C. J. Hoefsloot, and P. D. Iedema, J. Chem. Phys. , 1 (2005).30. R. Kapral, Adv. Chem. Phys., to appear (2008).31. E. Allahyarov and G. Gompper, Phys. Rev. E , 036702 (2002).32. N. Noguchi, N. Kikuchi, and G. Gompper, Europhys. Lett. , 10005 (2007).33. T. Ihle, E. T¨uzel, and D. M. Kroll, Europhys. Lett. , 664 (2006).34. E. T¨uzel, T. Ihle, and D. M. Kroll, Math. Comput. Simulat. , 232 (2006).35. C. M. Pooley and J. M. Yeomans, J. Phys. Chem. B , 6505 (2005).36. A. Lamura, G. Gompper, T. Ihle, and D. M. Kroll, Europhys. Lett. , 319(2001).37. H. Noguchi and G. Gompper, Phys. Rev. Lett. , 258102 (2004).38. I. O. G¨otze, H. Noguchi, and G. Gompper, Phys. Rev. E , 046705 (2007).39. J. F. Ryder, Mesoscopic Simulations of Complex Fluids , Ph.D. thesis, Univer-sity of Oxford (2005).40. I. O. G¨otze, private communication (2007).41. J. Erpenbeck, Phys. Rev. Lett. , 1333 (1984).42. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (ClarendonPress, Oxford, 1987).43. J. H. Irving and J. G. Kirkwood, J. Chem. Phys. , 817 (1950).44. H. T. Davis, Statistical Mechanics of Phases, Interfaces, and Thin Films (VCHPublishers, Inc., 1996).45. E. T¨uzel, G. Pan, T. Ihle, and D. M. Kroll, Europhys. Lett. , 40010 (2007).46. F. Reif, Fundamentals of Statistical and Thermal Thysics (Mc-Graw Hill,1965).47. R. Zwanzig,
Lectures in Theoretical Physics , vol. 3 (Wiley, New York, 1961).48. H. Mori, Prog. Theor. Phys. , 423 (1965).49. H. Mori, Prog. Theor. Phys. , 399 (1965).50. J. W. Dufty and M. H. Ernst, J. Phys. Chem. , 7015 (1989).51. T. Ihle, E. T¨uzel, and D. M. Kroll, Phys. Rev. E , 035701(R) (2004).52. C. M. Pooley, Mesoscopic Modelling Techniques for Complex Fluids , Ph.D.thesis, University of Oxford (2003).53. E. T¨uzel, M. Strauss, T. Ihle, and D. M. Kroll, Phys. Rev. E , 036701 (2003).54. E. T¨uzel, Particle-bases mesoscale modeling of flow and transport in complexfluids , Ph.D. thesis, University of Minnesota (2006).55. T. Ihle and E. T¨uzel,
Static and dynamic properties of a particle-based algo-rithm for non-ideal fluids and binary mixtures (2006), cond-mat/0610350.56. B. J. Berne and R. Pecora,
Dynamic Light Scattering: With Applications toChemistry, Biology and Physics (Wiley, New York, 1976). ulti-Particle Collision Dynamics , 056702 (2006).58. H. Noguchi and G. Gompper, Europhys. Lett. , 36002 (2007).59. H. Noguchi and G. Gompper, Phys. Rev. E , 016706 (2008).60. T. Ihle, submitted to J. Phys.: Condens. Matter (2008).61. H. B. Callen, Thermodynamics (Wiley, New York, 1960).62. Y. Hashimoto, Y. Chen, and H. Ohashi, Comp. Phys. Commun. , 56 (2000).63. Y. Inoue, Y. Chen, and H. Ohashi, Comp. Phys. Commun. , 191 (2004).64. Y. Inoue, S. Takagi, and Y. Matsumoto, Comp. Fluids , 971 (2006).65. T. Sakai, Y. Chen, and H. Ohashi, Comp. Phys. Commun. , 75 (2000).66. T. Sakai, Y. Chen, and H. Ohashi, Phys. Rev. E , 031503 (2002).67. T. Sakai, Y. Chen, and H. Ohashi, Colloids Surf., A , 297 (2002).68. A. Malevanets and J. M. Yeomans, Europhys. Lett. , 231 (2000).69. J. Elgeti and G. Gompper, in NIC Symposium 2008 , edited by G. M¨unster,D. Wolf, and M. Kremer (Neumann Institute for Computing, J¨ulich, 2008),vol. 39 of
NIC series , pp. 53–61.70. M. Doi and S. F. Edwards,
The Theory of Polymer Dynamics (Oxford Univer-sity Press, Oxford, 1986).71. E. Falck, O. Punkkinen, I. Vattulainen, and T. Ala-Nissila, Phys. Rev. E ,050102(R) (2003).72. R. G. Winkler, K. Mussawisade, M. Ripoll, and G. Gompper, J. Phys.: Con-dens. Matter , S3941 (2004).73. K. Mussawisade, M. Ripoll, R. G. Winkler, and G. Gompper, J. Chem. Phys. , 144905 (2005).74. M. A. Webster and J. M. Yeomans, J. Chem. Phys. , 164903 (2005).75. E. Falck, J. M. Lahtinen, I. Vattulainen, and T. Ala-Nissila, Eur. Phys. J. E , 267 (2004).76. M. Hecht, J. Harting, and H. J. Herrmann, Phys. Rev. E , 051404 (2007).77. S. H. Lee and R. Kapral, J. Chem. Phys. , 11163 (2004).78. Y. Inoue, Y. Chen, and H. Ohashi, J. Stat. Phys. , 85 (2002).79. J. T. Padding, A. Wysocki, H. L¨owen, and A. A. Louis, J. Phys.: Condens.Matter , S3393 (2005).80. S. H. Lee and R. Kapral, Physica A , 56 (2001).81. A. Lamura and G. Gompper, Eur. Phys. J. E , 477 (2002).82. J. T. Padding, private communication (2007).83. C. Pierleoni and J.-P. Ryckaert, Phys. Rev. Lett. , 2992 (1991).84. B. D¨unweg and K. Kremer, Phys. Rev. Lett. , 2996 (1991).85. C. Pierleoni and J.-P. Ryckaert, J. Chem. Phys. , 8539 (1992).86. B. D¨unweg and K. Kremer, J. Chem. Phys. , 6983 (1993).87. C. Aust, M. Kr¨oger, and S. Hess, Macromolecules , 5660 (1999).88. E. S. Boek, P. V. Coveney, H. N. W. Lekkerkerker, and P. van der Schoot,Phys. Rev. E , 3124 (1997).89. P. Ahlrichs and B. D¨unweg, Int. J. Mod. Phys. C , 1429 (1998).90. P. Ahlrichs, R. Everaers, and B. D¨unweg, Phys. Rev. E , 040501 (2001).91. N. A. Spenley, Europhys. Lett. , 534 (2000).92. C. P. Lowe, A. F. Bakker, and M. W. Dreischor, Europhys. Lett. , 397 (2004).93. G. K. Batchelor, J. Fluid Mech. , 245 (1972).94. A. J. C. Ladd, Phys. Fluids , 481 (1997).95. K. H¨ofler and S. Schwarzer, Phys. Rev. E , 7146 (2000).96. J. T. Padding and A. A. Louis, Phys. Rev. Lett. , 220601 (2004).8 G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler97. H. Hayakawa and K. Ichiki, Phys. Rev. E , R3815 (1995).98. H. Hayakawa and K. Ichiki, Phys. Fluids , 481 (1997).99. M. H. Ernst, E. H. Hauge, and J. M. J. van Leeuwen, Phys. Rev. Lett. ,1254 (1970).100. M. Hecht, J. Harting, M. Bier, J. Reinshagen, and H. J. Herrmann, Phys. Rev.E , 021403 (2006).101. B. V. Derjaguin and D. P. Landau, Acta Phys. Chim. , 633 (1941).102. W. B. Russel, D. A. Saville, and W. Schowalter, Colloidal dispersions (Cam-bridge University Press, Cambridge, 1995).103. S. H. Lee and R. Kapral, J. Chem. Phys. , 214916 (2005).104. J. G. Kirkwood and J. Riseman, J. Chem. Phys. , 565 (1948).105. J. P. Erpenbeck and J. G. Kirkwood, J. Chem. Phys. , 909 (1958).106. E. P. Petrov, T. Ohrt, R. G. Winkler, and P. Schwille, Phys. Rev. Lett. ,258101 (2006).107. B. H. Zimm, J. Chem. Phys. , 269 (1956).108. S. H. Lee and R. Kapral, J. Chem. Phys , 214901 (2006).109. M. H. Ernst, E. H. Hauge, and J. M. J. van Leeuwen, Phys. Rev. A , 2055(1971).110. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press,London, 1986).111. B. Liu and B. D¨unweg, J. Chem. Phys. , 8061 (2003).112. P. Ahlrichs and B. D¨unweg, J. Chem. Phys. , 8225 (1999).113. M. Fixmann, J. Chem. Phys. , 1594 (1983).114. M. Schmidt and W. Burchard, Macromolecules , 210 (1981).115. W. H. Stockmayer and B. Hammouda, Pure & Appl. Chem. , 1373 (1984).116. P. E. Rouse, J. Chem. Phys. , 1272 (1953).117. P. G. de Gennes, Scaling Concepts in Polymer Physics (Cornell University,Ithaca, 1979).118. J. des Cloizeaux and G. Jannink,
Polymer Solutions: Their Modelling andStructure (Clarendon Press, Oxford, 1990).119. D. Ceperly, M. H. Kalos, and J. L. Lebowitz, Macromolecules , 1472 (1981).120. K. Kremer and K. Binder, Comput. Phys. Rep. , 261 (1988).121. R. G. Winkler, L. Harnau, and P. Reineker, Macromol. Theory Simul. , 1007(1997).122. R. Chang and A. Yethiraj, J. Chem. Phys. , 7688 (2001).123. N. Kikuchi, A. Gent, and J. M. Yeomans, Eur. Phys. J. E , 63 (2002).124. N. Kikuchi, J. F. Ryder, C. M. Pooley, and J. M. Yeomans, Phys. Rev. E ,061804 (2005).125. I. Ali and J. M. Yeomans, J. Chem. Phys. , 234903 (2005).126. I. Ali and J. M. Yeomans, J. Chem. Phys. , 8635 (2004).127. I. Ali, D. Marenduzzo, and J. F. D. Yeomans, Phys. Rev. Lett. , 208102(2006).128. N. Watari, M. Makino, N. Kikuchi, R. G. Larson, and M. Doi, J. Chem. Phys. , 094902 (2007).129. F. Brochard-Wyart, Europhys. Lett. , 105 (1993).130. F. Brochard-Wyart, H. Hervet, and P. Pincus, Europhys. Lett. , 511 (1994).131. F. Brochard-Wyart, Europhys. Lett. , 210 (1995).132. L. Cannavacciuolo, R. G. Winkler, and G. Gompper, EPL , 34007 (2008).133. U. S. Agarwal, A. Dutta, and R. A. Mashelkar, Chem. Eng. Sci. , 1693(1994). ulti-Particle Collision Dynamics , 2513 (2004).135. O. B. Usta, J. E. Butler, and A. J. C. Ladd, Phys. Fluids , 031703 (2006).136. R. Khare, M. D. Graham, and J. J. de Pablo, Phys. Rev. Lett. , 224505(2006).137. D. Stein, F. H. J. van der Heyden, W. J. A. Koopmans, and C. Dekker, Proc.Natl. Acad. Sci. USA , 15853 (2006).138. O. B. Usta, J. E. Butler, and A. J. C. Ladd, Phys. Rev. Lett. , 098301(2007).139. G. S. Grest, K. Kremer, and T. A. Witten, Macromolecules , 1376 (1987).140. C. N. Likos, Phys. Rep. , 267 (2001).141. D. Vlassopoulos, G. Fytas, T. Pakula, and J. Roovers, J. Phys.: Condens.Matter , R855 (2001).142. M. Ripoll, R. G. Winkler, and G. Gompper, Phys. Rev. Lett. , 188302 (2006).143. G. S. Grest, K. Kremer, S. T. Milner, and T. A. Witten, Macromolecules ,1904 (1989).144. D. R. Mikulencak and J. F. Morris, J. Fluid Mech. , 215 (2004).145. A. Link and J. Springer, Macromolecules , 464 (1993).146. R. E. Teixeira, H. P. Babcock, E. S. G. Shaqfeh, and S. Chu, Macromolecules , 581 (2005).147. R. G. Winkler, Phys. Rev. Lett. , 128301 (2006).148. P. G. de Gennes, J. Chem. Phys. , 5030 (1974).149. P. LeDuc, C. Haber, G. Bao, and D. Wirtz, Nature (London) , 564 (1999).150. D. E. Smith, H. P. Babcock, and S. Chu, Science , 1724 (1999).151. C. Aust, M. Kr¨oger, and S. Hess, Macromolecules , 8621 (2002).152. Y. Navot, Phys. Fluids , 1819 (1998).153. R. Goetz and R. Lipowsky, J. Chem. Phys. , 7397 (1998).154. R. Goetz, G. Gompper, and R. Lipowsky, Phys. Rev. Lett. , 221 (1999).155. W. K. den Otter and W. J. Briels, J. Chem. Phys. , 4712 (2003).156. J. C. Shillcock and R. Lipowsky, J. Chem. Phys. , 5048 (2002).157. L. Rekvig, B. Hafskjold, and B. Smit, Phys. Rev. Lett. , 116101 (2004).158. M. Laradji and P. B. S. Kumar, Phys. Rev. Lett. , 198105 (2004).159. V. Ortiz, S. O. Nielsen, D. E. Discher, M. L. Klein, R. Lipowsky, and J. Shill-cock, J. Phys. Chem. B , 17708 (2005).160. M. Venturoli, M. M. Sperotto, M. Kranenburg, and B. Smit, Phys. Rep. ,1 (2006).161. H. Noguchi and M. Takasu, J. Chem. Phys. , 9547 (2001).162. H. Noguchi, J. Chem. Phys. , 8130 (2002).163. O. Farago, J. Chem. Phys. , 596 (2003).164. I. R. Cooke, K. Kremer, and M. Deserno, Phys. Rev. E , 011506 (2005).165. W. Helfrich, Z. Naturforsch. , 693 (1973).166. G. Gompper and D. M. Kroll, J. Phys.: Condens. Matter , 8795 (1997).167. G. Gompper and D. M. Kroll, in Statistical Mechanics of Membranes andSurfaces , edited by D. R. Nelson, T. Piran, and S. Weinberg (World Scientific,Singapore, 2004), pp. 359–426, 2nd ed.168. J. M. Drouffe, A. C. Maggs, and S. Leibler, Science , 1353 (1991).169. H. Noguchi and G. Gompper, Phys. Rev. E , 021903 (2006).170. J.-S. Ho and A. Baumg¨artner, Europhys. Lett. , 295 (1990).171. D. M. Kroll and G. Gompper, Science , 968 (1992).0 G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler172. D. H. Boal and M. Rao, Phys. Rev. A , R6947 (1992).173. G. Gompper and D. M. Kroll, J. Phys. I France , 1305 (1996).174. C. Itzykson, in Proceedings of the GIFT Seminar, Jaca 85 , edited by J. Abad,M. Asorey, and A. Cruz (World Scientific, Singapore, 1986), pp. 130–188.175. G. Gompper, in
Soft Matter – Complex Materials on Mesoscopic Scales , editedby J. K. G. Dhont, G. Gompper, and D. Richter (Forschungszentrum J¨ulich,J¨ulich, 2002), vol. 10 of