FULL P APERS

Journal of the Physical Society of Japan 82 (2013) 084003
http://dx.doi.org/10.7566/JPSJ.82.084003

Drag Coeﬃcient of a Rigid Spherical Particle
in a Near-Critical Binary Fluid Mixture
Ryuichi OKAMOTO1 , Youhei FUJITANI2 y, and Shigeyuki KOMURA3 z
2

1
Fukui Institute for Fundamental Chemistry, Kyoto University, Kyoto 606-8103, Japan
School of Fundamental Science and Technology, Keio University, Yokohama 223-8522, Japan3
Department of Chemistry, Tokyo Metropolitan University, Hachioji, Tokyo 192-0397, Japan(Received February 1, 2013; accepted June 3, 2013; published online July 22, 2013)

We calculate the drag coeﬃcient of a rigid spherical particle in an incompressible binary ﬂuid mixture. A weak
preferential attraction is assumed between the particle surface and one of the ﬂuid components, and the diﬀerence inthe viscosity between the two components is neglected. Using the Gaussian free-energy functional and solving thehydrodynamic equation explicitly, we can show that the preferential attraction makes the drag coeﬃcient larger as thebulk correlation length becomes longer. The dependence of the deviation from the Stokes law on the correlation length,
when it is short, turns out to be much steeper than the previous estimates.
KEYWORDS: colloid dynamics, wetting, two-component fluid, time-dependent Ginzburg–Landau equation, model H

1.

Introduction

Studies on colloid particles in a ﬂuid have contributed toprogress in fundamental sciences in not only the statics but
also the dynamics.1–3) A system of particles exhibits various
thermodynamic phases, while a particle exhibits Brownian
motion. If it moves translationally with a suﬃciently low
speed in a quiescent ﬂuid, a particle suﬀers a drag force, themagnitude of which is proportional to the particle speed. The
constant of the proportion  is called the drag coeﬃcient.
For a rigid spherical particle with the radius r0 in a threedimensional incompressible ﬂuid with the viscosity , thedrag coeﬃcient is given by  ¼ 6r0 , which is called the
Stokes law.4) In the theory of Brownian motion, the diﬀusion
constant of the particle is found to be kB T =ð6r0 Þ, where kB
and T are the Boltzmann constant and the temperature of theﬂuid, respectively.5,6)
It is probable, in a binary ﬂuid mixture, that the surface
of a colloid particle preferentially attracts one of the ﬂuid
components. This results in a concentration gradient (oradsorption layer) near the colloid particle. Several experimental groups have observed adsorption-induced an aggregation of colloid particles slightly outside the coexistence
curve of water-organic solvent mixtures.7–12) The magnitude
and range of the heterogeneity in the concentration tend to
become larger when the critical demixing point is approached. Furthermore, near the critical point, the concentration proﬁle and inter-particle interaction energy exhibit
universal behavior, i.e., it does not depend on the
microscopic details of the surface and solvent molecules if
the physical quantities are scaled appropriately.13–16)
The adsorption layer should aﬀect the ﬂuid ﬂow near thecolloid particle to modify the Stokes law. Light scatteringexperiments have revealed that the diﬀusion constant of acolloid particle becomes vanishingly small when the criticalpoint is approached.17–21) To be more precise, in Ref. 19, theratio = diverges as jT  Tc j , where Tc and  are thecritical temperature and the critical exponent characterizing
the divergence of the bulk correlation length c , respectively.
This result may be explained by the divergence of the‘‘eﬀective radius’’ of the particle, r0 þ c ; the adsorption
layer with the thickness c ampliﬁes the drag force.19) In

Refs. 20 and 21, the authors estimated dependence of thedrag coeﬃcient on the bulk correlation length without
solving the hydrodynamic equation explicitly. Recently,
Nakamura et al. have developed a microscopic theory toinvestigate the drag coeﬃcient of a rigid spherical particle
in a binary hard-sphere solvent in which two solvent
species have diﬀerent radii.22) They found that the dragcoeﬃcient grows as the molar fraction of the cosolvent(larger solvent spheres) increases, which stems from thediﬀerence between the radial distribution functions of thetwo solvent species near the particle surface. Camley and
Brown have examined the drag coeﬃcient of a circularinclusion in a two-dimensional ﬂuid membrane with anonconserved scalar order parameter surrounded by a threedimensional ﬂuid.23)
In this paper, we calculate the drag coeﬃcient of a rigid
spherical particle in an incompressible near-critical binary
ﬂuid mixture, which is not extremely close to the criticalpoint. A weak preferential attraction is assumed between thesurface and one of the ﬂuid components, and the diﬀerencein the viscosity between the two components is neglected.
Still, we can show that the adsorption layer near the surface
aﬀects the drag coeﬃcient. We formulate our problem in
Sect. 2 and explicitly solve the equations in a perturbative
way in Sect. 3, with some details relegated to Appendices. In
Sect. 4, up to the lowest order of strength of the preferential
attraction, we calculate the drag coeﬃcient, the expressionof which includes the exponential integral. Section 5 isdevoted to discussion.
2.

Formulation

Our system is composed of a rigid spherical particle
immersed in a two-component ﬂuid with a constant andhomogeneous temperature T . Let A be the mass density of
one component A and B that of the other component B.
They depend on the spatial position, represented by r, andvanish inside the particle. The region inside the particle isdenoted by C, its surface by @C, and the region outside theparticle by Ce . The whole region of our system is thus
represented by C [ Ce [ @C. We deﬁne the concentration
diﬀerence ’ as A  B and the total mass density tot as
A þ B .

084003-1#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

2.1 Statics
Let us consider the equilibrium proﬁle of ’ with theparticle ﬁxed. We set the polar coordinate system ðr; ; Þ sothat the origin coincides with the particle center. We write
er , e , and e for the unit vectors along the respective
coordinate curves. The unit vector along the polar axis
is denoted by ez ; we have er  ez ¼ cos  and e  ez ¼
 sin .
We assume the ﬂuid to be almost incompressible for thepresent. Its free energy density in the bulk can be separated
into the term free from the gradient, denoted by f ðtot ; ’Þ,and the term proportional to the square gradient of ’. The
free-energy functional of our system, F , can be considered
to be given by the sum of the volume integral of these terms
and the surface contribution, which is assumed to be givenby the surface integral of the potential fs determined by ’ at
the surface.24) Thus, we have


Z
1
2
dr f ðtot ðrÞ; ’ðrÞÞ þ Mjr’ðrÞj
F¼
2
Ce
Z
þ
dS fs ð’ðrÞÞ;
ð2:1Þ
@C

where M is a positive constant. Far from the particle, we canassume A and B to reach their respective plateau values,
which can be realized in experiments. Thus, a constant, ’1 ,can be deﬁned so that it satisﬁes’ðrÞ ! ’1

as jrj ! 1:

ð2:2ÞThe most probable ’ is considered to be macroscopicallyobserved in the static state, according to the mean-ﬁeld
approximation. In the limit of incompressibility, we can drop
the constant tot from the variables of f . The most probable ’
makes
Z
F
dr ’ðrÞ
ð2:3Þ
Ce

stationary with respect to ’, whereis the undetermined
multiplier. Below, the prime indicates the derivative withrespect to the variable, e.g., f 0 ð’Þ ¼ df ð’Þ=ðd’Þ. The
stationary condition gives0f ð’ðrÞÞ  M’ðrÞ ¼

e

in C ;

ð2:4Þ

R. OKAMOTO et al.

where implies the inﬁnitesimal change. Taking the limitr ! 1 in Eq. (2.4), we have ¼ f 0 ð’1 Þ.
We introduce a dimensionless correlation length in thebulk,sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
c
1
M
¼
:
ð2:8Þ
c 
r0 r0 f 00 ð’1 ÞHere, f 00 ð’1 Þ is positive for the thermodynamic stability. Linearizing Eq. (2.4) by regarding f 0 ð’Þ asþ
f 00 ð’1 Þð’  ’1 Þ approximately, we can obtain’ðrÞ ¼ ’1 þ

ð2:5Þwhere @Cþ implies that the derivative is evaluated justoutside the particle. For simplicity, we further assume thepotential function to be a linear function in ’,
fs ð’Þ ¼ h0  h’;ð2:6Þ

where h0 and h are assumed to be constant. The right-handside (rhs) of Eq. (2.5) thus equals h. Let us write A forthe chemical potential of the component A per unit mass and
B for that of the other. The functional derivative of F withrespect to A gives A . We thus ﬁnd Eq. (2.4), or , equal tohalf of the chemical potential diﬀerence, ð A  B Þ=2,
considering that we haveA A þ

B B ¼ð Aþ

B Þ tot2þð A  BÞ ’
; ð2:7Þ2

ð2:9Þ

where   r=r0 ¼ jrj=r0 is a dimensionless radial length.
The approximation would be valid when
jhf 000 ð’1 Þjc =M  f 00 ð’1 Þ;ð2:10Þ

considering that j’  ’1 j would not be much larger than the
product of the correlation length and jer  r’j at r ¼ r0 , i.e.,than jhjc =M.
2.2 Dynamics
In the ﬂuid dynamics out of the global equilibrium, we
assume the local equilibrium to regard Eq. (2.1) as the sumof the local free energies. As in the usual thermodynamics,
we can discuss an inﬁnitesimal reversible process of amacroscopically inﬁnitesimal comoving ﬂuid region by
assuming the environments with the temperature T , thepressure tensor ðrÞ, and the chemical potentials. We neednot consider a change in the temperature in this process
because we will consider the isothermal dynamics. The
tensor  amounts to the reversible part of the pressure tensor
in the dynamics.
We write 1 for the isotropic tensor. After the same
calculations that lead to the stress tensor in the hydrodynamic model of near-critical ﬂuids (the so-called
model H),25,26) we can ﬁnd ¼ p1 þ posm 1 þ grad ;ð2:11Þ

where p comes from the inﬁnitesimal change in tot and canbe regarded as dependent only on r in the steady state in thelimit of incompressibility, and we use

the left-hand side (lhs) of which is the functional derivativeof F with respect to ’ðrÞ, andMer  r’ðrÞ ¼ fs0 ð’Þ at @Cþ ;hr0 eð1Þ= c;Mð1 þ 1c Þ

posm  ’f 0 ð’Þ  f ð’Þ
and


1
2
jr’j þ ’’ 1 þ Mr’r’:
grad  M2

ð2:12Þ



ð2:13Þ

See Appendix A for details. The chemical potential
diﬀerence is found to be^ ðrÞ ¼ f 0 ð’ðrÞÞ  M’ðrÞ;ð2:14Þ

which satisﬁes the equation,’r ^ ¼ rposm þ r  grad ;ð2:15Þ

like the Gibbs–Duhem relation. We also ﬁnd Eq. (2.5) to bevalid in the inﬁnitesimal change we now consider. We
assume the local equilibrium even just outside the particle,
which means that the reversible part of the force exerted on
the particle by the ﬂuid is given by

084003-2#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

R. OKAMOTO et al.

Z

@Cþ

dSðp1 þ posm 1 þ grad Þ  er :

ð2:16Þ

The dissipative dynamics can be formulated in termsof the irreversible thermodynamics.27) As discussed inAppendix A, we assume that the dissipation occurs throughthe viscosity and the diﬀusive ﬂux between the twocomponents. Thus, Eq. (2.5) is valid in the dissipative
dynamics. Let v denote the velocity ﬁeld, and the dissipative
stress tensor is given by 2E, where the rate-of-strain tensor
E is deﬁned asE

1
frv þ ðrvÞT g;2

ð2:17Þ

with the superscript ‘‘T’’ indicating the transposition. At thegiven temperature T , the viscosity far from the particle has arenormalized value, , which depends weakly on c .26) We
assume for simplicity that the viscosity is given by  over
the region Ce , irrespective of ’. Thus, we use the steady
Stokes approximation to have
0 ¼ v  r’  L ^ ;ð2:18Þ

0 ¼ rðp þ posm Þ  r  grad þ v;

ð2:19Þ

where the Onsager coeﬃcient L is assumed to be a positive
constant, and the second term on the rhs of Eq. (2.18)
represents the divergence of the diﬀusive ﬂux, as is foundfrom Eq. (A13). We assume the ﬂuid to be incompressible
to have
0 ¼ div v:

ð2:20Þ

In the equilibrium state discussed in Sect. 2.1, ^ ðrÞ equals
the constant , and thus pðrÞ is a constant, referred to as pð0Þ
below, because of Eqs. (2.15) and (2.19).
3.

Calculation

We introduce a small dimensionless parameter for the
perturbation scheme formulated below. In a quiescent ﬂuid,
some external force is assumed to cause a slow translationalmotion of the particle. We consider the moment when theparticle center coincides with the origin of the coordinate
system introduced in Sect. 2.1. Let U be a nonzero constantscalar with the dimension of the speed, and we assume theparticle to move with the velocity Uez at the moment. It is
assumed that, far from the particle, the ﬂuid tends to be staticand the concentration diﬀerence is still given by ’1 . Thus,
we havevðrÞ ! 0and^ ðrÞ !

as r ! 1:

ð3:1Þ

and

er  r ^ ! 0

at r ! r0 þ;

ð3:2Þ

where r ! r0 þ means that r approaches r0 with r > r0
being kept.
We expand the ﬁelds as’ðrÞ ¼ ’ð0Þ ðrÞ þ ’ð1Þ ðrÞ þ oð Þ;vðrÞ ¼ vð1Þ ðrÞ þ oð Þ;^ ðrÞ ¼ þ ^ ð1Þ ðrÞ þ oð Þ:
ð3:3Þ
On each rhs of Eq. (3.3), the ﬁeld with the superscript‘‘(1)’’ is deﬁned so that it becomes proportional to after
being multiplied by , and oð Þ implies higher order terms.
The rhs of Eq. (2.9) should be referred to as ’ð0Þ ðrÞ. We
deﬁneðÞ  þ c ð1Þ= c
e
;
1þ c

ð3:4Þ

0

so as to have ’ð0Þ ðrÞ ¼ hðÞ=ðM2 Þ with  being r=r0 .
Using Eqs. (2.4) and (2.12)–(2.14), we have
^ ð1Þ ¼ f 00 ð’ð0Þ Þ’ð1Þ  M’ð1Þ ;ð3:5Þ

00 ð0Þ ð0Þ ð1Þ
pð1Þ
osm ¼ f ð’ Þ’ ’ ;ð0Þð1Þ
 r’ð1Þ þ ’ð1Þ ’ð0Þ þ ’ð0Þ ’ð1Þ Þ1grad ¼ Mðr’
ð0Þ
ð1Þ
ð1Þ
ð0Þ

þ Mðr’ r’

þ r’ r’ Þ:

ð3:6Þ
ð3:7Þ

From Eqs. (2.18) and (2.19), we havevð1Þ  r’ð0Þ ¼ L ^ ð1Þ ;
ð1Þ

ð0Þ

ð1Þ

ð3:8Þ
ð1Þ

ð3:9Þ
rp ¼ ’ r ^ þ v :As mentioned before, we still have Eqs. (2.2) and (2.5) inthe dynamics. Thus, also noting Eqs. (3.1) and (3.2), wehave’ð1Þ ; pð1Þ ; ^ ð1Þ ! 0

and vð1Þ ! 0

as r ! 1;ð3:10Þ

while
er  r’ð1Þ ; er  r ^ ð1Þ ! 0

andvð1Þ ! Uez

as r ! r0 þ:

ð3:11Þ
We expand the ﬁelds of the order of with respect tothe vector spherical harmonics. For the integers j and msatisfying j  0 and j  m  j, they are28)
P jm ð; Þ ¼ er Y jm ð; Þ;B jm ð; Þ ¼ N 1j ðe @ þ e cosec @ ÞY jm ð; Þ;C jm ð; Þ ¼ N 1
ð3:12Þ
j ðe cosec @  e @ ÞY jm ð; Þ;pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
where N j  jð j þ 1Þ is used, and Y jm denotes
spherical
pthe
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ
harmonics. In particular, we have Y10 ð; Þ ¼ 3=ð4Þ cos .
Although P jm is deﬁned for j  0, B jm and C jm are deﬁned
only for j  1. Let ðr; ; Þ be the components of r. We thendeﬁneX
^ ð1Þ ðrÞ ¼
Q jm ðrÞY jm ð; Þ;ð3:13Þ
jm

Far from the particle, irrespective of , we ﬁnd grad tovanish and posm to remain the same, and we can assume p toremain pð0Þ . The no-slip boundary condition is imposed at
the particle surface. The radial component of the diﬀusive
ﬂux should vanish at the moving surface, as mentioned in
Eq. (A11). Thus, we havevðrÞ ! UezpðrÞ ¼ pð0Þ þ pð1Þ ðrÞ þ oð Þ;

ð1Þ

p ðrÞ ¼X

p jm ðrÞY jm ð; Þ;ð3:14Þ

jm

vð1Þ ðrÞ ¼

X

fR jm ðrÞP jm ð; Þ þ T jm ðrÞB jm ð; Þjm

ð3:15Þ
þ S jm ðrÞC jm ð; Þg;
where Q jm ; p jm ; . . . ; S jm are the r-dependent expansion
coeﬃcients. We assume T00 and S00 to be equal to zero,and the sum is over all the values of j and m consideredcjust above Eq. (3.12). Because Y jm
¼ ð1Þm Y j;m , with thesuperscript ‘‘c’’ indicating the complex conjugate, the coeﬃcients should share the corresponding property, suchas084003-3#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

Rcjm ¼ ð1Þm R j;m ;ð3:16Þ

so that the ﬁelds are real.
After performing some algebra shown in Appendix B, we
ﬁnd S jm to vanish for any ð j; mÞ, andp jm ¼ R jm ¼ T jm ¼ Q jm ¼ 0

except for ð j; mÞ ¼ ð1; 0Þ:ð3:17Þ
Thus, we can concentrate on the mode of ð1; 0Þ, andEqs. (3.13)–(3.15) are respectively reduced to^ ð1Þ ¼ Q10 Y10 ;pð1Þ ¼ p10 Y10 ;ð3:18Þ

andvð1Þ ¼ R10 P10 þ T10 B10 :ð3:19Þ

Equation (2.20) relates T10 with R10 , as shown by Eq. (B1).
We write @r for the diﬀerentiation with respect to r, and @2r for@r @r . Substituting Eqs. (3.18) and (3.19) into Eq. (3.8) gives
Lð@2r þ 2r 1 @r  2r 2 ÞQ10 ¼ R10 @r ’ð0Þ :ð3:20Þ

The rhs of Eq. (3.9) is irrotational, which is combined with
(3.18) and (3.19) to giveðr 2 @4r þ 8r@3r þ 8@2r  8r 1 @r ÞR10 ¼ 2Q10 @r ’ð0Þ ;ð3:21Þ

with the aid of Eq. (B1). Equations (3.20) and (3.21) arerespectively the mode of ð1; 0Þ of Eq. (B3) and that ofEq. (B8). Equations (3.10) and (3.11) give
p10 ; R10 ; Q10 ! 0 as r ! 1;
ð3:22Þ
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
R10 ! 4=3U and @r Q10 ; @r R10 ! 0 as r ! r0 þ;
ð3:23Þ
as mentioned in the third paragraph of Appendix B.
From Eqs. (3.20), (3.22), and (3.23), we ﬁnd that Q10
vanishes if h vanishes. We introduce dimensionless functions, QðÞ  3MLQ10 ðrÞ=ðUhr02 Þ and RðÞ  R10 ðrÞ=U,
where we use   r=r0 as before. We also introduce a
dimensionless factorhr02
pﬃﬃﬃﬃﬃﬃﬃﬃﬃ :3M 5L

ð3:24Þ

If h vanishes, we ﬁnd RðÞ to be given byrﬃﬃﬃﬃﬃﬃ

4
1
3
 3þ
Rh¼0 ðÞ 
ð3:25Þ
32
2
from Eqs. (3.21)–(3.23). In general, as shown inAppendix B, applying the method of variation of parameters
to Eqs. (3.20)–(3.23), we can obtainZ1QðÞ ¼d ðÞQ ð; ÞRðÞ;
ð3:26Þ
1
Z1
RðÞ ¼ Rh¼0 ðÞ þ  2d ðÞR ð; ÞQðÞ: ð3:27Þ
1

Here,



Q ð; Þ andR ð; Þ fðÞ2 =2g þ 2
2fðÞ =2g þ 2

if 1   < 
if   Rcom ð; Þ þ  2 3  51
Rcom ð; Þ þ 2  3  5 1

ð3:28Þ

if 1   < ;
if   
ð3:29Þ

R. OKAMOTO et al.

where the common term isRcom ð; Þ 

3  5 2 5ð3 2  1Þ
þ
:
2 3
2 3 3

ð3:30Þ

The lowest-order correction due to a nonzero h value isgiven by replacing R by Rh¼0 on the rhs of Eq. (3.26) andsubstituting the resultant approximate solution of Q into therhs of (3.27). Thus, up to this order, we haveZ1QðÞ ¼d ðÞQ ð; ÞRh¼0 ðÞ;
ð3:31Þ
1
Z1
RðÞ ¼ Rh¼0 ðÞ þ  2d ðÞR ð; Þ1
Z1d 0 ð 0 ÞQ ð;  0 ÞRh¼0 ð 0 Þ:
ð3:32Þ
1

For later convenience, we write pð1Þ and vð1Þ in terms of Q
and R. Let r have the components ðr0 ; ; Þ in the polar
coordinate system. We can use Eqs. (B1) and (B6) torewrite the second equation of Eq. (3.18) as 1 3 3
ð1Þ
2 2
 @ þ 3 @ þ 2@ RðÞp ðrÞ ¼ U
r0 2
hr02 ð0Þ’ ðrÞQðÞ Y10 ð; Þ;ð3:33Þ
3ML
and use Eq. (B1) to rewrite Eq. (3.19) as@Y10 ð; Þ @ þ 1 RðÞ:
vð1Þ ðrÞ ¼ U er Y10 ð; Þ þ e@2
ð3:34Þ
4.

Results
For simplicity, we assume the Gaussian modelf ð’Þ ¼að’  ’1 Þ2 þ ð’  ’1 Þ;
2

ð4:1Þ

where a is a positive constant. This constant, being thereciprocal susceptibility, can be assumed to be proportional
to the temperature measured from the critical point. From
Eqs. (3.5) and (3.18), we haveða  MÞ’ð1Þ ðrÞ ¼ Q10 Y10 ;ð4:2Þ

and G10 ðrÞ can be deﬁned so as to satisfy’ð1Þ ðrÞ ¼ G10 ðrÞY10 ð; Þ:ð4:3Þ

From these two equations, we can obtainfa  Mð@2r þ 2r 1 @r  2r 2 ÞgG10 ðrÞ ¼ Q10 ðrÞ:ð4:4Þ
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Here we have c ¼ M=a=r0 from Eq. (2.8). Introducing adimensionless variable ~  = c ¼ r=ðr0 c Þ, we ﬁnd linear
independent solutions of the homogeneous equation associated with Eq. (4.4) to be given byrﬃﬃﬃﬃ

pﬃﬃﬃ 1 2
1~cosh ~  sinh ~ ;I3=2 ðÞ=
ð4:5Þ
~ ¼
~ 
~
rﬃﬃﬃﬃ

pﬃﬃﬃ 1 
1 ~
~
ð4:6Þ
1þ
K3=2 ðÞ=
~ ¼
e ;
~ 2
~
where I3=2 and K3=2 represent modiﬁed Bessel functions.
Equations (3.10) and (3.11) lead to

084003-4

G10 ðrÞ ! 0@r G10 ðrÞ ! 0

as r ! 1;
as r ! r0 þ:

ð4:7Þ
ð4:8Þ

#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

Applying the method of variation of parameters, we canobtain G10 so that it satisﬁes these boundary conditions. Bymeans of ~  = c and ~  = c , we ﬁrst deﬁne~ 3=2 ðÞ~ ;~
Gcom ð; Þ  wð c ÞK3=2 ðÞK

ð4:9Þ

where
wð c Þ 

0
1
I3=2
ð 1
c Þ  c I3=2 ð c Þ=2;
1
K03=2 ð 1c Þ  c K3=2 ð c Þ=2

and then deﬁne~ 3=2 ðÞ~ ~
Gcom ð; Þ þ K3=2 ðÞIG ð; Þ 
~ 3=2 ðÞ~ ~
Gcom ð; Þ þ I3=2 ðÞK

ð4:10Þ

if 1   < :
if   
ð4:11Þ

After performing some algebra, we can obtainZ1
pﬃﬃﬃ
Uhr04 cd  G ð; ÞQðÞ:
G10 ðrÞ ¼
pﬃﬃﬃ
3M 2 L  1

ð4:12Þ

We can write Drag ez for the drag force exerted on theparticle. Noting Eq. (2.16), we obtainZdS er  fðp þ posm Þ1 þ grad  2Eg  ez :Drag ¼ @Cþð4:13Þ
The unperturbed terms in Eq. (3.3) do not contribute to
the drag force because they are only dependent on r.
Equations (3.6) and (3.7) yield
ð0Þ
pð1Þ
osm ¼ a’ G10 Y10 ;ð4:14Þ

ð0Þ 0
ð1Þ
@r G10 þ G10 ’ð0Þ ÞY10 þ ’ð0Þ ’ð1Þ g1
grad ¼ Mfð’0
þ M’ð0Þ ½er rðG10 Y10 Þ þ frðG10 Y10 Þger ;ð4:15Þ
ð1Þ

where ’ can be rewritten as ða’
Eq. (3.5), and we can use
rðG10 Y10 Þ ¼

ð1Þ

^

ð1Þ

Þ=M because ofdG10 ðrÞ
G10 @Y10
Y10 er þ
e :
dr
r @

ð4:16Þ

As mentioned in Eq. (2.6) of Ref. 28, we use Eq. (B1) toobtain2er  E

ð4:18Þ
Because of Eq. (2.5), ’ equals h=M in the surface integral of Eq. (4.13). Thus, using Eqs. (3.23), (3.33), (4.8), and(4.14)–(4.18), we can replace the integrand of Eq. (4.13) by1 1 3 3
2 2
r @r þ 3r @r R10 Y10 cos 
r
2
2

ð4:20Þ
where j¼1 implies that the derivative should be evaluated at
 ¼ 1, and2þcc

ð4:21Þ

:pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
We deﬁne XðÞ so that Rh¼0 ðÞ ¼ 4=3XðÞ. Substituting
Eqs. (3.31), (3.32), and (4.12) into Eq. (4.20), we can obtain
the drag coeﬃcient   Drag =ð UÞ as
h2 r04
 ¼ 6r0 1 þðÞ
ð4:22Þ
27M 2 L c
up to the order of  2 / h2 , where

Z1 
3
1

d
ð cÞ ¼
ðÞ
 3
1
Z1
d ðÞQ ð; ÞXðÞ  2 c Wð c Þ
Z 11
Z1
pﬃﬃﬃ
d G ð1; Þd ðÞQ ð; ÞXðÞ11

ð4:23Þ
represents how the deviation from the Stokes law depends
on c . Equation (4.22) shows that the drag coeﬃcient is
independent of ’1 . The independence of Eq. (4.22) from thesign of h can be expected because neither component hasa special bulk property in our formulation. Equation (4.22)
would have a term proportional to h and would depend on its
sign if the viscosity were assumed to depend on ’.
By means of Mathematica version 8 (Wolfram Research),29)
we can calculate the integrals of Eq. (4.23) to obtain1ð cÞ ¼½A1 ð c Þ þ 256A2 ð c Þe2= c Eið2= c Þ1280ð1 þ c Þ2
 10A3 ð c Þe1= c Eið1= c Þþ 5fA4 ð c Þe1= c Eið1= c Þg2 ;where we use the exponential integral
Zx
et
dt ;
EiðxÞ t
1

ð4:24Þ

ð4:25Þ

and introduce
A1 ðzÞ  3840z3  6156z2 þ 168z þ 568
þ 32z1 þ 161z2  10z3 þ 5z4 ;ð4:26Þ

3

ð4:27Þ

A2 ðzÞ  90z  45z þ 2z ;2A3 ðzÞ  1152z  24  120z
2

þ 22z

3

4

þz

1

 18z

2

5

z ;

ð4:28Þ

A4 ðzÞ  12z1 þ z3
ð4:29Þ
to write Eq. (4.24) concisely. As c increases, the ﬁrst term
A1 is the most relevant in the braces of Eq. (4.24), making
ð c Þ almost proportional to c . Using the asymptotic
expansion
ex EiðxÞ

ðr 2 @2r þ 2r@r Þ’ð0Þ Y10 cos þ MG10 r


 @Y
 2 2
10
1
r
sin  : ð4:19Þ
hG10 þ r @r R102@
Thus, we can rewrite the rhs of Eq. (4.13) as)rﬃﬃﬃﬃ(
 1 3
hWð c Þ2
@ þ 2@ RðÞ
G10 ðr0 Þ ; 2 r0 U

3
2 
U
¼1

1

Wð c Þ  2 þ

ð4:17Þ
¼ ðer rot vð1Þ þ 2@r vð1Þ Þpﬃﬃﬃ 1¼ 2@r R10 P10 ð; Þ þ 2 r@r þ 1 @r R10 B10 ð; Þ :2
ð0Þ 0

R. OKAMOTO et al.

1
X

ð1Þn

n¼1

ðn  1Þ!
xn

as x ! 1ð4:30Þ

in Eq. (4.24), we obtain27 6cð8  111 c þ 1231 2c Þ
8ð1 þ c Þ2
as c ! 0þ;ð4:31Þ
which tells thattends to 0 as c ! 0þ, i.e., as c
approaches 0 along the real axis from the positive side. It is
to be noted that the correlation length is not assumed to be
microscopic in Eq. (4.1). In Fig. 1, we show that Eq. (4.24)084003-5ð cÞ#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

R. OKAMOTO et al.

2

0

10

δ

1.5
-2

10

1
-4

10

0.5
-6

10

0
-2

-8

-1

0

1

2

10

Fig. 2. The velocity ﬁeld for h ¼ 0 is drawn in a region with x  0 on thexz plane when viewed from the frame comoving with the particle. The
dashed curve represents the cross section of the particle surface.

-10

10

2

3 4 562

0.1

3 4 56ζc

2

1

3 4 56

10

Fig. 1. Plot of against c . The solid and dashed curves respectively
represent Eqs. (4.24) and (4.31).

vz(1)/U-1

0.0

is positive and increases with c . According to Eq. (4.22),
this means that the drag coeﬃcient increases from the value
given by the Stokes law as the strength of the preferential
attraction jhj increases and/or as the bulk correlation length
increases.
5.

-0.2

-0.4

-0.6

Discussion

Following the framework of the irreversible thermodynamics, we formulate the ﬂuid dynamics to examine the dragcoeﬃcient of a rigid particle in a near-critical binary ﬂuid
mixture. We take into account the preferential attractionof one of the ﬂuid components to the particle surface,
represented by Eq. (2.6), and the dynamical coupling
between the concentration heterogeneity and the hydrodynamic ﬂow, represented by Eqs. (2.13), (2.18), and (2.19).
The mean-ﬁeld approximation we use holds when thecorrelation length, c ¼ r0 c , is so small that higher-orderterms, such as the ’4 term, are negligible in Eq. (2.1).
We use the Gaussian model by assuming Eq. (4.1), which
is substituted into Eq. (2.4) to yield Eq. (2.9) without
approximation. We assume a weak preferential attraction
between the surface and one of the ﬂuid components( 2  1), and assume the viscosity  to be independent of
the composition, ’. Still, the drag coeﬃcient deviates fromthe value of the Stokes law, as shown by Eq. (4.22) with
Eq. (4.24). We should note the weak dependence of  onc , as mentioned above Eq. (2.18). When the correlation
length is small, we ﬁnd from Eq. (4.31) that the deviation of
the drag coeﬃcient from the value of the Stokes law is
proportional to 6c . This dependence is much steeper than theprevious estimate of the 2c dependence due to the surface
energy eﬀect, which dependence is obtained without solving
the hydrodynamic equation explicitly in Ref. 21. The slopeof the plots in Fig. 1 gradually decreases from six towards
the unity as c increases.
The velocity ﬁeld for h ¼ 0 can be calculated by means ofEqs. (3.25) and (3.34). This velocity ﬁeld is well known andleads to the Stokes law. Subtracting the particle velocity
from the ﬁeld, we draw the resultant ﬁeld on the xz plane in
Fig. 2 by means of Mathematica. The ﬁgure thus shows thevelocity ﬁeld for h ¼ 0 viewed from the frame comoving

-0.8

-1.0
0

2

4

6

x/r0 8

10

Fig. 3. Plot of vð1Þz =U  1 along the x-axis for h ¼ 0.

with the particle. There, an arrow represents the velocity
direction with its length proportional to the magnitude.
Although similar ﬁgures are given below in Figs. 4 and 6,
each ﬁgure has its constant of proportion, and thus acomparison of the lengths of two arrows in diﬀerent ﬁgures
is meaningless. In harmony with Fig. 2, we plot vð1Þz =U  1for h ¼ 0 along the x-axis on the positive side in Fig. 3. It
approaches 1 far from the particle and vanishes at theparticle surface, because of the last equations of Eqs. (3.10)
and (3.11).
As we obtained Eq. (4.24) from Eq. (4.23), we cancalculate the double integral on the rhs of Eq. (3.32). The
result, not shown here, contains the exponential integral. By
using this result and Eq. (3.34), we can obtain the deviationof vð1Þ caused by a nonzero h value. For c ¼ 0:1, thedeviation divided by U 2 on the xz plane is used for Fig. 4,
and its z-component on the x-axis is plotted in Fig. 5.
Comparing Fig. 2 with Fig. 4, we ﬁnd in the comoving
frame the preferential attraction brings about a decrease in
the magnitude of the ﬂuid velocity around the particle. Thistendency is also shown by Figs. 3 and 5, where vð1Þz =U  1for h ¼ 0 and the z-component of the deviation ﬁeld do notshare the same sign. Noting that the quotient of the deviation
divided by U 2 / h2 is plotted in Fig. 5, we ﬁnd thedecrease to be larger as h2 is larger. This means that theparticle drags the ﬂuid more signiﬁcantly then.

084003-6#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

4

4

3

3

2

2

1

1

0
-4

-2

0

2

R. OKAMOTO et al.

0
-4

4

Fig. 4. For c ¼ 0:1, the deviation of vð1Þ from vð1Þ for h ¼ 0 isrepresented by the arrows in a region with x  0 on the xz plane. We
calculate the deviation by using the second term on the rhs of Eq. (3.32),
and use the result divided by U 2 for the arrows. In a darker-gray region,
the dimensionless diﬀerence ð’ð0Þ  ’1 ÞM=ðhr0 Þ is larger. The white
semicircular region surrounded by the dark adsorption layer represents the
cross section of the particle.

-2

0

2

Fig. 6. Same as in Fig. 4, except forc ¼ 1.

-2

2

Δvz(1) /(Uλ )

8x10

4

-5

1.4x106

2

Δvz(1)/(Uλ )

4
1.2
1.0

2

0.8
0

0.6

0
0.4

2

4

6

x/r0 8102
Fig. 7. For c ¼ 1, the quotient of the deviation of vð1Þz divided by U isplotted by means of the solid curve, as in Fig. 5. The dashed curve
represents ð’ð0Þ  ’1 ÞM=ð10hr0 Þ for c ¼ 1. Here, the divisor 10 isintroduced for convenience of the graphical presentation.

0.2
0.0
0

2

4

6

x/r0 8

10

Fig. 5. For c ¼ 0:1, the deviation of vzð1Þ from its value for h ¼ 0 is
divided by U 2 , and the quotient is plotted along the x-axis by means of thesolid curve. The symbol  in the axis label means the deviation. The dashedcurve is drawn for reference to represent ð’ð0Þ  ’1 ÞM=ð104 hr0 Þ forc ¼ 0:1.

The second term on the rhs of Eq. (2.9) is proportional toh, and the function ’ð1Þ ðrÞ, induced by the imposed ﬂow,
also has a term proportional to h. This term can be calculated
from Eqs. (3.31), (4.3), and (4.12). The calculation can be
performed by Mathematica, like that of the integral in
Eq. (3.32), although the result is not shown here. The graylevel in Fig. 4 simply represents ð’ð0Þ  ’1 ÞM=hr0 ; itsvalue, divided by 104 for convenience, is plotted along thex-axis by means of the dashed curve in Fig. 5.
Similar ﬁgures for c ¼ 1 are shown in Figs. 6 and 7. The
solid curve of the latter ﬁgure takes larger values than that ofFig. 5. This implies that the particle should drag the ﬂuid
more signiﬁcantly not only as the preferential attraction
becomes stronger but also as the correlation length increases.
Then, in other words, the particle becomes apparently sticky
and eﬀectively larger. This is also deduced from Eq. (4.22)
and Fig. 1. The 6c dependence in Eq. (4.31) tempts us toassociate the factor h with 3c , representing the volume ofthe spherical region with the eﬀective radius. This steep
dependence cannot be explained simply by replacing r0 withr0 þ c in the Stokes law.
Finally, we make some remarks. (i) Near the critical point,
a mixture solvent is very sensitive to selective impurities
such as ions which interact with each solvent species
asymmetrically.30–32) Theoretical works have shown thateven a small amount of salt can play a signiﬁcant role in thestatic behavior of the adsorption layer.33–35) The ion eﬀect on
the drag coeﬃcient would be signiﬁcant in experiments.
(ii) Several groups have devised numerical schemes for
simulating the dynamics of colloid particles in a binary ﬂuid
mixture in a two-phase state (T < Tc ).36–39) Araki andTanaka have examined the dynamic eﬀect on the adsorptioninduced force between colloid particles.38,39) In Refs. 40 and
41, Cates and coworkers have utilized the scheme ofRefs. 36 and 37 to investigate bicontinuous interfaciallyjammed emulsion gels (‘‘bijels’’), in which colloid particles
are trapped at the liquid–liquid interface of the phaseseparated mixture solvent. In the present paper, we haveassumed that the ﬂuid is in a one-phase state. However, ourformulation can also be applied to these problems if the freeenergy density f ð’Þ is chosen so that it describes the phaseseparated ﬂuids (e.g., ’4 model). (iii) When the mixture
is very close to the critical point, the eﬀect of critical

084003-7#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

ﬂuctuations would be important. Because of the criticalﬂuctuations, the concentration proﬁle near the surface
becomes universal, as mentioned in Sect. 1. The ﬂuctuationeﬀect (or the renormalization eﬀect) on the drag coeﬃcient
can be taken into account by replacing the free-energyfunctional Eq. (2.1) by the ‘‘renormalized local functional’’
in Ref. 16.
Acknowledgement
The authors thank Professor A. Onuki for helpfulcomments on the critical dynamics.
Appendix A:

Colloid Dynamics in the Model H

The irreversible thermodynamics, based on the local
equilibrium and the conservation laws, provides us with aframework of phenomenological dynamics.27) The model
H dynamics for a bulk ﬂuid mixture with a constant andhomogeneous temperature is also formulated in terms of thisframework.26) Below we show that the formulation stillworks well even if colloid particles are immersed in the
mixture. We assume a single particle in the following
discussion for simplicity, but the extension to a system
containing many particles is straightforward.
To consider an inﬁnitesimal reversible change of a local
region, we deﬁne r as a small change in the positional
vector. The displacement is represented by r 7 ! r0  r þ r,
where r should be homogeneous in C. We can assume theno-slip condition at @C because it is assumed in the
dynamics we consider. The volume element dr is changed
into dr0 by the accompanying deformation, and theirrelationship is given by
dr ¼ ð1  r  rÞ dr0 :ðA1Þ

We can regard the local region as an open system so that theconcentration diﬀerence ’ is changed not only by the
deformation above but also by the ﬂux into or from theparticle bath although the ﬂux is assumed not to change tot .
The total change in ’ is thus represented by
0

0

0

0

’ðrÞ 7 ! ’ ðr Þ  ’ðrÞ þ ’d ðr Þ þ ’p ðr Þ;0’d ðr Þ ¼ ’ðrÞr  r:

surface. We thus ﬁnd Eqs. (2.5) and (2.14) to hold under the
local equilibrium. Also, considering the contribution due tothe inﬁnitesimal change of tot , we can obtain Eq. (2.11).
We have derived Eq. (A4) to discuss the dynamics, but
we can also use it to discuss the statics. Requiring Eq. (A4)
to vanish for arbitrary ’p , we can derive Eqs. (2.4) and(2.5). When ’ satisﬁes them, requiring Eq. (A4) to vanish
for an arbitrary displacement with the particle ﬁxed, we ﬁndr  osm to vanish in the equilibrium ﬂuid containing theparticle ﬁxed externally, which is consistent with Eqs. (2.4)
and (2.15). This result is obtained by applying the
integration by parts to the second term in the volume
integral of Eq. (A4). We can obtain Eq. (2.16) byconsidering the inﬁnitesimal translational displacement of
the particle in the equilibrium ﬂuid after taking into accountthe contribution involving tot .
The local equilibrium is assumed in the irreversible
thermodynamics.27,42) Thus, we have
!
!
D
D
’
f~
¼  : rv þ ^ tot
tot
;
ðA6Þ
Dt tot
Dt tot
where f~ is the free-energy density and D=Dt indicates the
Lagrangian time derivative. Here we consider the dynamics
in more general than considered in the text, allowing theunsteady dynamics. If we regard r= t as the velocity, witht denoting the inﬁnitesimal change in time t, we can also
ﬁnd the rhs of Eq. (A4) to be equal to Eq. (A5) byintegrating Eq. (A6) over the whole region of our system. A
similar procedure is mentioned in Sect. 3 of Ref. 43 andbelow Eq. (2.9) of Ref. 44. Introducing the dissipative stresstensor  and the diﬀusive ﬂux jirr , we can assume theconservation laws,
 
D ’
ðA7Þ
tot
¼ r  jirr ;Dt totDv¼ r  ð þ Þ;
ðA8Þ
tot
Dt
and

ðA2Þ

where ’d represents the change due to the deformation and’p due to the ﬂux. With the aid of Eq. (A1), we haveðA3Þ

After some calculations similar to those in Appendix 6Aof Ref. 26, we can deriveZ
F¼
dr½f f 0 ð’ðrÞÞ  M’ðrÞg ’p ðrÞ  osm ðrÞ : r r
CeZdSf fs0 ð’ðrÞÞ  Mer  r’ðrÞgf ’d ðrÞ þ ’p ðrÞg;
þ
@C

ðA4Þ
where osm is deﬁned as the sum of posm 1 and grad given
by Eqs. (2.12) and (2.13), respectively. According to thethermodynamics, Eq. (A4) should be equal to the sum ofZ
drf ^ ðrÞ ’p ðrÞ  ðrÞ : r rg
ðA5Þ
Ce

and similar contributions from the colloid surface andinterior. These contributions are neglected in this paper;
neither A nor B is assumed to chemically bond to the

R. OKAMOTO et al.

dV
¼
m
dt

Z
dS n  ð þ Þ;

ðA9Þ

@Ct

where m and V respectively denote the mass and the velocityof the particle, the subscript t of Ct means the timedependence of the comoving region, and n is the unit normal
vector of @Ct directed outside the particle. Here, forsimplicity, we assume no rotation of the particle, no external
conservative force exerted on the particle, and no dissipation
at the surface. The total free energy is given byZtot 2m
jvj dr þ jVj2 :F tot ¼ F þðA10Þ
e2
2Ct
The normal component of the diﬀusive ﬂux should vanish at
the surface, i.e.,n  jirr ¼ 0

at @Cþt :ðA11Þ

Substituting Eqs. (A7)–(A9) into Eq. (A6), we integrate
the result over the whole region of our system to obtainZd
F tot ¼
drð jirr  r ^   : rvÞðA12Þ
dt
Cetwith the aid of Eq. (2.5).

084003-8#2013 The Physical Society of Japan

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 084003

Noting Eq. (2.20), we can derive the linear phenomenological laws,
jirr ¼ Lr ^and  ¼ 2E;ðA13Þ

from Eq. (A12) in terms of the irreversible thermodynamics. Positive L and  ensure that the time derivative ofF tot is not positive. Substituting Eq. (A13) into Eqs. (A7)
and (A8) yields Eqs. (2.18) and (2.19), respectively. The
second equation of Eq. (3.2) follows Eq. (A11).
Appendix B: Details in the Perturbative Calculation
Substituting Eq. (3.15) into Eq. (2.20) gives
1
2
T jm ¼ N 1
j r @r r R jm

for j  1;

ðB1Þ

2

and R00 / r , which leads to
R00 ¼ 0ðB2Þ
ð1Þ

because of the boundary condition for v in Eq. (3.11).
Substituting Eqs. (3.13) and (3.15) into Eq. (3.8) gives
Lr 2 ðr 2 @2r þ 2r@r  N 2j ÞQ jm ¼ R jm @r ’ð0ÞðB3Þfor j  0. From Eqs. (B2) and (B3), we ﬁnd Q00 ¼ 0 withthe aid of the boundary condition for ^ in Eqs. (3.10) and(3.11).
Substituting Eqs. (3.13)–(3.15) into Eq. (3.9) gives
X
0¼
fO1jm ðrÞP jm þ O2jm ðrÞB jm þ O3jm ðrÞC jm g; ðB4Þ
jm

where Oijm ’s are deﬁned below. They vanish because of
the orthogonality of the vector spherical harmonics, and wehave

R. OKAMOTO et al.

@r Q jm ! 0

as r ! r0 þ:ðB11Þ

Equation (3.10) gives
p jm ; Q jm ; R jm ; S jm ! 0 as r ! 1:ðB12Þ

We ﬁnd S jm to vanish for any ð j; mÞ because of Eq. (B7)
together with the boundary conditions shown above. We
also ﬁnd Eq. (3.17) to satisfy Eqs. (B1)–(B4) and all theboundary conditions shown above.
We can regard Eqs. (3.20) and (3.21) as inhomogeneousdiﬀerential equations with respect to Q and R, respectively.
The homogeneous equations associated with them are both
equidimensional in r.45) A set of linear independent solutions
for the former equation is fr 2 ; rg, while a set for the latter
is fr 3 ; r 1 ; 1; r 2 g. Applying the method of variation ofparameters to the inhomogeneous equations, we obtain


Z
2
c1 þ
d ðÞRðÞQðÞ ¼ 
1


Z
2
d  ðÞRðÞ ; ðB13Þþ  c2 
1where c1 and c2 are constants, and


Z
3
2
2
d1 þ d  ðÞQðÞRðÞ ¼ 
1


Z
1
2
þ
d2  5d ðÞQðÞ1Zd  1 ðÞQðÞþ d3 þ 5 2
1


Z
2
2
3
d  ðÞQðÞ ;þ  d4  ðB14Þ

1

0 ¼ O1jm ðrÞ  @r p jmþ r2 fð@r r 2 @r  ðN 2j þ 2ÞÞR jm þ 2N j T jm g  ’ð0Þ @r Q jmðB5Þ

where d1 ; . . . ; d4 are constants. After determining theconstants by means of the boundary conditions given by
Eqs. (3.22) and (3.23), we obtain Eqs. (3.26) and (3.27).

for j  0,
0 ¼ O2jm ðrÞ  N j r 1 p jmþ r 2 f2N j R jm þ ð@r r 2 @r  N 2j ÞT jm g  N j r 1 ’ð0Þ Q jmðB6Þ

for j  1, and
0 ¼ O3jm ðrÞ  r 2 ð@r r 2 @r  N 2j ÞS jm

for j  1:

ðB7Þ

From Eqs. (B2) and (B5), we ﬁnd p00 ¼ 0 because of theboundary condition for p in Eq. (3.10). Deleting p jm from
Eqs. (B5) and (B6) for j  1, we can substitute Eq. (B1)
into the resultant equation to obtainr 2 fr 4 @4r þ 8r 3 @3r þ 2ð6  N 2j Þr 2 @2r
 4N 2j r@r þ N 2j N~ 2j gR jm ¼ N 2j Q jm @r ’ð0ÞðB8Þ
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ~
for j  1, where N j  ð j  1Þð j þ 2Þ.
The boundary condition at the surface is given by
Eq. (3.11). As r ! r0 þ, we thus have S jm ! 0 for anyð j; mÞ, while R jm ! 0 and T jm ! 0 except for ð j; mÞ ¼
ð1; 0Þ. We also have Eq. (3.23), andpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
T10 ! 2 2=3U as r ! r0 þ:ðB9Þ
Substituting these two equations into Eq. (B1) gives@r R jm ! 0

as r ! r0 þ:

The condition for ^ ð1Þ in Eq. (3.11) givesðB10Þ



okamoto r@fukui.kyoto-u.ac.jp
youhei@appi.keio.ac.jpzkomura@tmu.ac.jp
1) B. V. Derjaguin and L. D. Landau: Acta Physicochim. (USSR) 14(1941) 633.
2) E. J. W. Verwey and J. Th. G. Overbeek: Theory of the Stability of
Lyophobic Colloids (Elsevier, Amsterdam, 1948) Chaps. 8–12.
3) W. B. Russel, D. A. Saville, and W. R. Schowalter: Colloidal
Dispersions (Cambridge University Press, Cambridge, U.K., 1989)
Chaps. 1–5.
4) G. G. Stokes: Trans. Cambridge Philos. Soc. 9 (1851) 8.
5) W. Sutherland: Philos. Mag. 9 (1905) 781.
6) A. Einstein: Ann. Phys. (Leipzig) 322 (1905) 549.
7) D. Beysens and D. Estève: Phys. Rev. Lett. 54 (1985) 2123.
8) B. M. Law, J.-M. Petit, and D. Beysens: Phys. Rev. E 57 (1998) 5782.
9) J.-M. Petit, B. M. Law, and D. Beysens: J. Colloid Interface Sci. 202(1998) 441.
10) D. Beysens and T. Narayanan: J. Stat. Phys. 95 (1999) 997.
11) H. Guo, T. Narayanan, M. Sztuchi, P. Schall, and G. H. Wegdam:
Phys. Rev. Lett. 100 (2008) 188303.
12) D. Bonn, J. Otwinowski, S. Sacanna, H. Guo, G. Wegdam, and P.
Schall: Phys. Rev. Lett. 103 (2009) 156101.
13) M. E. Fisher and P. G. de Gennes: C. R. Seances Acad. Sci., Ser. B 287(1978) 207.
14) C. Hertlein, L. Helden, A. Gambassi, S. Dietrich, and C. Bechinger:
Nature 451 (2008) 172.
15) A. Gambassi, A. Maciołek, C. Hertlein, U. Nellen, L. Helden, C.
Bechinger, and S. Dietrich: Phys. Rev. E 80 (2009) 061143.
y084003-9#2013 The Physical Society of Japan

J. Phys. Soc. Jpn. 82 (2013) 084003

FULL PAPERS

16) R. Okamoto and A. Onuki: J. Chem. Phys. 136 (2012) 114704.
17) Ya. A. Bal’tsevich, V. G. Martynets, and E. V. Matizen: Zh. Eksp.
Teor. Fiz. 51 (1966) 983 [Sov. Phys. JETP 24 (1967) 654].
18) V. G. Martynets and E. V. Matizen: Zh. Eksp. Teor. Fiz. 58 (1970) 430[Sov. Phys. JETP 31 (1970) 228].
19) S. P. Lee: Phys. Rev. Lett. 36 (1976) 1319.
20) K. B. Lyons, R. C. Mockler, and W. J. O’Sullivan: Phys. Rev. Lett. 30(1973) 42.
21) K. B. Lyons, R. C. Mockler, and W. J. O’Sullivan: Phys. Rev. A 10(1974) 393.
22) Y. Nakamura, A. Yoshimori, and R. Akiyama: J. Phys. Soc. Jpn. 81(2012) SA026.
23) B. A. Camley and F. L. H. Brown: Phys. Rev. E 85 (2012) 061921.
24) J. W. Cahn: J. Chem. Phys. 66 (1977) 3667.
25) P. C. Hohenberg and B. I. Halperin: Rev. Mod. Phys. 49 (1977) 435.
26) A. Onuki: Phase Transition Dynamics (Cambridge University Press,
Cambridge, U.K., 2002) Chap. 6.
27) S. R. de Groot and P. Mazur: Non-Equilibrium Thermodynamics(Dover, New York, 1984) Chap. 4.
28) Y. Fujitani: J. Phys. Soc. Jpn. 76 (2007) 064401.
29) To calculate the integral with respect to , we replaced 1= c with avariable  in the integrand, took e =ð1 þ 1 Þ, which lies in Eq. (3.4),
out of the integral sign, and added the option ‘‘Assumptions !   1&&  > 0’’.

R. OKAMOTO et al.

30) A. Onuki, R. Okamoto, and T. Araki: Bull. Chem. Soc. Jpn. 84 (2011)569.
31) A. Onuki and R. Okamoto: Curr. Opin. Colloid Interface Sci. 16(2011) 525.
32) R. Okamoto and A. Onuki: Phys. Rev. E 82 (2010) 051501.
33) R. Okamoto and A. Onuki: Phys. Rev. E 84 (2011) 051401.
34) S. Samin and Y. Tsori: Europhys. Lett. 95 (2011) 36002.
35) U. Nellen, J. Dietrich, L. Helden, S. Chodankar, K. Nygard, J. F.
van der Veen, and C. Bechinger: Soft Matter 7 (2011) 5360.
36) J.-C. Desplat, I. Pagonabarraga, and P. Bladon: Comput. Phys.
Commun. 134 (2001) 273.
37) K. Stratford, R. Adhikari, I. Pagonabarraga, and J.-C. Desplat: J. Stat.
Phys. 121 (2005) 163.
38) T. Araki and H. Tanaka: Phys. Rev. E 73 (2006) 061506.
39) T. Araki and H. Tanaka: J. Phys.: Condens. Matter 20 (2008) 072101.
40) M. E. Cates and P. S. Clegg: Soft Matter 4 (2008) 2132.
41) E. Kim, K. Stratford, R. Adhikari, and M. E. Cates: Langmuir 24(2008) 6549.
42) H. Fukagawa and Y. Fujitani: Prog. Theor. Phys. 127 (2012) 921.
43) Y. Fujitani: J. Phys. Soc. Jpn. 70 (2001) 1556.
44) Y. Fujitani: J. Phys. Soc. Jpn. 79 (2010) 074002.
45) C. Bender and S. A. Orszag: Advanced Mathematical Methods for
Scientists and Engineers (Springer, New York, 1999) Chap. 1.

084003-10#2013 The Physical Society of Japan