Journal of the Physical Society of Japan
Vol. 77, No. 7, July, 2008, 074007#2008 The Physical Society of Japan

A Numerical Model for Brownian Particles Fluctuatingin Incompressible Fluids
Takuya I WASHITA1 , Yasuya N AKAYAMA2 , and Ryoichi Y AMAMOTO1;3
1
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510
Department of Chemical Engineering, Kyushu University, Fukuoka 819-0395
3
CREST, Japan Science and Technology Agency, Kawaguchi, Saitama 332-0012
2(Received April 4, 2008; accepted May 19, 2008; published July 10, 2008)

We present a numerical method that consistently implements thermal ﬂuctuations and hydrodynamic
interactions to the motion of Brownian particles dispersed in incompressible host ﬂuids. In this method,
the thermal ﬂuctuations are introduced as random forces acting on the Brownian particles. The
hydrodynamic interactions are introduced by directly resolving the ﬂuid motions with the particle motion
as a boundary condition to be satisﬁed. The validity of the method has been examined carefully bycomparing the present numerical results with the ﬂuctuation–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 polymericchain in a solvent.
KEYWORDS: Brownian motion, hydrodynamic interaction, thermal fluctuation, particle dispersions, direct
numerical simulation, long-time tails, polymer dynamics
DOI: 10.1143/JPSJ.77.074007

1.

Introduction

The dynamics of solid particles dispersed in host ﬂuids
is very complicated. Although computer simulations havebeen 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 ﬂuid is Newtonian. The maindiﬃculty comes from the consistent treatment of the socalled hydrodynamic interaction (HI) between particles
mediated by ﬂuid motions, which are induced by theparticle’s motion.
The mathematical expression for HI is greatly simpliﬁed
if the following assumptions are made: 1) the dispersed
particles are all spherical; 2) HI acting among particles ispair-wise additive; and 3) the motions of the host ﬂuid
instantaneously follows the motions of the particles (theStokes 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 ﬂuid motions. The Oseen
and the Rotne–Prager–Yamakawa (RPY) tensors are probably the most well known forms of such simpliﬁed 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) method1) is a widely usednumerical method along the lines of solving the tensor
equations. It is based on the Langevin-type equations forparticles 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 diﬀerent numerical approaches have beendeveloped recently,2–6) where the motions of host ﬂuids are
explicitly solved along with the motion of the particles sothat the HI is directly computed without using the threeassumptions described above. We refer those approaches asE-mail: iwashita@cheme.kyoto-u.ac.jp

direct numerical simulation (DNS) approaches. An apparent
beneﬁt of using DNS approaches for particle dispersions isthat the long-time behavior of particle motions is reproduced
accurately. For example, the velocity autocorrelation function (VACF) of a ﬂuctuating particle is expected to show anon-exponential power-law relaxation — known as the hydrodynamic long-time tail — if the thermal ﬂuctuations (TF)
are taken into account as well.
Recently, we also proposed an eﬃcient DNS scheme
for particle dispersions called the Smoothed Proﬁle (SP)method,7,8) where the discontinuous boundaries between
particles and a ﬂuid are smoothed out by using a continuous
proﬁle function, thereby achieving greater computationaleﬃciency. The SP method has been successfully applied
to some problems, such as the electrophoresis of charged
colloidal particles;9) however, the eﬀects of TF are neglected. This is not a bad approximation for large/heavy
particles, but it is insuﬃcient 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 newexperimental techniques that enable direct examination of
the properties of Brownian particles ﬂuctuating in a hostﬂuid, several interesting phenomena have recently beenreported, including the non-diﬀusive behavior of Brownian
particles10) and the rotational-translation coupling of a pairof spherical particles11) 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 havebeen well established for dilute dispersions, we ﬁrst 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 thencompared with analytical solutions based on the ﬂuctuation–
dissipation theorem. After the validity of the method has

074007-1

J. Phys. Soc. Jpn., Vol. 77, No. 7

T. IWASHITA et al.

been conﬁrmed 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 thedynamics of a single polymeric chain ﬂuctuating in a good
solvent. The time-correlation functions are calculated foreach Rouse mode of the chain and compared with thepredictions of the Rouse and Zimm models.
2.

Simulation Method

We now present our implementation of TF based on theLangevin type approach. The equation governing the solventwith a density f and a shear viscosity  is a modiﬁed
Navier–Stokes equation
f ð@t v þ v  rvÞ ¼ r p þ r2 v þ f  f p

ð1Þ

with the incompressible condition r  v ¼ 0, wherevðx; tÞ; pðx; tÞ is the velocity and pressure ﬁeld of thesolvent. A smooth proﬁle function 0  ðx; tÞ  1 distinguishes between ﬂuid and particle domains, namely  ¼ 1 in
the particle domain and  ¼ 0 in the ﬂuid 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 ﬂuid/particle interface. The mathematical expressions of  and  f p are givenin earlier papers7,8) in detail.
We consider dispersions composed of Np sphericalparticles with a radius a. The motion of the ith particle isgoverned by Newton’s equations of motion with stochastic
forces:Mi V_ i ¼ FH þ FC þ GV ; R_ i ¼ V i ;
ð2Þ
i

i

i

I i  _ i ¼ N H þ G ;
i

i

ð3Þ

where Ri , V i , and i are the position, translational velocity,and rotational velocity of particles, respectively. Mi and I i
H
are the mass and the moment of inertia, and FH
i and N i are
the hydrodynamic force and torque exerted by solvent on theparticle.7,8) FCi is direct interparticle interaction such as
Coulombic or the Lennard-Jones potential. In the present
study, we used the truncated Lennard-Jones interaction:
8 " 12  6 #


>
< 4

þ  ðrij < 21=6 Þ,rijrij
UL J ðrij Þ ¼>:
0ðrij > 21=6 Þ,ð4Þwhere rij ¼ jRi  Rj j. The parameter  characterizes thestrength of interactions, and  ¼ 2a represents the diameterof particles. GVi and G
i are random forces and torques due
to TF, which has the following properties
hGni ðtÞi ¼ 0;hGni ðtÞGnj ð0Þi ¼ n I ðtÞ ij ;ð5Þwhere the angle brackets denote taking an average over an
equilibrium ensemble. n represents the noise intensity foreach degree of freedom of the translation (n ¼ V) androtation (n ¼ ) of the particles. Each noise intensity iscontrolled so that the variance of the translational androtational velocity has a constant value; that is, hV 2i i ¼ C1and h 2i i ¼ C2 , where C1 and C2 are constant numbers. The
time evolution of the noise intensity is described by V ðt þ22

tÞ ¼ V ðtÞe1V i ðtÞ=C1 and  ðt þ tÞ ¼  ðtÞe1 i ðtÞ=C2 ,where t is the discrete time interval which plays a roleof thermostat.
The temperature of the system is deﬁned by the diﬀusive
motion of the dispersed particles. The translational particle
temperature kB T V is determined by the long-time diﬀusion
coeﬃcient DV of a spherical particle in the Stokes–Einstein
relation kB T V ¼ 6 aDV where DV is obtained from computer simulation. Similarly, the rotational particle temperature kB T  can be determined by the rotational diﬀusion
coeﬃcient D .
There are several advantages to the Langevin approachcompared with the ﬂuctuating hydrodynamics approach, forwhich TF is introduced in the Navier–Stokes equation as
stochastic stresses to be deﬁned on N d grid points of ﬂuid
simulations. One important advantage is the computationaleﬃciency: while a N d spatial grid requires generating OðN d Þrandom numbers for the ﬂuctuating hydrodynamics, ourmethod requires OðNp Þ random numbers for a dispersion
composed of Np ( N d ) particles. Second, if we consider thesolvent as a complex ﬂuid — with arbitrary constitutive
equation —, our method has another merit: in the ﬂuctuatinghydrodynamics approach, the friction tensor for complex
ﬂuids is required, but not here. Derivations of the friction
tensor for complex ﬂuids often have theoretical diﬃculties,
and even if the friction tensor is obtained, there may be alarge computational cost.
3.

Test of the Simulation Method

In order to test our simulation method, we consider thetranslational motion of a particle dragged with a constantexternal force F0 in a Newtonian ﬂuid. With the force
turned on, a dragged particle and the solvent have a steadystate solution: the particle has a constant velocity along
x-direction until t ¼ 0. Then, at t > 0, the external force isremoved 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 latticespacing  is taken to be the unit of length. Other units
are deﬁned so that we can set  ¼ 1 and f ¼ 1 in eq. (1),namely the units of time, mass, pressure are given byf 2 =, f 3 , and 2 =f 2 , respectively. For a sphericalparticle, the radius a ¼ 5, the thickness  ¼ 2, and theparticle density p ¼ 1 was used. In Fig. 1, we plottedu_ðtÞ=F0 for system size 323 , 643 , and 1283 where uðtÞdenotes the velocity in the x-direction of a dragged particle.
The response function almost coincides with the analyticalsolution based on a time-dependent friction12) when thesystem size is larger. The analytical solution is provided
in detail in Appendix. For the long-time region, theregression of the velocity shows a power-law decay Bt3=2where B ¼ 1=12f ð f Þ3=2 and the kinematic viscosity f ¼
=f , which depends only on the hydrodynamic property ofthe solvent.
A basic relation between the relaxation response of a
dragged particle and the VACF of a Brownian particle is theﬂuctuation–dissipation theorem (FDT)

074007-2

1 duðtÞ
V¼
hV i ðtÞ  V i ð0Þi;
F0 dt
3

ð6Þ

J. Phys. Soc. Jpn., Vol. 77, No. 7

T. IWASHITA et al.

-3

-3

10

βV<V(t) V(0)>/3-1/F0 du(t)/dt

10

-4

10

10-5
L =32
L =64
L =128analytical solution
Bt -1.5

10-6
-7

10

10-1

100

101

10-4

10-5
analytical solution
Markovian VACF
Bt -1.5
βV <V(t) V(0)>/3

-6

10

10-1

102

10 0

101

10 2

t

t

where hV i ðtÞ  V i ð0Þi=3 is the translational VACF for a
Brownian particle and V ¼ 1=kB T V . The relation holds thatthe response function of a dragged particle with external
force F0 is equal to the VACF of a particle in thermal
equilibrium. Under the same computational conditions as therelaxation experiment, the VACF in thermal equilibrium has
been calculated at the volume fraction  ¼ 0:002 andsystem size 643 . Figure 2 shows the VACF for a Brownian
particle. We can calculate the long-time diﬀusion coeﬃcient
DV by the mean-square displacement of a Brownian particle;
that is, limt!1 hjRi ðtÞ  Ri ð0Þj2 i ¼ 6DV t.
When simulating a Brownian particle with the HI, thediﬀusion coeﬃcient is aﬀected by ﬁnite-size eﬀects andis given by DV ¼ kB T V =6 aKðÞ where KðÞ represents
the eﬀect of the periodic boundary condition.13,14) Taking
the ﬁnite-size correction into account, the diﬀusion coefﬁcient at inﬁnite dilution is obtained as DVinf ¼ DV KðÞ.
The translational particle temperature is estimated to be
kB T V ¼ 6 aDVinf ’ 0:83. In Fig. 2, V hV i ðtÞ  V i ð0Þi=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 isrelated to the equipartition law of energy; that is, hV 2i i ¼
3kB T V =Me where Me is the eﬀective mass. The eﬀective
mass can be obtained via a hydrodynamic calculation as
Me ¼ Mi þ 0:5m0 , where m0 ¼ 4 f a3 =3 is the mass of
ﬂuid displaced by a spherical particle with radius a. In thissimulation, the eﬀective mass is estimated to be Me ’
4:2Mi , which is notably greater than the hydrodynamic
eﬀective mass. The disagreement may be due to theinﬂuence of the artiﬁcial smoothed boundary used between
a particle and ﬂuid.
Furthermore, we investigate Brownian particles in harmonic potentials. The harmonic potential with a springeq
constant k is introduced by adding Fex
i ¼ kðRi  Ri Þ,eqwhere Ri is its equilibrium position, to the equations of
motion of particles. The mean-square displacement of a

Fig. 2. (Color online) Pluses (þ): Translational velocity autocorrelation
function for Brownian particles at V ¼ 1:2 and  ¼ 0:002. The solidline 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=12f ð f Þ3=2 . The dashed line is the exponential decay obtained
from Markovian Langevin equation.

1
0.8<Δ r 2 > k / 2

Fig. 1. (Color online) The response functions, the translational velocity
of a dragged particle at a constant external force F0 , for system size 323 ,
643 , and 1283 . The response function is aﬀected by the ﬁnite-size eﬀect.
Dotted line is the analytical solution of a drag problem with the timedependent friction.12) Dash-dotted line is the algebraic power Bt3=2where B ¼ 1=12f ð f Þ3=2 .

0.6
0.4
0.2
0
0

simulation
V
k BT
0.5

1

1.5

2

2.5

3

3.5

4

4.5

k
Fig. 3. (Color online) The spring constant dependence of thelimt!1 hr 2 ðtÞik=2. kB T V is the translational particle temperature, which
is obtained by the diﬀusion coeﬃcient for a Brownian particle.

Brownian particle trapped in a harmonic potential is givenby:lim hr 2 ðtÞi ¼

t!1

2kB T
;
k

ð7Þ

where hr 2 ðtÞi ¼ hjRi ðtÞ  Ri ð0Þj2 i=3. In the simulation forthe particle number Np ¼ 8, each particle is trapped at thegrids of the fcc lattice by harmonic potentials and the directinteraction between particles is ignored. For k ¼ 0:5, 1.0,
2.0, and 4.0, the limt!1 hr 2 ðtÞik=2 is plotted in Fig. 3.
The result approaches kB T V for k ! 0 and is consistent withthe results obtained from the diﬀusion coeﬃcient for a
Brownian particle.
We investigate the response function of the rotationalmotion of a particle with a constant external torque N0 in asimilar fashion as the investigation of the particle’s translational motion. Figure 4 shows the relaxation response!_ =N0 where !ðtÞ denotes the particle’s rotational velocity.
The response function almost coincides with the analyticalsolution, and a power-law decay Ct5=2 is observed in thelong-time region where C ¼ =32f ð f Þ5=2 . The FDT for

074007-3

10

-3

10

-4

10

-5

10

10

T. IWASHITA et al.

10-6βV <V(t) V(0)>/3

βΩ <Ω(t) · Ω(0)>/3, -1/N0 dω /dt

J. Phys. Soc. Jpn., Vol. 77, No. 7

analytical solution
Ct -2.5-1/N 0 dω(t)/dt
βΩ <Ω(t) · Ω(0)>/3-7

Markovian RVACF

10-8
10

-1

100

101

-2

(a)
10-3

10-4

10

analytical solution
Φ=0.02
Φ=0.1
Φ=0.3
Φ=0.4
Φ=0.5
Bt -1.5

-5

10-6

102

10-1

100

102

t

t

the rotational velocity of a particle is also investigated. The
rotational velocity autocorrelation function h i ðtÞ  i ð0Þi=3(RVACF) is given in Fig. 4, where i ðtÞ denotes therotational 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
kB T  ¼ 8 a3 D ’ 1. The RVACF also shows a powerlaw decay Ct5=2 in the long-time region. The eﬀective
moment of inertia I e can be estimated by the value of thesame time correlation as I e ’ 2:6I i .
Applications

As a demonstration of a DNS incorporating TF ofparticles, our method is applied to a many-particle system
and a dilute polymeric chain.
4.1 Many particles system
In Fig. 5, the translational VACF for each volume fraction
 is presented. As the volume fraction is increased, it isfound 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 thatfor   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 ofthe correlation for a high volume fraction is much slower
with HI, probably due to the lubrication force between
particles.
4.2 A dilute polymeric chain in good solvent
The role of the HI is important in the dynamics of a dilutepolymeric chain in a good solvent, and the simple theoreticalmodel is known as the Zimm model.15) This model treats the
HI between beads as a hydrodynamic mobility matrix, suchas the Oseen tensor. Some groups have studied a single
polymeric chain with the HI using similar hybrid simulation

0.0001

βV <V(t) V(0)>/3

Fig. 4. (Color online) 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 constantexternal torque N0 . The solid line is the analytical solution of a dragproblem for the rotational particle velocity. The dashed line is thealgebraic power Ct5=2 where C ¼ =32f ð f Þ5=2 . The dashed-dottedline is the exponential decay obtained from Markovian Langevinequation.

4.

101

Φ=0.02
Φ=0.1
Φ=0.3
Φ=0.4
Φ=0.5

(b)
8e-05
6e-05
4e-05
2e-05
0
-2e-05

0

50

100

150

200

250

t
Fig. 5. (Color online) The translational velocity autocorrelation function
for each volume fraction. The parameters of the simulation are f ¼ 1,
 ¼ 1, p ¼ 1, a ¼ 5,  ¼ 2, V ¼ 1, and  ¼ 1. (a) Log–log plot. (b)Normal plot.

methods, and their results are in agreement with Zimm
theory.3,16,17) Here, we reexamine the validity of the Zimm
model using our present DNS method since it is supposed tobe more accurate than other methods used previously.
As a model of a polymeric chain, we study a bead-springmodel with a truncated Lennard-Jones potential and aﬁnitely extensible non-linear elastic (FENE) potential:18)
1
UF ðrÞ ¼  kc R20 ln½1  ðr=R0 Þ2 ;
2ð8Þ

where kc ¼ 30= 2 , R0 ¼ 1:5, and r is the distance
between the neighboring beads. The position vector of abead is described by Rn ðtÞ with n ¼ 0; 1; . . . ; Nch  1 where
Nch denotes the total number of beads.
The static property of a polymeric chain is characterized
2 2
by the static exponent , which is deﬁned as hR2G i / Nchbfor large Nch , where RG is the radius of gyration and b is theaverage bond length. The static exponent is related to thesize of a polymeric chain, which is ¼ 0:5 for a Gaussian
chain and ’ 0:6 for a self-avoided chain. The static
exponent of a polymeric chain can be calculated via thestatic structure factor SðkÞ. Figure 6 displays the static
structure factor for the bead numbers Nch ¼ 10, 15. In the
1range of R1
G  k  b , SðkÞ obeys the scaling relation
1=SðkÞ / kand can determine the static exponent
0:62
by ﬁtting a power law to our data.

074007-4

J. Phys. Soc. Jpn., Vol. 77, No. 7

T. IWASHITA et al.
1

S(k)

2

< Xp (t) · Xp (0)>/< X p (0)>

101

10

0

Nch =10
Nch =15
10-2

p=1
p=2
p=3
p=4
p=5
p=6
p=7
p=8
p=9
0

10-1

100

200

300

400

100

500

600

700

800

900

t

k
Fig. 6. (Color online) The static structure factor of a chain, for Nch ¼ 10(pluses), 15 (crosses). The parameters of the simulation are f ¼ 1,  ¼ 1,
p ¼ 1, a ¼ 5,  ¼ 2, V ¼ 1, and  ¼ 1. System size is 643 grids.

Fig. 7. (Color online) Normalized autocorrelation functions of the Rousemode Xp for various mode p. The chain length is Nch ¼ 10.

Within the approximations of the Zimm theory, theautocorrelation function of the Rouse mode decays exponentially as hX p ðtÞ  X p ð0Þi=hX2p i ¼ 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 forthe continuous model, and p p23 sin2 ðp =2Nch Þ forthe discrete model.
Figure 7 shows the normalized autocorrelation function of
the Rouse mode for Nch ¼ 10. The relaxation times are
obtained by a ﬁtting in the exponential short-time regimet 2 ½50 : 1000. Figure 8 shows the mode (p-) dependence ofthe relaxation time. By ﬁtting for p  5, the p-dependenceof the relaxation time is estimated as p p1:87 . In
Fig. 8, the prediction of the discrete Zimm model for thep-dependence of p is also plotted using obtained fromSðkÞ. The numerical results show a good agreement with theprediction of the Zimm model.
5.

Summary

We have developed a numerical method for consistently
implementing thermal ﬂuctuations and hydrodynamic interactions into models of the motions of Brownian particles
dispersed in incompressible host ﬂuids. We represented thethermal ﬂuctuations by random forces acting on Brownian

5

10

4

10

3

τp

To analyze the relaxation dynamics of a polymeric chain,
the real space motion is decomposed into a set of the Rousemodes (p ¼ 0; 1; . . . ; Nch  1),Nch
1 X
Rn ðtÞ cos½ðn  1=2Þp =Nch :ð9Þ
X p ðtÞ ¼
Nch n¼1

10102

present simulation
p -1.87
discrete Zimm model
100

101
p

Fig. 8. (Color online) Rouse mode p-dependence of the relaxation time
for chain length Nch ¼ 10. The dashed line is a ﬁtting power law tosimulation data for p  5. The theoretical prediction of the discrete Zimm
model is indicated as a dotted line, p p23 sin2 ðp =2Nch Þ.

particles and the hydrodynamic interactions by directly
resolving the ﬂuid motions. The validity of the method has
been examined carefully by comparing the present numerical results with the ﬂuctuation–dissipation theorem for adispersion of a single spherical particle. Simulations arethen performed for dispersions of many spherical particles,
and also for a polymeric chain in a ﬂuid. 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 theparticles. In the latter case, we found that our numericalresults coincide quite well with the theoretical predictions of
the Zimm model.

Appendix: A Dragged Particle with the Time-Dependent Friction
From the hydrodynamics, we can obtain the time-dependent friction12)


Z
2
u_ðsÞ
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ t
f a3 u_ðtÞ þ 6a2
ðtÞ ¼  6 auðtÞ þ
ds pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ :
f
3
ts
1

ðA:1Þ

The ﬁrst term is the standard Stokes resistance. The second term represents the additional mass, which is related to theacceleration of the particle. The third term represents the memory eﬀect, which addresses the temporal decay of the ﬂuid’s
momentum.
The equation of motion of a dragged particle with a time-dependent frictional force  and a constant external force F0 is
074007-5

J. Phys. Soc. Jpn., Vol. 77, No. 7

T. IWASHITA et al.

described by M u_ðtÞ ¼ ðtÞ þ F0 , where M is the mass of a particle. After the particle has reached a steady state (a constantvelocity), the external force is removed. Then, the regression of the velocity can be written as
Z
pﬃﬃﬃ
duðtÞ
F0 1 dy 0 y exp½yðjtj= B Þ
;
ðA:2Þ
¼
dt
Meff 0j1  yj2 þ 02 y
where 0 ¼ ð9f =ð2p þ f ÞÞ1=2 , Meff ¼ M þ 0:5m0 , p the particle density, and B ¼ M=6 a. From the ﬂuctuation–
dissipation relation ½1=ð3kB T V ÞhV i ðtÞ  V i ð0Þi ¼ ð1=F0 ÞduðtÞ=dt, a hydrodynamic velocity autocorrelation function of a
Brownian particle can be calculated exactly.
Similarly, the regression of the rotational velocity ! of a spherical particle with a moment of inertia I can be written as
3

2
Z
d!ðtÞ
N0 1 1 dy
yt
y3=2
¼
exp

ðA:3Þ

2 7;

 2
6
dt
8 a3 f 0 3
1
f 4
5
r
r
þ
1
y þy 1 y
3ff
where N0 is a constant external torque,

1)
2)
3)
4)
5)
6)
7)
8)
9)
10)

2
f ¼ a = f and

r ¼ I=8J. F. Brady and G. Bossis: Annu. Rev. Fluid Mech. 20 (1988) 111.
A. Malevanets and R. Kapral: J. Chem. Phys. 112 (2000) 7260.
P. Ahlrichs and B. Dünweg: J. Chem. Phys. 111 (1999) 8225.
J. T. Padding and A. A. Louis: Phys. Rev. E 74 (2006) 031402.
P. J. Atzberger, P. R. Kramer, and C. S. Peskin: J. Comput. Phys. 224(2007) 1255.
N. Sharma and N. A. Patankar: J. Comput. Phys. 201 (2004) 466.
Y. Nakayama and R. Yamamoto: Phys. Rev. E 71 (2005) 036707.
Y. Nakayama, K. Kim, and R. Yamamoto: to be published in Eur.
Phys. J. E (2008) (cond-mat/0601322).
K. Kim, Y. Nakayama, and R. Yamamoto: Phys. Rev. Lett. 96 (2006)
208302.
B. Lukić, S. Jeney, C. Tischer, A. J. Kulik, L. Forró, and E.-L. Florin:
Phys. Rev. Lett. 95 (2005) 160601.

a3 .

11) S. Martin, M. Reichert, H. Stark, and T. Gisler: Phys. Rev. Lett. 97(2006) 248301.
12) E. M. Lifshitz and L. D. Landau: Fluid Mechanics (Addison-Wesley,
Reading, MA, 1959).
13) A. A. Zick and G. M. Homsy: J. Fluid Mech. 115 (1982) 13.
14) H. Hashimoto: J. Fluid Mech. 5 (1959) 317.
15) M. Doi and S. F. Edwards: The Theory of Polymer Dynamics (Oxford
University Press, Oxford, U.K., 1986).
16) K. Mussawisade, M. Ripoll, R. G. Winkler, and G. Gompper: J. Chem.
Phys. 123 (2005) 144905.
17) G. Giupponi, G. De Fabritiis, and P. V. Coveney: J. Chem. Phys. 126(2007) 154903.
18) K. Kremer and G. S. Grest: J. Chem. Phys. 92 (1990) 5057.

074007-6