Phase separation in a binary mixture confined between symmetric parallel plates: Capillary condensation transition near the bulk critical point

1.8 MB·English

by  Shunsuke Yabunaka

Phase separation in a binary mixture confined between symmetric parallel plates:
Capillary condensation transition near the bulk critical point
Syunsuke Yabunaka1, Ryuichi Okamoto2, and Akira Onuki1
1Department of Physics, Kyoto University, Kyoto 606-8502
2Fukui Institute for Fundamental Chemistry, Kyoto University, Kyoto 606-8103
(Dated: December 11, 2013)
Weinvestigatephaseseparationofnear-criticalbinarymixturesbetweenparallelsymmetricwalls
inthestrongadsorption regime. Wetakeintoaccounttherenormalization effectduetothecritical
fluctuationsusingtherecentlocalfunctionaltheory[J.Chem. Phys. 136,114704(2012)]. Instatics,
avanderWaalsloopisobtainedintherelationbetweentheaverageorderparameterhψiinthefilm
3 and thechemical potential when thetemperature T is lower than thefilm critical temperatureTca
c
1
(inthecaseofanuppercriticalsolutiontemperature). Indynamics,welowerT belowthecapillary
0 condensationlinefromaboveTca. Wecalculatethesubsequenttime-developmentassumingnomass
2 c
exchangebetweenthefilmandthereservoir. Intheearlystage,theorderparameterψchangesonly

n in the direction perpendicular to the walls. For sufficiently deep quenching, such one-dimensional
a profiles become unstable with respect to the fluctuations varying in the lateral directions. The
J late-stage coarsening is then accelerated by the hydrodynamic interaction. A pancake domain of

1 thephase disfavored by thewalls finally appears in themiddle of thefilm.

] PACSnumbers: 64.75.St,64.70.qj,68.03.Fg
t
f
o
s I. INTRODUCTION thecapillarycondensationtransitionusingaMonteCarlo
.
t method. Foramicroscopicmodelof2dIsingstripes,Ma-
a ciol ek et al [16] found a (pseudo) CCL using a density-
m Thephasebehavioroffluidsconfinedinnarrowregions
matrix renormalization-group method. For square well
has been studied extensively [1–3]. It strongly depends
d- on the geometry of the walls and on the molecular inter- fluids inslitpores,Singhet al.[17]numericallyexamined
n the crossoverfrom 3d to 2d.
actions between the fluid and the walls. Its understand-
o ing is crucial in the physics of fluids in porous media. It Recently, two of the present authors [18] calculated
c
[ is also needed to study the dynamics of confined fluids. CCL near the bulk critical point using the local func-
tional theory [19, 20], which accounts for the renormal-
In particular, the liquid phase is usually favored by
1 ization effect due to the critical fluctuations. The low-
the walls in one-component fluids, while one component
v ering of the film critical temperature Tca from the bulk
2 ispreferentiallyattractedtothewallsinbinarymixtures c
critical temperature T was shown to be proportional to
31 [w4i–t1h0]t.heInpshuacshe sfaitvuoarteidonbsy, ntahrerowwalrlsegoiornmsamyahyoblde sfiollmede D−1/ν (where ν ∼= 0.6c3) in accord with the scaling the-
0 fraction of the disfavored phase. Hence, in the film ge- ory [11]. Along CCL, our calculations [18] and those by
1. ometry, there appears a first-order phase transition be- Maciol eket al[16]showedstrongenhancementoftheso-
called Casimir amplitudes [21]. Similar first-order tran-
0 tween these states, which forms a line (CCL) ending at
sitions were found between plates [22] and colloids [23]
3 a film critical point outside the bulk coexistence curve
1 in the T-µ∞ plane [1–3, 11–14], where µ∞ is the reser- in binary mixtures containing salt.
v: voirchemical potential[15]. We callit the capillarycon- The aim of this paper is to investigate the phase sep-
i densation transition even for binary mixtures, though aration in near-critical binary mixtures between paral-
X
this name has been used for the gas-liquid phase tran- lel plates using model H and model B [24, 25]. Here,
r sition in porous media [2]. Around CCL, the reservoir phase separation takes place around CCL and the hy-
a
is rich in the component disfavored by the walls for bi- drodynamic interaction is crucial in the late-stage phase
nary mixtures. With increasing the wall separation D, separation. It is worth noting that near-critical fluids in
the film critical point approaches the bulk critical point. porous media exhibit history-dependent frozen domains
Crossoverthenoccursbetweentwo-dimensional(2d)and andactivateddynamicswithnon-exponentialrelaxations
three-dimensional (3d) phase transition behaviors. [26,27]. Togaininsightintosuchcomplicatedeffects,we
For Ising films near the bulk criticality, Fisher and may start with near-critical fluids in the film geometry.
Nakanishi [11] presented the scaling theory of CCL in Treating near-critical fluids, we may construct a univer-
theT-hplane,wherehrepresentsappliedmagneticfield. saltheory with a few materials-independent parameters,
They also calculated CCL in the mean-field φ4 theory. where D much exceeds microscopic spatial scales.
Evans et al. used the density functional theory to cal- In the literature, much attention has been paid to the
culate the inhomogeneous structures in pores [12]. For interplay of wetting and phase separation[28–31], which
a Lennard-Jones fluid in cylindrical pores, Peterson et is referredto as surface-directed phase separation. How-
al.[13] obtained steady gas-liquid two-phase patterns. ever,simulationsincludingthehydrodynamicinteraction
For a lattice gas model, Binder and Landau [14] studied have not been abundant [30, 32–34]. We mention that
′
Tanaka andAraki[33] integratedthe model H equations with the ratio R =ξ /ξ being a universalnumber. We
ξ 0 0
in the semi-infinite situation and Jaiswal et al. [34] per- write ψ in the coexisting two phases as ±ψ with
cx
formedmoleculardynamics simulationto investigatethe
hydrodynamic flow effect between parallelplates. In our ψcx =bcx|τ|β, (2.3)
simulation, the order parameter ψ changes in the direc-
where b is a constant.
tion perpendicular to the walls in the strong adsorption cx
We assume that the bulk free energy F including the
regime. Then, the dynamics is one-dimensional in an
gradient part is of the local functional form [18–20],
early stage but the fluid flow in the lateraldirections ac-
celeratesthelate-stagecoarseningevenundertheno-slip
1
boundary condition on the walls [25, 30, 32–34]. F = dr[f + T C|∇ψ|2]. (2.4)
c
Z 2
On the other hand, Porcheron and Monson [35] nu-
merically studied the dynamics of extrusion and intru- Inthefollowing,wegiveasimpleformforthefreeenergy
sion of liquid mercury between a cylindrical pore and a density f = f(ψ,τ). In our theory, the critical fluctua-
reservoir. Such a process is crucial in experiments of ad- tionswithsizessmallerthanthecorrelationlengthξhave
sorption and desorption between a porous material and already been coarse-grainedat the starting point.
a surrounding fluid [2]. In our simulation we assume no
massexchangeimposingtheperiodicboundarycondition
inthelateraldirections,asintheprevioussimulationsof B. Coexistence-curve exterior
surface-directed phase separation.
The organizationof this paper is as follows. In Sec.II, Outside CX, f is of the Ginzburg-Landau form,
wewillsummarizethe results ofthe localfunctionalthe-
r u
ory of near-critical binary mixtures in the film geome- f/T = ψ2+ ψ4. (2.5)
c
try. We will newly present some results on the phase 2 4
behavior, which will facilitate understanding the phase
Here, we have omitted the free energy contribution for
separationnearCCL.InSec.III,wewillpresentoursim- ψ = 0, whose singular part is proportional to |τ|2−α
ulation results of the phase separation with the velocity
yielding the specific heat singularity. The coefficients r
field (model H) and without it (model B).
and u in f and C in F are renormalized ones in three
dimensions. As in the linear parametric model [36], we
useanonnegativeparameterw representingthedistance
II. THEORETICAL BACKGROUND from the critical point in the τ-ψ plane to obtain
r/τ = C ξ−2wγ−1, (2.6)
Thissectionprovidesthetheoreticalbackgroundofour 1 0
simulation on the basis of our previous paper [18]. The u/u∗ = C2ξ−1w(1−2η)ν, (2.7)
1 0
Boltzmann constant kB will be set equal to unity. C = C w−ην, (2.8)
1
∗
where C and u are constants. We may set C = 1
1 1
A. Ginzburg-Landau free energy by rescaling C1/2ψ → ψ without loss of generality.
1
In the present case, (C ξ )1/2ψ is dimensionless. The
1 0
∗ ∗
We suppose near-critical binary mixtures with an up- constant u is a universal number and we set u =
per critical solution temperature (UCST) T at a con- 2π2/9. ThefractionalpowersofwinEqs.(2.6)-(2.8)arise
c
stantpressure. Theorderparameterψ isproportionalto fromthe renormalizationofthe criticalfluctuations with
c−c , where c is the composition and c is its critical wavenumbers larger than the inverse correlation length
c c
value. The reduced temperature is written as ξ−1. We determine w as a function of τ and ψ by
w =τ +(3u∗C ξ )w1−2βψ2, (2.9)
τ =(T −T )/T . (2.1) 1 0
c c
which is equivalent to wγ = (r + 3uψ2)ξ2/C . Thus,
In our numerical analysis, the usual critical exponents 0 1
w =τ for ψ =0 and τ ≥0, while |ψ|∝wβ for τ =0.
take the following values[25]:
Thederivativeµ=∂f/∂ψatfixedτ denotesthechem-
icalpotential difference betweenthe two components[10,
α=0.110, β =0.325, γ =1.240,
15,19],butitwillbesimplycalledthechemicalpotential.
ν =0.630, η =0.0317, δ =4.815. (2.2) In terms of the ratio S =τ/w, it reads
At the critical composition with τ > 0, the correlation µ 2−α+4(1−α)S+5αS2
length is written as ξ =ξ0τ−ν, where ξ0 is a microscopic Tc = 6[2β+(1−2β)S]ξ02 C1wγψ. (2.10)
length. The coexistence curve in the region τ < 0 is
denoted by CX. The correlationlength on CX is written On CX, we require µ = 0, which yields the equation
as ξ = ξ′|τ|−ν, where ξ′ is another microscopic length 2−α+4(1−α)S+5αS2 = 0 for S. On CX, this gives
0 0
S = −1/σ or w = −στ with σ = 1.714. Together with
Eqs.(2.3) and (2.9), we obtain
b2 =(1+σ)σ2β−1/3u∗C ξ . (2.11)
cx 1 0
W introduce the susceptibility χ=χ(τ,ψ) defined by
T /χ=∂µ/∂ψ=∂2f/∂ψ2. (2.12)
c
For ψ = 0 and τ > 0, we simply obtain χ(τ,0) =
C−1ξ2τ−γ. OnCX,wewriteχ =χ(τ,ψ ). Intermsof
1 0 cx cx
the criticalamplitude ratio R =χ(|τ|,0)/χ for τ <0,
χ cx
the susceptibility on CX reads
χ =R−1C−1ξ2|τ|−γ. (2.13)
cx χ 1 0
Some calculations give R = 8.82 [37]. In terms of χ,
χ
thecorrelationlengthisexpressedasξ =(Cχ)1/2,which
′
yieldsthecriticalamplituderatioR =ξ /ξ =2.99[37].
ξ 0 0
For τ =0, we have ξ ∝|ψ|−ν/β.
C. Coexistence-curve interior
The interior of CX is given by |ψ| < ψ and τ < 0,
cx
where we need to define the free energy density f to ex-
amine two-phase coexistence. We assume a ψ4-theory
with coefficients depending only on τ, where ∂f/∂ψ=µ
and ∂2f/∂ψ2 = T /χ are continuous across the coexis-
c
tence curve. We then obtain
(f −fcx)/Tc =(ψc2x/8χcx)(ψ2/ψc2x−1)2, (2.14) F(lIeGft.) a1n:d ω(Clooc(lozr)Do3n/lTince()rigNhotr)mvaslizz/edD f1odr (pτr/oτfiDle,sµ∞ψ/(µz)D/)ψ=D
where f is the free energy density on CX and χ is (−2,−14.1) (top), (−20,−30.8) (middle), and (−20,−17.8)
cx cx
defined by Eq.(2.13). We also set (bottom). Top: Adsorption-dominatedprofileAwith ψ(z)>
0. Middle: Two profiles B and C on the capillary condensa-
C =C =C |στ|−ην, (2.15) tion line with the same grand potential Ω. In (b’), the two
cx 1
curvesofωloc(z)enclosetworegionswiththesameareaclose
whichisthevalueofC onCX.Therenormalizationeffect to the surface tension σ. Bottom: Three profiles D,E, and F
inside CX is assumed to be unchanged from that on CX with the same τ and µ∞ (see Fig.3). In (c’), Ω increases in
theorder of D, E, and F.
with the same τ. The µ inside CX then reads
µ/T =(ψ2/ψ2 −1)ψ/2χ (2.16)
c cx cx
The surface tension σ between coexisting bulk two is supposed to much exceed D. The fluid is close to the
phases is given by the standard expression, bulk criticality and above the prewetting transition line
[1, 4, 6, 7]. We use our local functional theory, neglect-
σ = 2T (C /χ )1/2ψ2 /3 ing the two-dimensional thermal fluctuations with sizes
c cx cx cx
= A T /ξ2, (2.17) exceeding D in the xy plane.
s c We scale τ and ψ in units of τ ∝ D−1/ν and ψ ∝
D D
where ξ = ξ′|τ|−ν is the correlation length on CX. The D−β/ν, respectively, defined by
0
universal number A is estimated to be 0.075 in our
s
model, while its reliable value is about 0.09 [18]. τD = (ξ0/D)1/ν, (2.18)
ψ = (24β/ν/3u∗C ξ )1/2τβ. (2.19)
D 1 0 D
D. Near-critical fluids between parallel plates
In equilibrium theory, it is convenient to assume that
the fluid between the walls is in contact with a large
We suppose a near-criticalfluid between parallel sym- reservoircontaining the same binary mixture, where the
metric walls in the region 0 < z < D, where D is much order parameter is ψ∞ and the chemical potential is
longer than any microscopic lengths. To avoid the dis-
cussion of the edge effect, the lateral plate dimension L µ∞ =µ(ψ∞,τ). (2.20)
where the space integral dr is within the film, the sur-
face integral dS is onthRe walls atz =0 and D, andh1
isasurfacefieRldsymmetricallygivenonthetwowalls. In
addition, we neglect the surface free energy of the form
dSλ−1ψ2 assumedintheliterature[4–7](orweconsider
tRhe limit λ→∞).
In Eq.(2.21) ω is the local grand potential density
loc
including the gradient part,
1
ω =ω + T C|∇ψ|2, (2.22)
loc s c
2
whereω is the excessgrandpotentialdensity writtenas
s
ωs =f(ψ)−f(ψ∞)−µ∞(ψ−ψ∞). (2.23)
Now minimization of Ω yields the bulk equation,
FIG. 2: (Color online) Phase diagram of a near-critical fluid δF T
in a film for large adsorption in the µ∞/µD- τ/τD plane. δψ =µ− 2cC′|∇ψ|2−TcC∇2ψ =µ∞, (2.24)
The bulk coexistence line is given by τ < 0 and µ∞ = 0.
′
On its left, there appears a first-order capillary condensa- where C = ∂C/∂ψ. The boundary conditions at z = 0
tion line (red bold line) ending at a film critical point, which and D are given by
is calculated from 1d profiles. Displayed also are values of ′ ′
ψ (x,y,0)=−ψ (x,y,D)=−h /C. (2.25)
µ∞/µD in steady two-phase coexistence in our simulation of 1
a2D×2D×D system (×). Startingpoint of oursimulation where ψ′ =∂ψ/∂z.
(t<0) is also shown (◦).
The role of h in this paper is simply to assure
1
the strong adsorption regime ψ /|τ|β ≫ (C ξ )−1/2
0 1 0
[10, 18, 19], where ψ is the boundary value of ψ. This
0
regime is eventually realizedon approachingthe critical-
ity (however small h is). In our simulation, the profile
1
of ψ in the region 0 < z < ξ is nearly one-dimensional
depending only onz evenin two phase states (see Figs.5
and 6). It decays slowly as (ℓ +z)−β/ν for 0 < z < ξ
0
[8–10, 18], where ℓ is a short microscopic length intro-
0
ducedbyRudnickandJasnow[9]. Withthegradientfree
energy in the form of Eq.(2.22), ℓ is expressed as
0
ℓ ∼ξ (C ξ )−ν/2βψ−ν/β ∼D(ψ /ψ )−ν/β, (2.26)
0 0 1 0 0 D 0
in terms of ψ . The excess surface adsorption of ψ in
0
the region 0 < z < ℓ is of order ψ ℓ ∼ ψ1−ν/β and is
0 0 0 0
negligible for large ψ from β/ν ∼ 2, while that in the
0
region ℓ < z < ξ is of order ξ1−β/ν for ξ(τ,ψ ) < D/2
FIG. 3: (Color online) Isothermal curves in the µ∞/µD- 0 m
hψi/ψD plane,which arecalculated for1d profilesatτ/τD = [8–10]. In the strong-adsorption regime we calculate the
−2, -10, and -20. For τ less than its film critical value average of ψ along the z axis,
(= −3.14τD), a van der Waals loop appears. Dotted parts D
of the curves for τ/τD = −10 and -20 are not stable in con- hψi= dzψ/D. (2.27)
tactwithareservoir. PointsA,B,C,D,E,andFcorrespondto Z0
thecurves in Fig.1.
In two-phase states, hψi depends on (x,y).
From Eq.(2.25) it follows the estimation h ∼
1
C(ψ )ψ /ℓ . As h /|τ|βδ−ν →∞, we find
0 0 0 1
Here, µ∞ corresponds to magnetic field h for films of h =ψδ−ν/β(B +B τψ−1/β +···), (2.28)
1 0 1 2 0
Ising spin systems. We are interested in the case µ∞ <
0 (or ψ∞ < 0) and ψ0 > 0, where ψ0 is the value of where B1 and B2 are positive constants. This is the ex-
ψ at the walls. If equilibrium is attained in the total pression for D →∞. In this regime, the surface free en-
system including the reservoir, we should minimize the ergyinEq.(2.21)isgivenby−2B ψ2ν/βA,whereA=L2
1 0
filmgrandpotentialΩ. Includingthesurfacefreeenergy, is the surface area. In our previous paper[18], we ex-
we assume the form, amined the film phase behavior at fixed large ψ , treat-
0
ing the surface free energy as a constant. In our sim-
ulation, we assume the boundary condition (2.25) with
Ω=Z dr ωloc−TcZ dS h1ψ, (2.21) h1/C =1011ψD/D to obtain ψ0/ψD ∼=14.8.
E. Capillary condensation transition
We consider the capillary condensation transition on
thebasisofone-dimensional(1d)profilesψ =ψ(z). From
Eq.(2.25), we have the symmetry ψ(z) = ψ(D−z). In
the region 0<z <D/2, Eq.(2.24) is integrated to give
ψ0 C(ψ)/2 1/2
z = dψ . (2.29)
Z (cid:20)ω (ψ)+Π(cid:21)
ψ s
Here, Π = −A−1∂Ω/∂D is the osmotic pressure. It is
theforcedensityperunitareaexertedbythefluidtothe
plates. In our case, Π < 0, indicating attractive inter-
wall interaction. In the 1d case, it is also written as
Π = f(ψ∞)−f(ψm)−µ∞(ψ∞−ψm)
= −ω (ψ ). (2.30)
s m
At the midpoint z = D/2, we set ψ = ψ(D/2). The
m
fluid at the midpoint can be in the phase favoredby the
walls with ψ ∼ ψ due to the strong adsorption on
m cx
the walls or in the disfavored phase with ψm ∼=ψ∞ <0.
Equation(2.30)indicatesΠ∼=2µ∞ψcx intheformercase
and Π∼=−Tc(ψm−ψ∞)2/2χcx ∼=0 in the latter case, so
Π can be very different in these two cases [38].
Figure 1 displays typical 1d profiles of ψ(z) from
Eq.(2.29) and ω (z) in Eq.(2.22) in the range 0 <
loc FIG. 4: (Color online) Phase diagrams in the τ/τD-ψm/ψD
z < D/2, which will be needed to explain our simu- plane in (a) and in the τ/τD-hψi/ψD plane in (b), where
lation results. Here we set (τ/τD,ψ∞/ψD,µ∞/µD) = ψm is the midpoint value. In (a) and (b), the capillary con-
(−2,−1.14,−14.1) (top), (−20,−1.85,−30.8) (middle), densation curves (bold lines) are calculated from 1d profiles,
and(−20,−1.81,−17.8)(bottom),whereµ∞ismeasured where points for five τ/τD are those along the z axis with
in units of (x,y)=(D,D)and(D,0)inthefinaltwo-phasestatesinour
simulation (see Figs.4 and 5). Bulk coexistence curve is in
µ =T /D3ψ ∝Dν/β−3. (2.31) brokenline. In(b)ourphaseseparationprocessisillustrated
D c D
byarrows, where thetotal order parameter is conserved.
Salient features in Fig.1 are as follows. (i) In Fig.1(a), ψ
is positive in the whole region with < ψ > /ψ = 1.20.
D
(ii) In Figs.1(b) and (b’), the fluid is on the capil-
In Fig.2, we show the capillary condensation line
lary condensation line, where we give two equilibrium
(CCL) from 1d profiles located on the left of the bulk
profiles B and C with the same Ω. Here, B repre-
coexistence line in the τ-µ∞ plane. In our previous pa-
sents an adsorption-dominatedstate with ψ ∼ψ and
m cx per [18], the correspondingphase diagramwas displayed
< ψ > /ψ = 1.99, while for C the film center is oc-
D in the τ-ψ∞ plane. The discontinuities of the physical
cupied by the disfavored phase with ψ ∼ −ψ and
m cx quantities across CCL increase with increasing |τ| van-
<ψ >/ψ =0.187. In the right panel (b’), the integral
D ishing at a film critical point. At this film criticality, τ,
of ω (z) in the region 0 < z < D/2 is the same for B
loc ψ∞, and µ∞ are calculated as
andC.Theenclosedtworegionshavethesamearea24.7
in units of Tc/D2, which is close to the surface tension τ ψ∞ µ∞
σ = 29.2T /D2 at this τ. In addition, Π ∼ −σ/D for B , , =(−3.14,−1.27,−16.3), (2.32)
c (cid:18)τ ψ µ (cid:19)
and Π ∼= 0 for C [38]. (iii) In Figs.1(c) and (c’), the pa- D D D
rametersarethose slightlybelowthe capillarycondensa- where we also have hψi/ψ = 0.989 and ψ /ψ =
D m D
tionline (inFig.2below). Here,therearethreesolutions −0.173. Hereafter, the chemical potential µ∞ on this
with the common τ and µ∞, but < ψ > /ψD is 2.03 for CCL will be written as µca(τ). Our numerically calcu-
cx
(D), 0.028 for (E), and 1.20 for (F) (see Fig.3 below). If lated CCL is well fitted to the linear form,
we perform simulation in contact with a reservoir with
these τ and µ∞, the profile D is realized at long times. µccax(τ)/µD +16.3∼=0.86(τ/τD+3.14). (2.33)
In the very early stage of our simulation, the dynamics
is one-dimensional and the profile F is approached after In Fig.1(b’), the two areas enclosed by the two curves
quenching from A in Fig.1(a) (see Fig.7 below). of ω D3/T are the same (∼ σD2/T ). Thus, for
loc c c
TABLE I: Values of µ∞/µD from our simulation in a finite
2D×2D×Dsystemforfiveτ/τD. Thecorrespondingvalues
of µccax/µD on theCCL from 1d profiles are also shown.
τ/τD −8.5 −10 −13.3 −16.6 −20
µ∞/µD (finite) -25.8 -27.8 -32.2 -36.6 -41.4
µccax/µD (1d) -20.6 -21.8 -24.9 -27.7 -30.8
|τ|/τ ≫1,thesurfacetensionσ andthefreeenergydif-
D
ference per unit area −2µcaψ D are of the same order.
cx cx
SeethesentencesbelowEq.(2.30)andthe explanationof
Fig.1(b’). For |τ|≫τ , it follows the relation,
D
µca ∼−σ/ψ D ∼−|τ/τ |2ν−βµ . (2.34)
cx cx D D
Since 2ν − β ∼= 0.94, the theoretical formula (2.34) is
consistent with the numerical formula (2.33). Note that
Eq.(2.34)is equivalent to the Kelvin equation known for
the gas-liquid transition in pores [1, 2].
We have already presented a special case of three 1d
profiles in Fig.1(c) for (τ/τD,µ∞/µD) = (−20,−17.8).
InFig.3,weshowisothermalcurvesintheµ∞-hψiplane,
which are calculated from 1d profiles with τ/τ = −2,
D
−10, and −20. The relation between µ∞ and hψi is
monotonic for τ/τ ≥ −3.14 (above the film critical
D
temperature), while it exhibits a van der Waals loop for
τ/τ < −3.14 with three 1d states in a window range
D
µ∞1 < µ∞ < µ∞2 [39]. Here, µ∞1 and µ∞2 coincide at
the film criticality. The isothems consist of stable and FIG.5: (Coloronline)Equilibriumtwo-phasestateatτ/τD =
unstable parts characterized by the sign of the film sus- −10 in a 2D×2D×D system: ψ/ψD (top) and ωlocD3/Tc
ceptibility defined by (middle)inthexyplaneatz =D/2(left)andintheyzplane
atx=D(right). Bottom: ψ/ψD (left)andωlocD3/Tc (right)
alongthezaxisfor(x,y)=(D,D)(blueboldline)and(D,0)
χfilm =Tc(∂hψi/∂µ∞)τ. (2.35)
(red bold line), while dotted lines represent 1d profiles from
Eq.(2.29).
In Fig,3, points A,B,C,D,E, and F correspond to the
curves in Fig.1. Dotted parts of the two curves of
τ/τ =−10 and −20 are not stable in the presence of a
D
mass current from a reservoir with common µ∞. III. PHASE SEPARATION DYNAMICS
Previously,some authors [13, 14] calculatedthe stable
partsofisothermsoftheaveragedensityinthefilmversus
We performed simulation of phase separation in a
thechemicalpotential. Inourlocalfunctionaltheory,the
L×L×DcellwithL=2Dimposingtheperiodicbound-
three 1d profiles can be calculated since a unique profile
ary condition along the x and y axes. In this section,
follows for any given set of τ and hψi. In equilibrium
we describe phase separation realized for deep quench-
fluctuation theory of films [1], χ is proportional to
film ing. However, it was not realized for shallow quenching
the variance of the order parameter fluctuations, so its
(|τ|/τ .7), for which χ in Eq.(2.35) is positive.
D film
negativity indicates thermodynamic instability.
Furthermore, Fig.4 gives the phase diagrams in the τ-
ψ and τ-hψi planes. Bold lines represent the capillary
m
condensationcurvefrom1dprofilesasinFig.2. Insteady A. Dynamic equations and simulation method
two-phase states in our simulation, ψ and hψi depend
m
on x and y, so points for five τ represent ψ (x,y) = Supposinganincompressiblefluidbinarymixturewith
m
ψ(x,y,D/2) and hψi(x,y) with (x,y) = (D,D) and a homogeneous temperature, we use the model H equa-
(D,0)(seeFigs.5and6). Theformerlinepassesthrough tions [24, 25, 40]. The order parameter ψ is a conserved
a domain of the phase disfavored by the walls and the variable governed by
latter through the favored phase only. Phase diagrams
similar to Fig.4(b) have been obtained in experiments of ∂ψ δF
=−∇·(ψv)+λ∇2 , (3.1)
the capillary condensation in porous media [2]. ∂t δψ
where χ is the susceptibility on CX in Eq.(2.13) and
cx
D is the mutual diffusion constant of the Stokes form,
ξ
D =T /6πη¯ξ, (3.4)
ξ c
with ξ = ξ′|τ|−ν being the correlation length on CX. In
0
oursimulation, |ψ| isoforderψ atz =D/2(see Figs.5
cx
and 6 below), which supports Eqs.(3.3) and (3.4). The
viscosity η¯ exhibits a very weak critical singularity and
may be treated as a constant independent of τ.
In this paper, we also performed simulation for model
B without the hydrodynamic interaction [24], where ψ
obeys the diffusive equation,
∂ψ δF
=λ∇2 . (3.5)
∂t δψ
The kinetic coefficient λ is assumed to be given by
Eqs.(3.3) and (3.4) as in the model H case. Then, com-
paring the results from the two models, we can examine
the role of the hydrodynamic interaction in phase sep-
aration. Model B has been used to investigate surface-
directed phase separation in binary alloys [31].
InintegratingEqs.(3.1)and(3.5),themeshlengthwas
∆x = D/32 and the time interval width was ∆t = 2×
10−6t . The initial state was the 1dprofile A in Fig.1(a)
0
at τ/τ = −2 with small random numbers (∼ 10−4)
D
superimposedatthemeshpoints. Att=0,wedecreased
FIG.6: (Coloronline)Equilibriumtwo-phasestateatτ/τD = τ toafinalreducedtemperature. Fort>0,therewasno
−20 in a 2D×2D×D system: ψ/ψD (top) and ωlocD3/Tc massexchangebetweenthefilmandthereservoirsothat
(middle) in the xy plane at z = D/2 (left) and in the yz thetotalorderparameter drψwasfixedat1.20ψ DL2.
plane at x=D (right). Bottom: ψ/ψD (left) and ωlocD3/Tc We will measure time afteRr quenching in units ofD
(right) along the z axis for (x,y) = (D,D) (blue bold line)
and (D,0) (red bold line), which are closer to the 1d profiles
t =D2/D , (3.6)
(dotted lines) from Eq.(2.29) than in Fig.5. 0 ξ
which is the mutual diffusion time in the film assumed
to be much longer than the thermal diffusion time. We
note that the natural time unit in bulk phase separation
where λ is the kinetic coefficient and the functional
hasbeentheorderparameterrelaxationtimet =ξ2/D
derivative δF/δψ may be calculated from Eq.(2.4) with ξ ξ
[25, 40, 41]. Here, t /t = (D/ξ)2 = R2|τ/τ |2ν ∼ 100
the aid of Eqs.(2.10), (2.15), and (2.16) outside and in- 0 ξ ξ D
side CX. We neglect the random source term originally for |τ/τD|∼10.
present in critical dynamics [24, 25], because we treat
thedeviationsmuchlargerthanthethermalfluctuations.
The velocity field v satisfies ∇·v = 0 and vanishes at B. Steady two-phase states
z = 0 and D. In the Stokes approximation[40], v is de-
termined by For sufficiently deepquenching, we realizedphase sep-
aration to find a steady two-phase state at long times
η¯∇2v =∇p0+ψ∇(δF/δψ), (3.2) both for model H and model B. We could also calcu-
late this final state more accurately from the following
whereη¯istheshearviscosityandtheroleofapressurep0 relaxation-type equation,
istoensure∇·v =0. SeeAppendixfortheexpressionof
the stresstensorinnear-criticalfluids andthe derivation ∂ψ δF δF
=−L − , (3.7)
of Eq.(3.2). ∂t 0(cid:20)δψ (cid:28)δψ(cid:29) (cid:21)
t
The kinetic coefficients λ and η¯ should be treated as
renormalized ones [24, 25, 40] (see the last sentence of whereL isaconstantandhδF/δψi isthespaceaverage
0 t
Subsec.IIA).Inthevicinityofthebulkcoexistencecurve, ofδFδψattimet. Becauseofitssimplicity,weintegrated
λ may be approximated by Eq.(3.7)withafinemeshlengthof∆x=D/64. Thedata
points in Figs.2 and 4 and the snapshots in Figs.5 and 6
λ=χ D /T , (3.3) are those from the steady states of Eq.(3.7).
cx ξ c
the deviation of ψ /ψ or hψi/ψ by of A , we obtain
m D D χ
that of µ∞/µD.
In Figs.5 and 6, we display the final profiles of ψ and
ω in the xy plane at z = D/2 (left) and in the zy
loc
plane at x = D (right) for τ/τ = −10 and −20. In
D
these cases, ξ onCX is 0.078D and0.051D,respectively,
whichisoftheorderoftheinterfacethickness. Displayed
in the bottom panels are 1d profiles of ψ and ω along
loc
the z axis for the two lateral points (x,y) = (D,D) and
(D,0). These profiles are rather close to the 1d profiles
from Eq.(2.29) in accord with Fig.4.
C. Time evolution
Both for model H and model B, early-stage time-
evolution proceeds as follows. Just after quenching, ψ
changes only along the z axis to approach the 1d profile
atthe final τ with fixedhψi (see Fig.3). If this 1dprofile
satisfies the instability condition χ < 0, it follows 3d
film
spinodaldecomposition. Ontheotherhand,ifitsatisfies
FIG.7: (Coloronline)Early-stagetime-evolutionofψ/ψD for
model H (left) and model B (right) after quenchingfrom the the stability conditionχfilm >0, it remainsstationaryin
profile A in Fig.1(a) to τ/τD = −20. Shown are the profiles simulation without thermal noise.
along the z axis for (x,y) = (D,D) at t/t0 = 0, 0.1, 0.2, Figure 7 displays ψ after quenching to τ/τD = −20.
and 0.3 (top) and along the x axis for (y,z) = (1.5D,0.5D) In the top panels, it is plotted along the z axis with
(bottom) at later times. In the top panels, the profile F in (x,y)=(D,D)att/t =0,0.1,0.2,and0.3. Thevelocity
0
Fig.1(c)isapproachedwithoutnoticeabledifferencesbetween field nearly vanishes for model H, so there is almost no
thetwo models. Inthebottom panels, fluctuationsin thexy
difference between the results of these models. However,
plane appear and coarsening is much quickened for model H
in the bottom panels, the 1d profile becomes unstable
than formodel B, wheret/t0 = 0.5, 0.6, and0.7 for model H
with respect to the fluctuations varying in the xy plane
and 1.0, 1.1, and 1.2 for model B.
for t/t & 0.5. The velocity field grows gradually for
0
modelH.Inthissecondtimerange,thedomainformation
is much quicker for model H than for model B.
In the left panel of Fig.8, we show the free energy de-
crease ∆F =F(t)−F(+0) at τ/τ =−20 as a function
D
oftformodelHandmodelB.Here,F(t)isthetotalbulk
free energy in Eq.(2.4) with F(+0) being its value just
afterquenching. Itsdecreaseisacceleratedwithdevelop-
mentof the fluctuations in the xy plane. The coarsening
isslowerformodelBthanformodelHbyabout5times.
In the right panel of Fig.8, we show the characteristic
velocity amplitude v¯(t) for model H, which is defined by
FIG. 8: (Color online) Time-evolution of the normalized free v¯(t)2 = dr|v|2/DL2. (3.8)
energy decrease ∆F(t)/Tc for model H and model B (left) Z
and v¯(t)t0/D for model H (right) after quenching from the
For τ/τ =−20, v¯(t) is equal to 0.021, 0.106, and 0.485
profile A in Fig.1(a). At points (a), (b), (c), and (d) (right) D
snapshots of ψ and v will be given in Fig.9. at t/t0 =0.5, 0.6, and 0.7, respectively, increasing up to
0.8, in units of D/t .
0
In Fig.9, we show late stage snapshots of ψ and v
in the xz plane (top) and in the xy plane (middle) for
In Fig.4, the deviations of ψ and hψi from those on τ/τ =−20 for model H, where t/t is equal to (a) 0.7,
m D 0
CCLaresurprisinglysmall,thoughthelateraldimension (b) 0.9, (c) 1.1, and (d) 1.3 after quenching. In (a) we
L is only 2D. However, in Table 1, the final two-phase can see a network-like domain of the disfavored phase.
values of µ∞ considerably deviate from those on CCL. In (b) three domains can be seen, where the middle one
This may be ascribed to the relatively small size of the is being absorbed into the bottom one, soon resulting in
susceptibilityχ forthesecases. Thatis,ifwesetχ = two domains at t/t ∼ 1.0. This process gives rise to a
cx cx 0
A−1T ψ /µ on CX, we obtain A = 120 and 283 for dipinv¯(t)inFig.8(b)sincethesetwodomainsareconsid-
χ c D D χ
τ/τ = −10 and −20, respectively. Here, if we multiply erably apart. In (c) and (d), furthermore, coalescence of
D
FIG. 9: (Color online) Cross-sectional velocity field v (arrows) and order parameter ψ (in gradation according the color bar)
in phase separation for τ/τD = −20 in model H, where t/t0 is equal to (a) 0.7, (b) 0.9, (c) 1.1, and (d) 1.3 after quenching.
Displayedare(vx,vz)andψ inthexz planeaty=D (top),thoseaty=1.5D (middle),and(vx,vy)andψ inthexy planeat
z =D/2 (bottom). Arrows below panels indicate typical magnitudes of the velocities, where σ is the surface tension and η¯is
theshear viscosity. Final state is given in Fig.6.
thesetwodomainsistakingplace. The arrowsbelowthe plates in the strong adsorption regime around the capil-
panels indicate the typical velocity, 0.15σ/η¯ or 0.3σ/η¯, larycondensationline (CCL).Using model Handmodel
where σ is the surface tension and η¯ is the shear viscos- B,simulationhasbeenperformedina 2D×2D×D cell.
ity. Note that the typical velocity in the late-stage bulk We summarize our main results.
spinodal decomposition is given by v =0.1σ/η¯[25, 41], (i) In Sec.II, we have presented the singular free energy
c
which follows from the stress balance σ/R ∼ 6πη¯v /R with the gradient part outside and inside the bulk coex-
c
with R ∼ v t being the typical domain length. In our istence curve. Applying it to near-criticalfluids between
c
case, from Eq.(2.17) these velocities are related as parallel plates, typical 1d profiles have been given in
Fig.1. CCL has been plotted in the τ-µ∞ plane in Fig.2.
D/t =T /6πη¯ξD =(ξ/6πA D)(σ/η¯). (3.9) Thepointsforthe steadytwo-phasestatesfromoursim-
0 c σ
ulation are located in the left side of CCL. In Fig.3, we
Thus, D/t =0.036σ/η¯for τ/τ =−20 in Fig.9. have also found the van der Waals loop of isothermal
0 D
curves in the hψi-µ∞ plane, where τ is smaller than the
film critical value. The phase diagrams have been plot-
ted in the τ-ψ and τ-hψi planes in 3d (bulk) and 2d
IV. SUMMARY AND REMARKS m
(film) in Fig.4. The Kelvin relation (2.34) has also been
obtained, since the osmotic pressure Π is of order −σ/D
Insummary,wehaveexaminedthephaseseparationin
right below CCL [38].
anear-criticalbinarymixturebetweensymmetricparallel
(ii) In Sec.III, we have first displayed the cross-sectional 5) The static part of this work is applicable to any
profilesofψ andω insteadytwo-phasestatesinFigs.5 Ising-like near-critical systems and can readily be
loc
and 6. The profiles along z axis for (x,y) = (D,D) generalized to n-component spin systems. In the
and (D,0) closely resemble the corresponding 1d pro- dynamics, we have used model H with a homogeneous
files. For quenching to τ/τ = −20, we have examined temperature and incompressible flows. On the other
D
time-evolutionofψ. Itoccursonlyalongthez axisinthe hand, in one-component near-critical fluids, the latent
veryearlystageinthetoppanelsofFig.7,wherethereis heat released or absorbed at the interfaces gives rise to
no difference between the results of model H and model significant hydrodynamic flow because of the enhanced
B.Subsequently,inhomogeneitiesappearinthexy plane. isobaric thermal expansion [25, 42]. Also promising in
The free energydecrease∆F(t)=F(t)−F(+0)andthe future should be extension of this work to near-critical
typical velocity amplitude v¯(t) defined in Eq.(3.8) have fluids in porous media.
been plotted in Fig.8. The velocity field considerably
quicken the interface formation and the coarsening for
modelHthanformodelB.Profilesofψ andv inthelate
stage coarsening due to the flow have been presented in
Acknowledgments
Fig.9,wherethe domaincoalescencecanbe seenandthe
maximum velocity is of order 0.1σ/η¯as in bulk spinodal
This work was supported by Grant-in-Aid for Scien-
decomposition [41].
tific Research from the Ministry of Education, Culture,
We make some remarks.
Sports,ScienceandTechnologyofJapan. S.Y.wassup-
1) In the static part of our theory, we neglect the
ported by the Japan Society for Promotion of Science.
thermal fluctuations varying in the lateral directions
A.O.wouldliketothankSanjayPuriforinformativecor-
withwavelengthslongerthanD. Thus this 2dtransition
respondence.
exhibits mean-field behavior. In fact, the curves of ψ
m
vs τ and hψi vs τ are parabolic near the film criticality Appendix: Stress tensor in near-critical fluids
in Fig.4.
2) In our simulation, we soon have only two or three In near-critical fluids, we treat slow flows with low
domains in the cell as in the bottom panels of Fig.9. Reynolds numbers. The total (reversible) stress tensor
The lateral dimension L=2D in this paper is too short is givenby p δ +Π , wherep is nearlyhomogeneous
0 ij ψij 0
to investigate the domain growth law in the xy plane. throughout the film and the reservoir. The Π is the
ψij
Simulation with larger L/D should be performed in stress tensor due to the composition deviation[25],
future.
3) Fromthe van der Waals loopof the isothermalcurves
Π =p δ +C(∇ ψ)(∇ ψ), (A1)
ψij ψ ij i j
in the hψi-µ∞ plane in Fig.3, we may predict how phase
separation proceeds after quenching. We have examined
where ∇ =∂/∂x . The diagonal part p is written as
i i ψ
phase separation via spinodal decomposition. However,
in real experiments, phase separation may occur via nu- 1
cleation for metastable 1d profiles. Note that hysteretic pψ = ψµ−f − (C+C′ψ)|∇ψ|2−Cψ∇2ψ
2
behavior has been observedin phase-separating fluids in
= ψ(δF/δψ)−f −C|∇ψ|2/2, (A2)
pores and has not been well explained [1, 2, 26, 27].
4) In a number of experiments and simulations of ′
whereC =∂C/∂φ. ThesecondpartinEq.(A1)contains
surface-directed phase separation [29, 33, 34], compo-
off-diagonal components relevant for curved interfaces.
sition waves along the z axis have been observed near
In deriving Eq.(3.2), we use the relation,
the wall in the early stage. In these cases, the degree of
adsorption has changed appreciably upon quenching. In
the strong adsorption regime in this paper, 1d dynamics ∇jΠψij =ψ∇i(δF/δψ). (A3)
occurs in the initial stage, but there are no composition Xj
waves as in the top panels of Fig.7.
[1] R. Evans and U. M. B. Marconi, J. Chem. Phys. 86, [4] J. W. Cahn, J. Chem. Phys. 66 3667 (1977).
7138 (1987); R. Evans, J. Phys.: Condens. Matter 2, [5] K.Binder,inPhaseTransitionsandCriticalPhenomena,
8989 (1990). C. Domb and J. L. Lebowitz, eds. (Academic, London,
[2] L.D. Gelb, K.E. Gubbins, R. Radhakrishnan, and M. 1983), Vol. 8, p. 1.
Sliwinska-Bartkowiak,Rep.Prog.Phys.62,1573(1999). [6] P.G. deGennes, Rev. Mod. Phys. 57, 827 (1985).
[3] K. Binder, D. Landau, and M. Mu¨ller, J. Stat. Phys. [7] D.BonnandD.Ross,Rep.Prog.Phys.64,1085(2001).
110, 1411 (2003); M. Mu¨ller and K. Binder, J. Phys.: [8] B. M. Law, Prog. Surf. Sci. 66, 159 (2001).
Condens.Matter17, S333(2005); K.Binder,J.Horbach, [9] J. Rudnick and D. Jasnow, Phys. Rev. Lett. 48, 1059
R.Vink,andA.DeVirgiliis,SoftMatter,4,1555(2008). (1982); ibid. 49, 1595 (1982)