Phase separation in binary hard-core mixtures
Marjolein Dijkstra and Daan Frenkel
FOM Institute for Atomic and Molecular Physics, Kruislaan 407, 1098 SJ Amsterdam, The Netherlands
Jean-Pierre Hansen
Laboratoire de Physique, Uniti de recherchg 1325 Associke au CNRS, Ecole Normale Sup&ieiwe de Lyon,
46, All&e d’ltalie, 69364 Lyon Cidex 07, France
(Received 14 March 1994; accepted 1 April 1994)
We report the observation of a purely entropic demixing transition in a three-dimensional binary ’
hard-core mixture by computer simulations. This transition is observed in a lattice model of a binary
hard-core mixture of parallel cubes provided that the size asymmetry of the large and small particles.
is sufficiently large (33, in the present case). In addition, we have performed simulations of a single
athermal polymer in a hard-core solvent. As we increase the chemical potential of the solvent, we
observe a purely entropy-driven collapse of the polymer: the scaling of the radius of gyration R, of
the polymer with the number of segments N changes from that of a polymer in a good solvent to that
of a collapsed polymer. Both for the study of the hard-core demixing and of the polymer collapse,
it was essential to use novel collective Monte Carlo moves to speed up equilibration. We show that ;
in the limit ~t/oz-+O, the pair distribution function for an off-lattice binary hard-core mixture of
parallel cubes with side lengths ol and cr2 diverges at contact for the large particles. For the lattice
system, we
calculated the pair distribution functions g(r) up to the fourth virial coefficient. The
difference in g(r) at contact for a binary system and a pure system at the same packing fraction
gives a rough criterion, whether the mixture phase separates,
1. INTRODUCTION
One of the most striking phenomena in mixtures is phase
separation. Until recently, it was generally believed that
phase separation in simple mixtures is induced by energetic
effects. An important question is therefore whether a phase
separation can take place in an “athermal” mixture, i.e., a
mixture for which the energy of mixing depends linearly on
composition. The simplest model for an athermal mixture is
a binary system containing large and small hard spheres. In
1964, Lebowitz and Rowlinson3 showed that, within the
Percus-Yevick closure of the Ornstein-Zernike equation,
fluid hard sphere mixtures are stable with respect to phase
separation. Thus far, -computer simulations4-s found no evi-
dence for fluid-fluid phase separation in additive hard sphere
mixtures. In this
context, additive
means
~~~ = (o-& -t- (7&/2, where aij denotes the distance of clos-
est approach of particles of type i and j. In contrast, it is
usually possible to have phase separation in a mixture with
positive nonadditivity, i.e., a,,> (0;1~ f ~aB)/2,a-i’ be-
cause for such systems, the pure phases can till space more
effectively than the mixture. However, recently Biben and
Hansen” showed’that within the ‘Rogers-Young closure of
the Ornstein-Zernike equation, additive hard sphere mix-
tures will be unstable at high densities for diameter ratios
larger than 5. For hard spheres, this thermodynamically self-
consistent closure is known to be more accurate than the
Percus-Yevick approximation. Similar predictions have sub-
sequently been made using other approximations.‘3 In prin-
ciple, computer simulation is a suitable tool to investigate
whether phase separation can take place in purely athermal
mixtures. Unfortunately, direct simulation of phase separa-
tion in mixtures of very dissimilar spheres is difficult be-
cause
of
slow equilibration.14
The numerical’ difficulties are
less severe for lattice models of hard-core.mixtures. We have
therefore performed simulations of lattice models of additive
hard-core particles.
The results of grand canonical Monte Carlo simulations
of a mixture of large and small cubes’ with side ratios of 2 or
3 on a cubic lattice are reported in this paper. In Sec. II, the
simulation method is discussed and the results are presented.
The more technical aspects of the simulation method are ex-
plained in the Appendices. A question related to phase seps-
ration induced by purely entropic effects is whether an ather-
ma1 polymer in a hard-core solvent undergoes a collapse
transition. The results of a simulation of an athermal polymer
in a hard-core solvent are reported in Sec. III. In Sec. IV, we
show that in an off-lattice system of large and small parallel
cubes, the pair distribution function of the large cubes di-
verges at contact in the limit oi/oi+O, where cr, and cr2 are
the side lengths of the cubes. This tendency of the large
cubes to stick.together is strongly suggestive of phase sepa-
ration. For the lattice
system,
we computed the pair distribu-
tion functions of the binary hard-cube mixture and compared
them with the pair distribution functions for a pure system at
the same packing fraction.
II. COMPUTER SIMULATIONS OF A BINARY MIXTURE
OF HARD PARALLEL CUBES
A. Simulation method
The model that we consider is a mixture of large and
small parallel hard cubes on a lattice. The diameter of each
cube corresponds to an even number of lattice spacings. This
model is clearly additive: at close packing we can fill space
just as well in the mixed phase as in the pure phases. Hence
there is no trivial volume-driven demixing. In our simula-
tions, we considered mixtures of cubes with side ratios of 2
J. Chem. Phys. 101 (4), 15 August 1994
0021-9606/94/101(4)/3179/11/$6.00 8 1994 American Institute of Physics
3179
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
3180 Dijkstra, Frenkel, and Hansen: Phase separation in hard-core mixtures
or 3 on a lattice with, respectively, 28X28X28 and 30
X30X30 lattice points. In both cases, the side of the small
cubes is equal to 2. We performed grand-canonical Monte
Carlo (GCMC) simulations, where the independent variables
are the fugacities of the large and the small cubes zI and Z, .
Three types of trial moves can be performed:
(i) random displacement of a particle in the box;
(ii) changing the identity of a particle (large to small or
vice versa);
(iii) removal or insertion of a particle.
However, the acceptance ratio fdr a random displace-
ment of a large particle is small in a dense system of small
particles, as the displacement of .a large particle is ‘strongly
hindered by the small ‘particles. In order to speed up equili-
bration, we used collective particle moves that employed a
generalization of ‘the configurational-bias Monte Carlo
scheme of Ref. 15. In this -approach, the large particle was
moved to a random trial position. Typically several small
particles would
occupy
this region in space. These particles
were then moved to the volume vacated by the laige particle
and inserted using RosenblGth sampling.14 The trial move
was then accepted with a probability deterniined by the ratio
of the Rosenbluth weights of the new and old configurations.
Of course, a trial move would be rejected immediately if it
resuIted in overlap of two or more large particles. In order to
investigate the ipfluen.ce of the surface-to-volume ratios on
phase separation, we also, performed simulations of a three-
dimensional system of hard .parallel platelets (size 6X6)<2)
and cubes (size 2y2X2) on a cubic lattice of 30X30X30.
For the sake of. comparison, we also. simulated a two-
dimensional mixture gf hard squares (size 6X6 and 2X2)
and a two-dimensional system of hard squarqs. (size 2X2)
and parallel hard rods (18X2).
b. Results
~I . . _-*,_
GCMC simulations were performed on a mixture of
large and small cubes for a range of different values for the
fugacities z1 and zS. In each simulation, we computed the
volume fractions of the large and small cubes. To allow
faster equilibration, subsequent runs were started from pre-
viously equilibrated configurations at a higher or lower
fugacity. Most runs’donsisted of 1X105-1X106 cycles per
particle. In each cycle, we attempt a random displacement of
a particle in the box, and we try to change the identity of a
random particle. Once every ten cycles, -a removal or inser-
tion of a particle 5n the box is attempted. In Figs. 1 and 2, the
fugacity of the large cubes is plotted vs the volume fraction
of the large cubes and vs the volume fraction of the small
cubes for a mixture of side ratio 3. The different curves in the
figures are computed at constant fugacity of the small cubes.
In Fig. 3, we plot the fugacity of the large cubes versus the
volume fraction of the large cubes for a mixture of cubes
with size 4X4X4 and 2X2X2. For a side ratio of 3, we
observed that upon increasing the fugacity of the solvent
(i.e., of the small particles), the slope of the curves of con-
stant solvent fugacity tends to zero at the inflection point. For
still higher solvent fugacity, we find two. different volume
fractions of the large cubes for the
same
fugacity of the large
particles (see Fig. l)..In Fig. 2, .we observe the same flatten-
I
I I I
I
(
0 0.2 0.4 0.6 0.8 1
G
Large
FIG. 1. The fugacity of the large cubes (size 6X6X6))vs the volume fraction
of the large cubes at different fugacities (z,=O, 100, 500, 1X103, 1.5X 103,
2X103, and 5X103) of the small particles (size 2X2X2). The star denotes
the critical pdint (zl=6.02X-1C17; zb= 1.63X 103, and &se=0.57) and is
derived by a global fit of all isofugacity curves. The rest of the binodal
(dashed) curve could not be estimated as accurately as the critical point.
This curve should therefore reconsidered as a guide to the eye.
ing of the isofugacity curves. For still higher solvent fugac-
ity, we again find DVO different volume fractions of the small
cubes for the same fugacity for the large cubes. This is ex-
actly what we expect for a demixing transition. The critical
point can now be derived from a global fit of all isofugacity
curves. When we take this into account, we can sketch the
demixing region, as shown in Fig. 1. As the cubes are
equally likely to be found on all lattice sites in both phases
(i.e., there is no sublattice ordering), the present phase tran-
sition corresponds most closely to liquid-liquid coexistence
in an off-lattice system. The simulations in the demixing re-
gion are very time consuming and therefore we have not
attempted to locate the binodal curve more accurately. For a
side ratio of 2, we find no evidence for a similar flattening of
the isofugacity curves and we did not fiild two different vol-
ume fractions of the large cubes for the same fugacity of the
large particles. Thus we find no evidence for a demixing
transition in the less asymmetric system. However, if-instead
of cubes, .we consider platelets, with a comparable volume
(6X6X2, instead of 4X4X4), we again observe demixing
(see Figs. 4 and 5). This is plausible as the clustering of the
particles is controlled by depletion forces. When two large
particles are brought into contact, the volume accessible to
the small particles increases by an amount that is propor-
tional to the diameter of the small particles and the area of
contact of the large particles. The resulting gain in entropy of
J. Chem. Phys., Vol. 101, No. 4, 15 August 1994
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
Dijkstra, Frenkel, and Hansen: Phase separation in hard-core mixtures
3181
0
0.1 0.2 0.3 0.4
$
Small
PIG. 2. The fugacity of the large cubes (size 6X6x6) vs the volume fraction
of the small cubes
at diierent fugacities
(z,= 0, 100,
500, 1
X 103,
1.5X 103,
2X 103, and 5X IO’) of the small particles (size 2X2X2).
1
lOI
lo2
103
lo*
lo6
%I
2 lo6
N
lo7
1O8
IO9
lOi
1or1
1ol2
0
0.2 0.4 0.6 0.8
1
+
Large
PIG. 3. The fngacity of the large cubes (size 4X4X4) vs the volume fraction
of the large cubes at different fugacities (z,= 100, 1X 103, 1X104, and
1x10s) of the small particles (size 2X2X2).
lo1
lo2
lo3
lo4
lo5
h:
2 lo6
z
N
lo7
lo8
LO9
1o1O
loll
I-II.-. J.-.--l
0
0.2 0.4. 0.6
0.8 1
@
Plates
PIG. 4. The fugacity of the platelets (size 6X6X2) vs their volume fraction
at different fugacities (z,=SX lo’, 103, 1O4a 5X104, 10’. 2X105, and
5X
10”) of the small particles (cubes of size 2~2x2).
the solvent is the driving force that makes the large particles
cluster. The larger the surface-to-volume ratio ‘of the large
particles, the stronger is the tendency to demix.
In a two-dimensional ‘mixture of squares with size 6X6
and 2X2, we did not tind demixing. At thesame “volume
ratio,” but with rods of size 18X2 instead of the large
squares, the system appeared to approach a spinodal, but we
were not able to reach it.
111. SIMULATIONS OF AN ATHERMAL POLYMER
The driving force that makes the particles cluster de-
pends on the surface-to-volume ratio of the large particles.
This argument should apply not only to rigid particles, such
as rods and disks, but also to flexible particles, such as linear
polymers. We therefore also looked for entropic demixing in
an athermal polymer solution. In fact, in this case, we did not
study the demixing directly. Rather, we looked for a closely
related phenomenon, namely, the solvent-induced collapse of
an isolated polymer. This collapse signals the transition from
the good-solvent to the poor-solvent regime. There are com-
pelling theoretical arguments to assume that a polymer col-
lapse must occur in an athermal polymer solution when the
polymer-solvent interaction is nonadditive.“’ In order to in’-
vestigate if such a collapse can occur in an “additive” ather-
ma1 polymer solution, we performed simulations of a single
hard-core polymer in a solvent of cubes, where the size of
the cubic monomers. of the polymer is the same as the size of
the solvent molecules, namely, 2X2X2 in units of the lattice
J. Chem. Phys., Vol. 101, No. 4, 15 August 1994
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
3182
Dijkstra, Frehkel, and Hansen: Phase separation in hard-core mixtureS
lo1 -1---r-
lo2
lo3
10”
lo5
3
2 lo6
t?
lo7
lo8
10Q
1o1O
101
0
0.2
0.4. 0.6
G
Cubes
:.
-J
FIG. 5. The fugacity
of
the platelets
(size 6X6X2) vs volume fraction of the
small cubes at different fugacities (z;=SXlO’, 103, 104, 5X104, lo?,
2X
10’. and 5X 10’) of the small particles (cubes of size 2X2X2).
spacing. For the simulation of the hard-cork poljriner in a
solvent, we used the con@ura&onal bias Monte Carlo
(CBMC) method15 for generating polymer conformations
and the GkMC method for the solvent. Thus, we performed
1000 [----
t
. . .
.
‘ioo :
NEl
lx
10
i
--. . . . ..I ‘.I”’
5
0
,’ I- bb ..,...1 .,.A I
1 10
100 1000
N
._
FIG. 6. The mean square of the radius of gyration (Ri) of a polymer with
cubic monomers (size 2X2X2) with and without solvent (cubes of size:
2X2X2) vs the number of segments. Average solvent volume fraction 0.0
(circles), 0.3 (crosses), and 0.7 (open diamonds). Note that
R: scales
as d”
with the number of monomers (N). For the two low density runs, we End
~=.0.56t0.02 (for a polymer in a good solvent ~-0.58). The high density
run yields 1~0.34rtO.02. Fdr a collapsed polymer, we expect “Euclidian”
scaling Y= l/3.
FIG. 7. A snapshot of conformation of a hard-core lattice poiymer (N
= 100) without hard-core solvent (top) and in a solution of hard-core mono-
mers with a volume fraction of 0.7 (bottom). Note the solvent-induced col-
lapse.
the simulations on a system with fixed volume and tempera-
ture and at constant chemical potential of the solvent. The
following trial moves have been performed:
(i) random displacement of a solvent molecule in the
box;
(ii) removal or insertion of a solvent molecule;
(iii) regrowing a fraction of the polymer at either end of
th5 chain;
(iv) regrowing a fraction of the polymer that does not
include a free chain-end:
The last trial move was essential for a faster equilibra-
tion of the polymer and is explained in more detail in Ap-
pendices B and C. The dependence of the radius of gyration
R, on the number of monomers N is shown in Fig. 6 at
different chemical potentials of the solvent. Figure 6 shows
that when the chemical potential, and hence the volume frac-
J. Chem. Phys., Vol. 101, No. 4,15 August 1994
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
Dijkstra, Frenkel, and Hansen: Phase separation in hard-core mixtures
3183
tion of the solvent, is increased, the polymer undergoes a
collapse transition. We find that the square of the radius of
gyration
Ri
scales as N2’.
For the two low chemical poten-
tials that correspond to solvent volume fractions of 0 and
about 0.3, we find v=O.56?0.02, which corresponds to the
case of a polymer in a good solvent (~0.58). For high
chemical potential (average solvent volume fraction of 0.7),
we find v=0.34?0.02. For a collapsed polymer we expect
“Euclidian” scaling v= l/3. Figure 7 shows the very drastic
change in the polymer shape when we go to high chemical
potential of the solvent.
IV. VIRIAL EXPANSIONS FOR A BINARY MIXTURE OF
HARD PARALLEL CUBES
A. Continuous system-behavior of g**(r) for y+O
In order to gain a better understanding of the physical
origin of phase separation in a hard-cube mixture, we ana-
lyzed the density expansions of the pair-distribution function
of a closely related model, viz., parallel hard cubes off lat-
tice. For the one-component parallel hard-cube model, the
virial expansion of the pressure and the pair correlation func-
tion has been reported in the literaturei6-l8 up to terms of
order p’. Here we consider the virial expansion of the same
quantities for the binary mixture. For two cubes with side
length cl and era at rl and r.,, the Mayer function is defined
bY
f12hvr2)=-II
@(c+12-lrl,i-r2,il),
i=l
(1)
where a(x) is the Heavyside step function and
cri2=(oi+a2)/2. We computed the virial coefficients, up to
the fourth, for the pair distribution function of the large par-
ticles, i.e., g2a(x,y,z), as the latter quantity is expected to
show evidence of incipient clustering that precedes demix-
ing. If we define +=&x,y,z) as the potential energy of a
pair of particles, one at the origin and the other at (x,y ,z),
the virial expansion for g22(x,y,z)exp(P+) in the densities
p1 and p2 ist9
The diagrams consist of circles and bonds. The circles
represent the particle coordinates and the sizes of the circles
denote the sizes of the particles. Open circles correspond to
coordinates that are not integrated over; the black circles
represent the variables of integration. The bonds between the
circles are associated with the Mayer function. For example,
the first diagram in the pair distribution function denotes
I
dr3 f13(rl,r31f23(r22r3),
(3)
where particles 1 and 2 are large cubes and particle 3 is a
small cube. Note that the three-dimensional integral factor-
izes into a product of three one-dimensional integrals. We
calculated s22i~2,~2.~&q@#4 and g22~0,0,~2kq@4~
up to the fourth virial coefficient. The expressions for these
pair distribution functions in terms of the packing fractions
q 'ofpI
and R=
dp2 are
=1+771+772+777
3 3
-~‘~,+zT;z
+A
A2
11 19
+v1v2 q-2h+5+r +z
+9&
(5)
where X is~ equal to a1/02. For cubes facing each other
(O,O,gJ, the divergence in g22(x,y,z) for A+0 already ap-
pears in the third virial coefficient, while for cubes diagonal
to each other (c2,g2,02), the divergence only appears at the
fourth virial coefficient. This is plausible because depletion
forces are stronger when the area of contact of the large
particles is larger. Thus, in the first case, the depletion effect
is stronger.. In Table I, the values of the one-dimensional
integrals are listed that correspond to the graphs.
a.d5~/,Z)exp(P4) = 1 +PI A
+ Pz
A
2 T.
-q&J b +dsJ t’pq +.pg ).
+y-%n +dq +4&J +2&+2pq,
+cgn +4F;b +BB +pg)
+
. . . . . .
(2)
J. Chem. Phys., Vol. 101, No. 4, 15 August 1994
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp