A Numerical Model for Brownian Particles Fluctuating in Incompressible Fluids

We present a numerical method that consistently implements thermal fluctuations and hydrodynamic interactions to the motion of Brownian particles dispersed in incompressible host fluids. In this method, the thermal fluctuations are introduced as random forces acting on the Brownian particles. The hydrodynamic interactions are introduced by directly resolving the fluid motions with the particle motion as a boundary condition to be satisfied. The validity of the method has been examined carefully by comparing the present numerical results with the fluctuation-dissipation theorem whose analytical form is known for dispersions of a single spherical particle. Simulations are then performed for more complicated systems, such as a dispersion composed of many spherical particles and a single polymeric chain in a solvent.


Introduction
The dynamics of solid particles dispersed in host fluids is very complicated.Although computer simulations have been extensively used as a tool to investigate those systems, obtaining meaningful results is not yet straightforward even for the simplest case where the particles are monodisperse spheres and the host fluid is Newtonian.The main difficulty comes from the consistent treatment of the so-called hydrodynamic interaction (HI) between particles mediated by fluid motions, which are induced by the particle's motion.
The mathematical expression for HI is greatly simplified if the following assumptions are made: 1) the dispersed particles are all spherical; 2) HI acting among particles is pair-wise additive; and 3) the motions of the host fluid instantaneously follows the motions of the particles (the Stokes approximation).Then, the HI is expressed as a tensor that is a function of the particle's positions and velocities, without explicitly dealing with the fluid motions.The Oseen and the Rotne-Prager-Yamakawa (RPY) tensors are probably the most well known forms of such simplified HI functions; the former neglects the size of the particle, while the latter takes into account some corrections to the particle size.The Stokesian dynamics (SD) method 1) is a widely used numerical method along the lines of solving the tensor equations.It is based on the Langevin-type equations for particles implementing the RPY tensor and the lubrication correction.The latter is necessary when the distance between particles is small compared to the particle radius.
Completely different numerical approaches have been developed recently, [2][3][4][5][6] where the motions of host fluids are explicitly solved along with the motion of the particles so that the HI is directly computed without using the three assumptions described above. W refer those approaches as direct numerical simulation (DNS) approaches.An apparent benefit of using DNS approaches * E-mail: iwashita@cheme.kyoto-u.ac.jp for particle dispersions is that the long-time behavior of particle motions is reproduced accurately.For example, the velocity auto correlation function (VACF) of a fluctuating particle is expected to show a non-exponential power-law relaxation-known as the hydrodynamic longtime tail-if the thermal fluctuations are taken into account as well.
Recently, we also proposed an efficient DNS scheme for particle dispersions called the Smoothed Profile (SP) method, 7,8) where the discontinuous boundaries between particles and a fluid are smoothed out by using a continuous profile function, thereby achieving greater computational efficiency.The SP method has been successfully applied to some problems, such as the electrophoresis of charged colloidal particles; 9) however, the effects of thermal fluctuations (TF) are neglected.This is not a bad approximation for large/heavy particles, but it is insufficient for particles whose size is on the order of a micrometer or smaller; the coupling of TF and HI becomes crucial in these systems.Owing to some new experimental techniques that enable direct examination of the properties of Brownian particles fluctuating in a host fluid, several interesting phenomena have recently been reported, including the non-diffusive behavior of Brownian particles 10) and the rotational-translation coupling of a pair of spherical particles 11) where the coupling of TF and HI plays an essential role.
In the present study, our main goal is to implement TF and HI consistently within a DNS framework for particle dispersions.We aim to introduce TF to our original SP method to achieve this end.Since theoretical analyses have been well established for dilute dispersions, we first perform simulations for a dilute dispersion composed of a single spherical particle in a cubic box with the periodic boundary condition.The numerical results for the VACF are then compared with analytical solutions based on the fluctuation-dissipation theorem.After the validity of the method has been confirmed for this simple system, simulations are performed for dense dispersions composed of many spherical particles for which analytical solutions are unknown.We furthermore performed DNS simulations to examine the dynamics of a single polymeric chain fluctuating in a good solvent.The timecorrelation functions are calculated for each Rouse mode of the chain and compared with the predictions of the Rouse and Zimm models.

Simulation Method
We now present our implementation of TF based on the Langevin type approach.The equation governing the solvent with a density ρ f and a shear viscosity η is a modified Navier-Stokes equation with the incompressible condition ∇ • v = 0, where v(x, t), p(x, t) is the velocity and pressure field of the solvent.A smooth profile function 0 ≤ φ(x, t) ≤ 1 distinguishes between fluid and particle domains, namely φ = 1 in the particle domain and φ = 0 in the fluid domain.Those domains are separated by thin interfacial regions whose thickness is characterized by ξ.The body force φf p is introdued to ensures the rigidity of particles and the non-slip appropriate boundary condition at the fluid/particle interface.The mathematical expressions of φ and φf p are given in earlier papers 7,8) in detail.
We consider dispersions composed of N p spherical particles with a radius a.The motion of the ith particle is governed by Newton's equations of motion with stochastic forces: where R i , V i , and Ω i are the position, translational velocity, and rotational velocity of particles, respectively.M i and I i are the mass and the moment of inertia, and F H i and N H i are the hydrodynamic force and torque exerted by solvent on the particle. 7, 8)F C i is direct interparticle interaction such as Coulombic or the Lennard-Jones potential.In the present study, we used the truncated Lennard-Jones interaction: where The parameter ǫ characterizes the strength of interactions, and σ = 2a represents the diameter of particles.G V i and G Ω i are random forces and torques due to thermal fluctuations, which has the following properties where the square brackets denote taking an average over an equilibrium ensemble.α n represents the noise intensity for each degree of freedom of the translation (n = V ) and rotation (n = Ω) of the particles.Each noise intensity is controlled so that the variance of the translational and rotational velocity has a constant value; that is, where C 1 and C 2 are constant numbers.The time evolution of the noise intensity is described by α , where ∆t is the discrete time interval which plays a role of thermostat.
The temperature of the system is defined by the diffusive motion of the dispersed particles.The translational particle temperature k B T V is determined by the longtime diffusion coefficient D V of a spherical particle in the Stokes-Einstein relation k B T V = 6πηaD V where D V is obtained from computer simulation.Similarly, the rotational particle temperature k B T Ω can be determined by the rotational diffusion coefficient D Ω .
There are several advantages to the Langevin approach compared with the fluctuating hydrodynamics approach, for which TF is introduced in the Navier-Stokes equation as stochastic stresses to be defined on N d grid points of fluid simulations.One important advantage is the computational efficiency: while a N d spatial grid requires generating O(N d ) random numbers for the fluctuating hydrodynamics, our method requires O(N p ) random numbers for a dispersion composed of N p (≪ N d ) particles.Second, if we consider the solvent as a complex fluid-with arbitrary constitutive equation-, our method has another merit: in the fluctuating hydrodynamics approach, the friction tensor for complex fluids is required, but not here.Derivations of the friction tensor for complex fluids often have theoretical difficulties, and even if the friction tensor is obtained, there may be a large computational cost.

Test of the Simulation Method
In order to test our simulation method, we consider the translational motion of a particle dragged with a constant external force F 0 in a Newtonian fluid.With the force turned on, a dragged particle and the solvent have a steady state solution: the particle has a constant velocity along x-direction until t = 0.Then, at t > 0, the external force is removed and the particle and solvent relax toward a rest state due to dissipation.
Numerical simulation has been performed in three dimensions with periodic boundary conditions.The lattice spacing ∆ is taken to be the unit of length.Other units are defined so that we can set η = 1 and ρ f = 1 in eq.( 1), namely the units of time, mass, pressure are given by ρ f ∆ 2 /η, ρ f ∆ 3 , and η 2 /ρ f ∆ 2 , respectively.For a spherical particle, the radius a = 5, the thickness ξ = 2 and the particle density ρ p = 1 was used.In Fig. 1, we plotted − u(t)/F 0 for system size 32 3 , 64 3 and 128 3 where u(t) denotes the velocity in the x-direction of a dragged particle.The response function almost coincides with the analytical solution based on a time-dependent friction 12) when the system size is larger.The analytical solution is provided in detail in Appendix.For the long-time region, the regression of the velocity shows a power-law decay Bt −3/2 where B = 1/12ρ f (πν f ) 3/2 and the kinematic viscosity ν f = η/ρ f , which depends only on the hydrodynamic property of the solvent.
A basic relation between the relaxation response of a dragged particle and the velocity autocorrelation function of a Brownian particle is the fluctuation-dissipation theorem (FDT) where V i (t) • V i (0) /3 is the translational velocity autocorrelation function (VACF) for a Brownian particle and The relation holds that the response function of a dragged particle with external force F 0 is equal to the VACF of a particle in thermal equilibrium.
Under the same computational conditions as the relaxation experiment, the VACF in thermal equilibrium has been calculated at the volume fraction Φ = 0.002 and system size 64 3 .Figure 2 shows the VACF for a Brownian particle.We can calculate the long-time diffusion coefficient D V by the mean-square displacement of a Brownian particle; that is, lim When simulating a Brownian particle with the HI, the diffusion coefficient is affected by finite-size effects and is given by represents the effect of the periodic boundary condition. 13,14) aking the finite-size correction into account, the diffusion coefficient at infinite dilution is obtained as D V inf = D V K(Φ).The translational particle temperature is estimated to be k B T V = 6πηaD V inf ≃ 0.83.In Fig. 2, β V V i (t)•V i (0) /3 approaches to the analytical solution in the long-time region, shows the power-law decay Bt −3/2 , and gives the response function of a dragged particle.In the short-time region, a gap between simulation and the analytical solution is observed.The amplitude of a particle's velocity is related to the equipartition law of energy; that is, V 2 i = 3k B T V /M e where M e is the effective mass.The effective mass can be obtained via a hydrodynamic calculation as M e = M i + 0.5m 0 , where m 0 = 4πρ f a 3 /3 is the mass of fluid displaced by a spherical particle with radius a.In this simulation, the effective mass is estimated to be M e ≃ 4.2M i , which is notably greater than the hydrodynamic effective mass.The disagreement may be due to the influence of the artificial smoothed boundary used between a particle and fluid.Furthermore, we investigate Brownian particles in harmonic potentials.The harmonic potential with a spring constant k is introduced by adding where R eq i is its equilibrium position, to the equations of motion of particles.The mean-square displacement of a Brownian particle trapped in a harmonic potential is given by: where ∆r 2 (t) = |R i (t) − R i (0)| 2 /3.In the simulation for the particle number N p = 8, each particle is trapped at the grids of the fcc lattice by harmonic potentials and the direct interaction between particles is ignored.For k = 0.5, 1.0, 2.0 and 4.0, the lim t→∞ ∆r 2 (t) k/2 is plotted in Fig. 3.The result approaches k B T V for k → 0 and is consistent with the results obtained from the diffusion coefficient for a Brownian particle.
We investigate the response function of the rotational motion of a particle with a constant external torque N 0 in a similar fashion as the investigation of the par-ticle's translational motion.Figure 4 shows the relaxation response − ω/N 0 where ω(t) denotes the particle's rotational velocity.The response function almost coincides with the analytical solution, and a power-law decay Ct −5/2 is observed in the long-time region where C = π/32ρ f (πν f ) 5/2 .The FDT for the rotational velocity of a particle is also investigated.The rotational velocity autocorrelation function Ω i (t)•Ω i (0) /3 (RVACF) is given in Fig. 4, where Ω i (t) denotes the rotational velocity of a Brownian particle.In the Green-Kubo formula for the rotational velocity of the particle, the rotational particle temperature is estimated to be k B T Ω = 8πηa 3 D Ω ≃ 1.The RVACF also shows a power-law decay Ct −5/2 in the long-time region.The effective moment of inertia I e can be estimated by the value of the same time correlation as I e ≃ 2.6I i .

Applications
As a demonstration of a DNS incorporating thermal fluctuation of particles, our method is applied to a manyparticle system and a dilute polymeric chain.

Many particles system
In Fig. 5, the translational velocity autocorrelation function for each volume fraction Φ is presented.As the volume fraction is increased, it is found that the relaxation was more rapid than in low volume fractions, and the power-law long-time tail of the VACF gradually disappear in Fig. 5(a).Figure 5(b) shows that for Φ ≥ 0.4, the VACF has a negative overshoot, which represents an oscillative motion of a tagged particle due to a transient cage composed of surrounding particles.Compared with Brownian dynamics without HI, the decay of the correlation for a high volume fraction is much slower with HI, probably due to the lubrication force between particles.

A dilute polymeric chain in good solvent
The role of the HI is important in the dynamics of a dilute polymeric chain in a good solvent, and the simple theoretical model is known as the Zimm model. 15)This model treats the HI between beads as a hydrodynamic mobility matrix, such as the Oseen tensor.Some groups have studied a single polymeric chain with the HI using similar hybrid simulation methods, and their results are in agreement with Zimm theory. 3,16,17) Hre, we reexamine the validity of the Zimm model using our present DNS method since it is supposed to be more accurate than other methods used previously.
As a model of a polymeric chain, we study a beadspring model with a truncated Lennard-Jones potential and a finitely extensible non-linear elastic (FENE) potential 18) : where k c = 30ǫ/σ 2 , R 0 = 1.5σ and r is the distance between the neighboring beads.The position vector of a bead is described by R n (t) with n = 0, 1, ..., N ch − 1 where N ch denotes the total number of beads.
The static property of a polymeric chain is characterized by the static exponent ν, which is defined as R 2 G ∝ N 2ν ch b 2 for large N ch , where R G is the radius of gyration and b is the average bond length.The static exponent is related to the size of a polymeric chain, which is ν = 0.5 for a Gaussian chain and ν ≃ 0.6 for a selfavoided chain.The static exponent ν of a polymeric chain can be calculated via the static structure factor S(k).To analyze the relaxation dynamics of a polymeric chain, the real space motion is decomposed into a set of the Rouse modes (p = 0, 1, Within the approximations of the Zimm theory, the autocorrelation function of the Rouse mode decays exponentially as X p (t) • X p (0) / X 2 p = exp(−t/τ p ), where τ p denotes the relaxation time of the Rouse mode.The Zimm theory predicts the relation between the static exponent ν and the relaxation time τ p .The prediction is τ p ∼ p −3ν for the continuous model, and τ p ∼ p 2−3ν sin −2 (pπ/2N ch ) for the discrete model.
Figure 7 shows the normalized autocorrelation function of the Rouse mode for N ch = 10.The relaxation times are obtained by a fitting in the exponential shorttime regime t ∈ [50 : 1000].Figure 8 shows the mode (p-) dependence of the relaxation time.By fitting for p ≤ 5, the p-dependence of the relaxation time is estimated as τ p ∼ p −1.87 .In Fig. 8, the prediction of the discrete Zimm model for the p-dependence of τ p is also plotted using ν obtained from S(k).The numerical results show a good agreement with the prediction of the Zimm model.

Summary
We have developed a numerical method for consistently implementing thermal fluctuations and hydrodynamic interactions into models of the motions of Brownian particles dispersed in incompressible host fluids.We represented the thermal fluctuations by random forces acting on Brownian particles and the hydrodynamic interactions by directly resolving the fluid motions.The validity of the method has been examined carefully by com- paring the present numerical results with the fluctuationdissipation theorem for a dispersion of a single spherical particle.Simulations are then performed for dispersions of many spherical particles, and also for a polymeric chain in a fluid.In the former case, we found that the hydrodynamic long-time tail in the VACF-clearly observed for a single particle dispersion-becomes weak with increasing volume fraction of the particles.In the latter case, we found that our numerical results coincide quite well with the theoretical predictions of the Zimm model.

Appendix: A Dragged Particle with the Timedependent Friction
From the hydrodynamics, we can obtain the timedependent friction 12) ζ(t) = − 6πηau(t) + The first term is the standard Stokes resistance.The second term represents the additional mass, which is related to the acceleration of the particle.The third term represents the memory effect, which addresses the temporal decay of the fluid's momentum.
The equation of motion of a dragged particle with a

Fig. 1 .
Fig.1.The response functions, the translational velocity of a dragged particle at a constant external force F 0 , for system size 323 , 64 3 and 1283 .The response function is affected by the finitesize effect.Dotted line is the analytical solution of a drag problem with the time-dependent friction. 12)Dash-dotted line is the algebraic power Bt −3/2 where B = 1/12ρ f (πν f ) 3/2 .

Fig. 2 .
Fig. 2. Pluses(+): Translational velocity autocorrelation function for Brownian particles at β V = 1.2 and Φ = 0.002.The solid line is the analytical solution of a drag problem for the translational particle velocity.The dotted line is the algebraic power Bt −3/2 where B = 1/12ρ f (πν f ) 3/2 .The dashed line is the exponential decay obtained from Markovian Langevin equation.

Fig. 3 .
Fig. 3.The spring constant dependence of the limt→∞ ∆r 2 (t) k/2.k B T V is the translational particle temperature, which is obtained by the diffusion coefficient for a Brownian particle.

Fig. 4 .
Fig. 4. Crosses(×): Rotational velocity autocorrelation function for Brownian particles at β Ω = 1 and Φ = 0.002.Pluses(+): The relaxation response of rotational velocity of a particle with a constant external torque N 0 .The solid line is the analytical solution of a drag problem for the rotational particle velocity.The dashed line is the algebraic power Ct −5/2 where C = π/32ρ f (πν f ) 5/2 .The dashed-dotted line is the exponential decay obtained from Markovian Langevin equation.

Figure 6
displays the static structure factor for the bead numbers N ch = 10, 15.In the range of R −1 G ≪ k ≪ b −1 , S(k) obeys the scaling relation S(k) ∝ k −1/ν and can determine the static exponent ν ∼ 0.62 by fitting a power law to our data.

Fig. 7 .Fig. 8 .
Fig. 7. Normalized autocorrelation functions of the Rouse mode Xp for various mode p.The chain length is N ch = 10.