Table Of ContentInfluence of polydispersity on the critical parameters of an effective potential model
for asymmetric hard sphere mixtures
Julio Largo1 and Nigel B. Wilding1
1Department of Physics, University of Bath, Bath BA2 7AY, United Kingdom
6 (Dated: February 4, 2008)
0
We report a Monte Carlo simulation study of the properties of highly asymmetric binary hard
0
sphere mixtures. This system is treated within an effective fluid approximation in which the large
2
particlesinteractthroughadepletionpotential(R.Rothetal,Phys. Rev. E625360(2000))designed
n to capture the effects of a virtual sea of small particles. We generalize this depletion potential to
a include the effects of explicit size dispersity in the large particles and consider the case in which
J the particle diameters are distributed according to a Schulz form having degree of polydispersity
6 14%. Theresultingalteration (withrespect tothemonodisperselimit) ofthemetastablefluid-fluid
critical pointparametersisdeterminedfortwovaluesoftheratioof thediametersofthesmall and
] large particles: q ≡ σs/σ¯b = 0.1 and q = 0.05. We find that inclusion of polydispersity moves the
h
critical point to lower reservoir volume fractions of the small particles and high volume fractions
c
of the large ones. The estimated critical point parameters are found to be in good agreement with
e
m thosepredictedbyageneralized correspondingstatesargumentwhichprovidesalinktotheknown
critical adhesion parameter of the adhesive hard sphere model. Finite-size scaling estimates of the
-
clusterpercolation lineintheonephasefluidregion indicate thatinclusion of polydispersitymoves
t
a the critical point deeper into the percolating regime. This suggests that phase separation is more
t
s likely to be preempted bydynamical arrest in polydisperse systems.
. PACS numbers: 64.70Fx, 68.35.Rh
t
a
m
I. INTRODUCTION AND BACKGROUND veloped novel algorithms [6, 7] offer some hope of future
-
d progress in accessing greater size asymmetries, no direct
n evidence has (to date) been obtained for the existence
Highlyasymmetricmixturesofhardsphereshavelong
o
of a fluid-fluid phase separation in additive hard sphere
servedasaprototypemodelforsystemsoflargecolloidal
c
mixtures.
[ particles dispersed in a sea of smaller colloids. A key
physical feature of such systems is the mediation by the In view of these difficulties, a fruitful alternative to
1
small particles of so-called depletion forces between the simulations of the full two component mixture is to at-
v
9 largeones[1]. This force hasits origininentropiceffects tempt to map it onto an effective one-component sys-
1 associatedwiththe dependence ofthe free volumeofthe tem which can be simulated more easily. This “effec-
1 small particles on the degree of clustering of the large tive fluid” approachwas takenby Dijkstra, vanRoijand
1 ones. Although, the typical range of the forces is rather Evans [3], andby Almarza amd Enciso[4] who proposed
0 limited(oforderthediameterofthesmallparticles),they a model depletion potential by tracing outfrom the par-
6
can be very strong. tition function the degrees of freedom associated with
0
A longstanding issue in this contextconcerns the abil- the small particles. The resulting interparticle potential
/
at ity of depletion forces to engender phase transitions in for the large particles is parameterized by the size ratio
m binary hard sphere mixtures. Biben and Hansen [2] ad- q =σs/σb betweensmallandlargeparticles,andareser-
- dressed this matter using integral equation theory, and voirvolumefractionofthesmallparticlesη˜s. Simulation-
d predicted that for sufficient size asymmetry, depletion based free energy measurement of the resulting system
n forces engender a fluid-fluid spinodal instability. Other yielded,forvaluesofq <0.1,afluid-fluidseparationthat
o theoretical studies have arrived at often conflicting con- wasmetastablewithrespecttoabroadsolid-fluidcoexis-
c
clusions in this regard (see discussion in ref. [3]), but tence region. Explicit simulations of the two component
:
v experiments on colloidalsystems (eg. ref. [5]), do appar- mixture confirmed the accuracy of the model depletion
Xi ently confirm a transition, although in certain circum- potential as far as the location of solid-fluid and solid-
stances it is found to be metastable with respect to a solidtransitions wasconcerned,but could notaccess the
r
a broad fluid-solid coexistence region. likely region of fluid-fluid separation. Subsequent work
Ideally one would like to settle the matter of the exis- by Roth, Evans and Dietrich [8] yielded a more accurate
tence of a fluid-fluid transition(as well as its stability or depletion potential by fitting to accurate DFT predic-
otherwise with respect to the fluid-solid boundary), by tions. Howeverto ourknowledgethe latter potentialhas
computer simulation. Unfortunately, direct simulation to date not been used to study phase behaviour.
studiesofveryasymmetricadditivemixturesareseverely Realcolloidalfluidsarepolydisperse,thatistheircon-
hamperedby extremely slowrelaxation. Accordingly,all stituent particles exhibit an essentially continuous range
such studies to date have been restricted to mixtures of of size, shape or charge. Introducing polydispersity into
relatively low asymmetries, for which a fluid-fluid tran- modelfluidsisknowntoalterfluid-fluidcriticalpointpa-
sition is less likely to be observable. While recently de- rameters[9,10]aswellasfreezingboundaries[11,12]. It
2
is therefore pertinent to enquire as to the effect of poly-
dispersity on the location of the metastable fluid-fluid
R +R
transition of hard sphere mixtures. Indeed this question W =ǫ b sW¯ . (1)
2R
haspreviouslybeenaddressedinpartbyWarren[13]who s
applied the moment free energy method [14] to study a
Here R and R are the radii of the large and small
b s
polydisperse version of the equation of state of Boublik spheresrespectively. ThescaleddepletionpotentialW¯ =
and Mansoori [15] for binary hard sphere mixtures. The W¯(x,η )isafunctionofx=h/σ ,thedistancefromcon-
s s
results showed that for sufficiently large size ratio and
tactmeasuredinunits ofthe smallsphere diameter,and
polydispersity of the large particles, a fluid-fluid spin-
thevolumefractionη ofthe smallspheres.. The param-
s
odal appearsin the model. The transitionwas predicted
eter ǫ takes the value ǫ = 2 for wall-sphere interactions
to become morestablewith increasingdegreeofpolydis-
and ǫ=1 for sphere-sphere interactions.
persity.
Between contact at x=0 and the location of the first
Most other theoretical investigations of polydispersity
maximum x the scaled depletion potential is expressed
0
inhardspheremixtures[16–18]havefocussedontheform
in terms of a cubic polynomial:
of the depletion potential and did not explicitly consider
the consequences for phase equilibria. The sole study
of phase behaviour to date (of which we are aware) is βW¯(x,η )=a(η )+b(η )x+c(η )x2+d(η )x3, x x .
s s s s s 0
≤
that of Fasolo and Sollich [11] who applied the moment (2)
free energy method to the Asakura-Oosawa(AO) model The coefficients a,b,cand d were obtained by Roth et al
[19]. This model describes the limit of maximum non- by fitting to depletionpotentials calculatedwithinDFT.
additivity of the small particles and in contrast to addi- For x > x , they assume that the asymptotic regime al-
0
tive mixtures, the monodisperse AO model is known to readysetsin. Fortheinteractionbetweentwospheres(a
exhibit a stable fluid-fluid phase transitionfor size ratios different expression applies to the asymptotic behaviour
q . 0.5 [20, 21]. The introduction of size polydispersity between a sphere and a hard wall) this is
to the large spheres [11] was observed to disfavor both
fluid-fluid and fluid-solid phase separation, though the
effect was larger for the latter transition. The net result βW¯(x,ηs)=βW¯asym(x,ηs), x>x0 (3)
was a lowering of the q value necessary for occurrence of
where
stablefluid-fluidphaseseparation,andanincreaseinthe
stability of this transition with respect to freezing.
In the present work we apply specialized Monte Carlo A (η )
simulation techniques to an effective fluid model for ad- βW¯asym(x,ηs) = σ−1(pσ s+h)exp(−a0(ηs)σsx) (4)
ditive mixtures with a view to elucidating the effect of s b
cos[a (η )σ x Θ (σ )], x>x
largeparticle(colloidal)polydispersityontheparameters 1 s s p s 0
× −
of the fluid-fluid critical point. Our results indicate that
Here the denominator measures the separation between
polydispersity shifts the critical point to lower reservoir
the centers of the spheres in units of σ , and the prefac-
volumefractionsη˜ ofthesmallerparticles,andtohigher s
s torsa (η )anda (η )canbecalculatedfromthePercus-
volumefractionsη ofthelargeones. Adeterminationof 0 s 1 s
b Yevick bulk pair direct correlationfunction [8]. The am-
the percolation threshold using finite-size scaling meth-
plitude A (η ) and phase Θ (η ) are chosen such that
ods shows that the addition of polydispersity moves the p s p s
thedepletionpotentialanditsfirstderivativearecontin-
fluid-fluid critical point deeper into the percolation re-
uousatx . Theyareweaklydependentonthesizeratio.
gion. We further find that accurate predictions for the 0
Fig. 1 shows the form of the potential for the size ratio
critical reservior volume fraction η˜crit can be obtained
s q =0.1 at a selection of values of η .
by matching the secondvirialcoefficientofthe depletion s
TheREDpotentialwastestedinref. [8]bycomparing
potentialtothatoftheadhesivehardspheremodelatits
with computer simulation results for hard sphere mix-
(independently known) critical point.
turesandwasfoundtoperformwellforvolumefractions
of the small particles in the range 0 < η < 0.3. In
s
the simulations to be described below, we use a trun-
II. MODELS
cated version of the RED potential, cutoff at x = 0.3,
and with no correction. Furthermore we shall employ
A. Depletion potentials the potential at finite volume fractions η of the large
b
particles, but assume the potential remains two-body in
Two model depletion potentials are considered in this form, being parameterised by the reservoir volume frac-
work. Thefirst,onwhichweshallfocusprimarily,isdue tion of small particles η˜ , which plays a role similar to a
s
to Roth, Evans and Dietrich [8], and we shall refer to it chemical potential. Clearly η˜ η in the limit η 0.
s s b
→ →
as the RED potential. It derives from accurate density We note that approximate expressions exist which allow
functionaltheory(DFT) studies ofhardspheremixtures one to convert from η˜ to η at finite η [3], at least in
s s b
and takes the form: the monodisperse case.
3
B. Incorporating polydispersity
2
Thekeytoincorporatingpolydispersityintotheabove
0 framework is to generalize the form of the depletion po-
tential for the case of two interacting large particles of
W(x)-2 η=0.24 daipffpeeraelnttoratdhiei DRe1rjaangduiRn2a.ppTrhoxisimisatrieoandi[l2y4,ac2h5i]e,vwedhibchy
s
ηs=0.26 relates the depletion force between two spheres of differ-
η=0.28
-4 ηs=0.30 ent radii (R1 and R2) to the potential between two flat
s plates:
-6
R R
1 2
0 0.1 0x.2 0.3 0.4 Fss =2πR +R Uww(h). (6)
1 2
FIG. 1: The form of the RED potential for q = 0.1 at a Theapproximationalsoyieldstheforcebetweenasphere
selection of values of ηs. (of radius R) and a wall:
F =2πRU (h) (7)
The second potential that we have studied, albeit to sw ww
a lesser extent and solely in the monodisperse context,
Clearly, if the two spheres in the first case have equal
is that due to G¨otzelmann et al [22, 23]. In contrast to
radii (R =R ), then
the DFT-based RED potential, this was derived purely 1 2
within the framework of the Derjaguin approximation,
although it too is expressed as a series expansion. We F =πRU (h), (8)
ss ww
shallemployitinthetruncatedformstudiedbyDijkstra,
vanRoijandEvans[3],andrefertoastheDREpotential: giving the well known result F =2F .
sw ss
Returningtothedepletionpotentialofeq.1,theabove
considerations prompt one to write:
1+q
βV (r ) = 3λ2η˜ +(9λ+12λ2)η˜2 (5)
eff ij − 2q s s ′ ′
R R
+(36λ+(cid:2)30λ2)η˜s3 ]; −1<λ<0 W = R1′ 1+R22′ W¯ (9)
where λ=x 1. where
−
Although both potentials have a qualitatively similar
R +R
′ 1 s
format shortrange,the DRE potential neglects the cor- R = (10)
1 R
rect damped oscillatory decay at larger particle separa- s
R +R
tions. A comparison of the two potentials for size ratio R′ = 2 s (11)
q =0.1 and η˜s =0.3 is shown in fig. 2. 2 Rs
It is readily verifiable that in the limiting cases R
1
2 →∞
and R R , one recovers the expression of Roth et al
1 2
→
[8] (eq. 1) with ǫ=2 and ǫ=1 respectively.
With regard to the effect of this generalization on the
0
parameterized form of the potential, we note firstly that
the quantity x remains unaffected because it is simply
x) the distance fromcontactinunits ofσs. However,inthe
W(-2 asymptotic part (eq. 4), the denominator σs−1(σb + h)
measuring the separation between the centers of the
DRE potential spheres in σ units is given in the polydisperse case by
s
-4 RED potential σ−1[(σ +σ )/2+h]. Additionally,whileinthemonodis-
s 1 2
perse limit the location of the first maximum of the po-
tential (at which the asymptotic behaviour is presumed
-6 to set in) is given by x = σ−1(σ +σ ), in the polydis-
0 0.1 0.2 0.3 0.4 0 s s b
x perse case one has instead x = σ−1[σ +(σ +σ )/2].
0 s s 1 2
Fig. 3 gives some examples of the influence of dissimilar
FIG.2: Comparison oftheREDandDREpotentials for q=
particle sizes on the depletion potential.
0.1 and η˜s =0.3.
In the presentwork,we addressthe situationin which
thelargeparticlesexhibitacontinuousvariationofsizes.
4
3 3
2.5
1
2
-1
) )
W(x η=0.30, σ=0.5, σ=0.5 σ(b1.5
s 1 2 f
-3 η=0.30, σ=1.0, σ=1.0
ηs=0.30, σ1=1.5, σ2=1.5 1
s 1 2
η=0.30, σ=1.5, σ=0.5
-5 s 1 2 0.5
-70 0.1 0.2 0.3 0.4 00 0.5 σ1 1.5 2
x
b
FIG. 3: The form of the RED depletion potential for four FIG.4: Theimposedformoff(σb),correspondingtoaSchulz
combinations of the pair radii. In all cases we have set σs = distribution (eq. 12) with z = 50. The diameters σb of the
0.1 large particles are measured in units of σ¯b =1
In order to quantify the form of the polydispersity, we
dimensional phase diagram [27]. Although this parame-
label each particle by the value of its diameter σ . The
b terizationprovidestheoperationalbasisforscanningthe
systemcanthenbedescribedintermsofadensitydistri-
dilution line, we shall (in accordance with convention)
bution ρ(σ ) measuring the number density of particles
b quote our results in terms of the overall volume fraction
of each σ . Experimentally, the distribution of colloidal
b ofthe largeparticles. The latter is relatedto the density
particles sizes in a system is (in general) fixed by the
distribution via
synthesis of the fluid. To reflect this situation in our
simulations, we assign ρ(σ ) an ad-hoc prescribed func-
b
∞
tionalform,whichwechoosetobeoftheSchulztype[26] π
η = σ3ρ(σ )dσ . (15)
defined by the normalized distribution function: b 6 b b b
Z0
Finally in this section, we note that within the poly-
f(σ )= 1 z+1 z+1σzexp z+1 σ . (12) disperse context, the distribution of sizes of the large
b z! σ¯ b − σ¯ b particlesimplies thatthe size ratioq canonlybe defined
(cid:18) b (cid:19) (cid:20) (cid:18) b (cid:19) (cid:21)
in terms of an average. Accordingly we take q σ /σ¯ .
s b
Here z is a parameter which controls the width of the ≡
distribution, while σ¯ 1 sets the length scale. We have
b
≡
elected to study the case z = 50, for which the corre- III. COMPUTATIONAL METHODS
spondingformoff(σ )isshowninfig.4. Theassociated
b
degreeofpolydispersityisdefinedasthenormalizedstan-
Monte Carlo simulations were performed within the
dard deviation of the size distribution:
grand canonical ensemble GCE using the methods de-
scribed in refs. [10, 28–30]. Here we briefly outline prin-
(σb σ¯b)2 21 cipal elements of the strategy and refer the interested
δ = h − i (13) reader to those papers for a fuller description.
σ¯
b Within the GCE framework, the density distribution
ρ(σ ) is obtained as an ensemble averageoveran instan-
FortheSchulzdistributiononefindsδ =1/√z+1. With b
taneously fluctuating distribution. The form of ρ(σ ) is
z = 50, this formula yields δ 14%. Note that for b
≈ controlled by the conjugate chemical potential distribu-
computationalconvenience,loweranduppercutoffswere
tion µ(σ ), which was tuned (cf ref. [30]) at all points
imposed on the range of allowed particle sizes σ . These b
b
in the phase diagram such as to yield the desired Schulz
were chosen such that 0.5 σ 1.5.
≤ b ≤ shape f(σ) (eqs. 12, and scale ρ0 14). This tuning was
Theimposeddensitydistributionisrelatedtof(σb)by b
achievedbyjointuseofthenon-equilibriumpotentialre-
finement (NEPR) method [29], coupled with histogram
ρ(σb)=ρ0bf(σb) (14) extrapolation [31] in terms of µ(σb). It should be noted,
however, that the DRE and RED depletion potentials
where ρ0 is the average number density of large parti- do not lend themselves to histogram extrapolation with
b
cles. Since f(σ ) is fixed, the form of ρ(σ ) is parame- respect to the model parameter η˜ which controls the
b b s
terizedsolely by ρ0, variations ofwhich correspond(at a formoftheinteractionpotential. Thisisbecauseη˜ does
b s
given η˜ ) to traversing a “dilution line” in the full finite not appear as an overallscale factor in the Hamiltonian,
s
5
a situation which contrasts, for example, to tempera- RED potential
ture reweighting in simpler potentials such as Lennard- q η˜scrit ηbcrit
Jonesium. Consequently in order to scan the phase dia- 0.1 0.3195(5) 0.274(10)
gramwith respect to changes in η˜s, separatesimulations 0.05 0.1765(5) 0.289(15)
were utilized in each instance.
Our principal aim is a determination of the
polydispersity-inducedshifts inthe criticalpoint param-
DREpotential
eters of the model depletion potentials. To this end we
have employed a crude version of the finite-size scaling q η˜scrit ηbcrit
0.1 0.255(15) 0.286(15)
(FSS) analysis described in ref. [32]. The analysis in-
volvesscanningtherangeofρ0 andη˜ untiltheobserved [0.289] [0.223]
b s
probability distribution of the fluctuating instantaneous 0.05 0.151(1) 0.271(15)
volumefractionoflargeparticlesp(ηb),matchestheinde- [0.165] [0.235]
pendently known universal fixed point form appropriate
to the Ising universality class in the FSS limit. Owing TABLE I: Estimates of critical point parameters, obtained
to the relatively large depth of the interparticle poten- using the methods described in the text. Values estimated
from the dataof ref.[3] are given in square brackets.
tial well (see fig. 2) compared to eg. the Lennard-Jones
(LJ)potential,theacceptancerateforparticleinsertions
and deletions was found to be very low, resulting in ex-
and the corresponding estimates for the DRE potential
tendedcorrelationtimeforthedensityfluctuations. Con-
are presented in tab. I.
sequently we were able neither to study a wide range of
systemsizesnorobtaindataofsufficientstatisticalqual-
0.6
ity to permit a more sophisticated FSS analysis. Never- η=0.3200
s
theless it transpires that our estimated uncertainties on η=0.3190
0.5 s
the critical point parameters are sufficient to resolve the Fixed point
Ising form
polydispersity-inducedtrendsinthecriticalpointparam-
0.4
eters that we set out to identify.
)
b
η 0.3
(
P
IV. SIMULATION RESULTS
0.2
Our results are divided into three sections. Firstly we
0.1
locate the fluid-fluid criticalpointfor boththe RED and
DRE potentials in the monodisperse limit. Moving on
0
to the polydisperse case, we determine the effect of the -3 -2 -1 0 1 2 3
a (η -<η >)
addedpolydispersityonthecriticalpointparameters. Fi- 0 b b
nally we use finite-size scaling to estimate the locus of
the cluster percolation threshold in both the mono- and FIG.5: Estimatesoftheorderparameterdistributionp(ηb)at
η˜s =0.3200andη˜s =0.3190,forV =(5.2σ¯)3. Representative
poly-disperse cases.
error bars are shown. Also included is the fixed point Ising
magnetisation distribution. All distributions are scaled to
unit norm and variancevia thenon-universalscale factor a0.
A. Monodisperse limit
Withregardtotheresults,oftab.I,wenotethatfora
1. Critical region givenpotentialform,the estimatesofηcrit appearrather
b
insensitive to the value of q. We further note that for a
As outlined above, we have tuned the values of η˜ given q there is a substantial shift in η˜crit between the
s s
and µ until the measured probability distribution of the two forms of the depletion potential. The latter finding
volume fractions of the large particles matched (as far is perhaps not too surprising given the significant differ-
as possible given the computational complexity of this ence in the contact value and range of the well depth of
problem) the universal Ising fixed point form. Fig. 5 the RED and DRE potentials, as well as the rather rad-
shows distributions obtained in this way for the case of ical truncation made by the DRE potential, of the long
the RED potential with q = 0.1 at η˜ = 0.3200 and ranged oscillatory part of the interactions (cf. fig. 2).
s
η˜ = 0.3190. Although the statistical quality is not par- Indeed the sensitivity of phase behaviour to the deple-
s
ticularly good, comparison of the forms of the distribu- tionpotentialwelldepthandrangehasbeenemphasised
tions with that of the fixed point form indicates that by Germain et al [33], albeit in the context of fluid-solid
the given values of η˜ straddle criticality, permitting the coexistence.
s
estimate η˜crit = 0.3195(5). This estimate for the RED We also note significant discrepancies between our es-
s
potential critical point, together for that for q = 0.05, timates of the critical point and those of Dijkstra et al
6
[3] for the DRE potential. Although no error bars are metastable for a period of time sufficient for us to col-
quotedinref.[3],itseemlikelytousthatthisdiscrepancy lect useful data in the critical region. Unfortunately, the
is statistically significant. Its source may be traceable to freezingbecameunmanageablewhenweattemptedtoob-
theusein[3]ofindirectfreeenergymeasurementstoob- tain data in the fluid-fluid coexistence region. As noted
tainthephasediagram,incontrasttothegenerallymore by other authors [3, 35] the coexistence curve of deple-
accurate direct grand canonical FSS approach employed tion potentials appears to be rather flat near the critical
here. Interestingly, our estimates for the critical η˜ val- point. Thus even a modest excursion into the two phase
s
ues lie much closer than those of ref. [3] to the results of regionresultsin highliquiddensities, whichinourexpe-
acomputationusingintegralequationtheoryofboththe rience froze very rapidly.
depletion potential and its phase behaviour [34].
0.6
B. Polydisperse case
0.5
Turning now to the polydisperse case, we have ob-
tained the critical point parameters in a manner simi-
0.4
lar to that employed in the monodisperse limit. Fig. 7
ηb0.3 shows the comparison of our estimates for the critical
point distribution p(η ) for the RED potential in both
b
the mono and polydisperse cases with q = 0.1. Clearly
0.2
the distributionfor the polydispersecaseis substantially
shifted to higher volume fractions compared to that for
0.1
themonodispersecase. Owingtotheslowfluctuationsof
η ,thetrueaverageofthesedistributionscouldnotbede-
0 b
0 1×108 2×108 3×108 4×108 5×108 termined to high precision. However, the peak positions
Monte Carlo Sweeps
are rather insensitive to the fluctuations, and on the ba-
sis of critical point universality one expects that given
FIG. 6: Time evolution of the system volume fraction near
sufficient statistics, the form of the distributions should
the parameters of the metastable critical point. Eventually,
becomesymmetric. Onecanthereforeestimateηcrit from
thesystem freezes spontaneously. b
theaverageofthepeakpositions. Theresultsofsodoing
aresummarizedintab.II,fromwhichone discoversthat
In ref. [3] it was demonstrated via free energy mea-
the principal influence of polydispersity on the critical
surements that the fluid-fluid critical point for the DRE point parameters is a significant decrease in η˜crit with
s
potential is metastable with respect to freezing. While
respect to its monodisperse value, and a significant in-
we have not attempted to perform a systematic study creasein ηcrit. Forthe form anddegreeof polydispersity
b
of the freezing transition in the present work, our simu- we have studied, the decrease in η˜crit is about 6%, while
s
lations confirm the metastability in so far as some runs the concomitant increase in ηcrit is about 17%.
b
were observed to freeze into an f.c.c. crystal structure.
An example of the time evolution of the density in such
5
a run for the DRE potential is shown in fig. 6. Such monodisperse η∼=0.320
freezing was also observed for the RED potential indi- polydisperse η∼=s0.300
s
cating that here too the criticalpoint is metastable with 4
respecttocrystallization. Asanaside,wenotethatfrom
a computational point of view, our observation of freez-
3
ing within the grandcanonicalensemble is somewhatre- )
b
η
markable since the algorithm becomes very inefficient at (
p
crystal densities. The key factor in achieving this in the 2
presentcaseistheinclusion–alongsidethestandardinser-
tions, deletions and resizing moves–of particle displace-
1
ment moves. Without the latter, the systemwas not ob-
served to crystallize on simulation time scales. Another
apparent factor controlling the ease of freezing appears 00 0.1 0.2 0.3 0.4 0.5 0.6
η
tobewhetherthecrystallatticeparameteriscommensu- b
rate with the choice of system box size. We further note
that our frozen structures do not attain the near-close FIG.7: Comparisonofthedistributionp(ηb)attheestimated
packing densities observed in ref. [3]. This is due to the criticalparametersinthemonodisperseandpolydispersesys-
presence of defects in the frozen configurations. tems. Inbothcasesthesizeratioq=0.1andthesystemsize
is V =(6σ¯)3
Notwithstanding the eventual relaxation to a crys-
talline state, our systems were usually found to remain
7
RED potential C. Percolation threshold
q η˜scrit ηbcrit
0.1 0.300(1) 0.336(15) Percolation is a necessary, though not sufficient con-
0.05 0.1655(5) 0.345(5) dition for gelation and dynamical arrest in colloidal sys-
tems. Sincegelationcanaffecttheabilityofexperiments
TABLEII:Estimatedcriticalpointparametersforthemodel to observe equilibrium phase behaviour in general, and
polydisperse system described in sec. IVB specifically fluid-fluid phase separation, it is important
to determine the location of the percolation line in the
phase diagramandits relationshipto the fluid-fluid crit-
ical point. Additionally, it is of interest to ask to what
The influence of polydispersity on the near critical
extent this relationship is affected by polydispersity.
point phase behaviour is further observable in terms of
In order to locate the percolation threshold, it is nec-
particlesize fractionationeffects. Specifically,when fluc-
essarytoidentify pairsofparticlesthatare‘bonded’and
tuationsoftheinstantaneousvaluevolumefractionη ex-
b
checkforspanningofclustersofsuchparticles. However,
ceedtheiraveragevalue,the distributionofparticlesizes
in contrastto lattice models or fluid systems such as the
is shifted to largerdiameters; and converselyfor fluctua-
adhesive hard sphere model, the definition of a ‘bond’
tions ofη to values lowerthan the average. The scaleof
b
in systems with continuous potentials is somewhat am-
the effect is shown in fig. 8 for q = 0.1 at the estimated
biguous. We therefore adopt a criterion which derives
critical point parameters. The presence of such fraction-
from that used for cluster identification in spin models.
ation implies that the critical point need not lie at the
Specifically, we determine the interaction energy u be-
apexofthecloudcurvethatmarkstheonsetofphasesep-
tween all pairs of particles and assigna bond with prob-
aration [27]. Indeed we did observe some evidence for a
abilityp =1 exp(βu). Clustersofbondedparticles
weakseparationofcloudandshadowcurvesatη˜s =η˜scrit, are thenboinddentifie−d using the algorithm of Hoshen and
although precise quantification of the effect was compli-
Kopelman [36]. In ref. [37] Miller and Frenkel identi-
cated by a noticibly increased tendency of the polydis-
fied the percolation threshold with those values of the
perse system to relax to a high density state following a
model parameters for which the proportion of configu-
fluctuationtohighdensity. Thenatureofthisrelaxation
rations containing a spanning cluster is 50%. However,
closelyresembledthatobservedinthemonodispersecase
finite-size scaling arguments [38] show that better esti-
(cf. fig. 6) and indeed visualization of the arrested con-
mates may be obtained by examining the finite-size be-
figurations revealed some evidence of crystalline order,
haviour of plots of the fraction of spanning clusters as
albeit with a high concentration of defects. We caution,
a function of η . An example of such a plot is shown
however,thatthese findingsshouldnotbeinterpretedas b
in fig. 9 for the RED potential in the monodisperse case
providing strong evidence for a freezing of the polydis-
for η˜ = 0.28. Data are shown for 4 system sizes, and
perse system because it was not possible to ensure that s
indicate that there is a well defined intersection point at
the overalldensity distribution remained on the dilution
η 0.21. This intersection point provides a good mea-
line during the relaxation to the high density state. b ≈
sure of the percolation threshold in the thermodynamic
limit [38]. By contrast, application of the 50% criterion
2 to data for a single system size can considerably overes-
timate the percolation threshold.
parent distribution _
size-distribution ηb<_ηb Percolation lines were determined using this intersec-
size distribution η>η
1.5 b b tionmethodforthemonodisperseandpolydisperseRED
potential at q = 0.1. They are shown in fig, 10 together
with our estimates of the critical point parameters. One
σ) 1 sees that in both cases the critical point lies well within
(
f the percolation regime, though much more so for the
polydisperse system than for the monodisperse system.
0.5
V. LINKING TO THE ADHESIVE HARD
0 SPHERE MODEL
0.5 0.75 1 1.25 1.5
σ
Asimpleyetgeneralmethodforfindingtwopotentials
FIG. 8: The normalized distribution of particle sizes for
that are ‘equivalent’ in a corresponding states sense, is
instantaneous volume fractions ηb below (dashed line) and
to match their second virial coefficient B [37, 39, 40]:
above (dotted line ) the average value, close to the critical 2
point. Also shown (solid line) in the overall Schulz “parent”
distribution ∞
B = 2π e−βu(r) 1 r2dr. (16)
2
− −
Z0 (cid:16) (cid:17)
8
1
Monodisperse RED potential
V=(9σ)3 q η˜s predicted η˜s simulation
b
0.8 V=(7.5σ)3 0.1 0.320 0.3195(5)
n b
o V=(6σ)3 0.05 0.177 0.1765(5)
acti0.6 V=(4.5bσ)3 Monodisperse DREpotential
fr b
n q η˜s predicted η˜s simulation
o
olati0.4 0.1 0.256 0.255(15)
c 0.05 0.151 0.151(1)
r
e
P Polydisperse RED potential
0.2
q η˜s predicted η˜s simulation
0.1 0.310 0.300(1)
0
0 0.1 0.2 0.3 0.4 0.05 0.172 0.1655(5)
η
b
TABLE III: Comparison of the simulation estimates of the
FIG. 9: Fraction of percolating configurations as a function critical point parameters with the predictions arising by
of ηb for the RED potential in the monodisperse limit. The matching the second virial coefficient to that of the critical
potential parameters are η˜s = 0.28, q = 0.1, and data are AHSmodel.
shown for four system sizes.
The AHS model exhibits a fluid-fluid phase transi-
tion, the critical point of which has been estimated to
0.32
occur [42] at τ = 0.1133(5),ρ = 0.508(10). This
percolation c c
(monodisperse) monodisperse c.p. value of τ implies that for the AHS model at criticality,
c
0.3 Bcrit = 4.826v with v = (4/3)πσ3. It is therefore of
2 − 0 0
polydisperse c.p. interestto assesswhether, via the matchingof B values
2
∼ηs0.28 for the RED and DRE potentials to that of the critical
AHS model, reasonable predictions can be made for the
critical point parameters of the depletion potentials. To
0.26 percolation
(polydisperse) thisendwehavenumericallyevaluatedB acrossarange
2
of η˜ values for each depletion potential and q value of
s
0.24 interest. By so doing we could determine that value of
η˜ for which B matched the value Bcrit = 4.826v .
0.1 0.15 0.2 0η.25 0.3 0.35 0.4 Tsable. III show2s the resulting predicti2ons for−η˜crit fo0r
b s
two values of q, together with our simulation estimates.
FIG.10: Percolation linefor theREDpotential(q=0.1) for Clearly, in each instance, the agreement is remarkable.
boththemonodisperseandpolydispersecases,asdetermined One can attempt to extend the above approachto the
by the method described in the text. Statistical errors are polydisperse depletion potentials. To do so, we first ob-
comparablewiththesymbolsizes; linesareguidestotheeye. tain the contribution to the second virial coefficient for
ThesystemsizeinbothcaseswasL=(6σ¯)3. Alsoshownare
interactionsbetweenallpairsofspecies σ ,σ . The over-
i j
theestimated critical point parameters (cf. sec. IIB).
all coefficient for the mixture can then be approximated
asaweightedaverageofpairs[43],wherethe weightfac-
tor is the probability of interaction between a pair of
HerewecomparethevalueofB fortheREDandDRE
2 particles of size i and size j. In our case,this is givenby
depletionpotentialsinthe monodisperselimit, withthat
theproductofthecorrespondingvaluesofthenormalized
oftheadhesivehardsphere(AHS)model[41]. Thelatter
Schulz distribution (eq.12):
compriseshardparticleswhichexperienceafiniteattrac-
tion only at contact, the strength of which is controlled
∞ ∞
via a ‘stickiness’ parameter τ. The overall interaction
B = f(σ )f(σ )B (σ ,σ )dσ dσ (19)
can be written 2 i j 2 i j i j
Z0 Z0
Matching to Bcrit = 4.826v as before, one obtains
σ 2 − 0
e−βu(r) =Θ(r σ)+ δ(r σ), (17) for the two q values studied, the predictions for η˜scrit
− 12τ − shown in tab III. Here the agreement with the simula-
tionestimatesislessimpressivethaninthemonodisperse
withrtheseparationofparticlecentersandσtheparticle
case. Although the absolute value of the prediction still
diameter. The second virial coefficient follows as
agrees to within about 3% with the simulation estimate,
andthe signofthepolydispersity-inducedshiftinη˜crit is
s
2π 1 correctly predicted, its magnitude is underestimated by
BAHS = σ3 1 . (18)
2 3 − 4τ a factor of two.
(cid:18) (cid:19)
9
The larger relative discrepancy between the predicted (and hence also the polydisperse depletion potentials),
andmeasuredη˜crit maypointtoabreakdowninthepres- as B = 4.826 5.759/6.621 4.2v . The resulting
s 2 − × ≈ − 0
ence of polydispersity of the assumed model invariance predictions for the critical reservoir volume fraction of
of the critical B value. Indeed one might expect such a the small particles are η˜crit = 0.304, for q = 0.1, and
2 s
failurebecausethevalueofB isbasedsolelyonthepair η˜crit = 0.169 for q = 0.05, both of which agree within
2 s
potential and takes no account of the ability of a poly- error with our simulation results. Of course in view of
disperse fluid to exploitlocal size segregationin order to our assumptions, this accordmay be fortuitous, but it is
packmoreeffectivelythanacorrespondingmonodisperse nevertheless suggestive that a more detailed assessment
one. In order to address this issue directly, one would of the effect of polydispersity on critical point B values
2
require estimates of critical point B values for a poly- in other systems (specifically the AHS model) would be
2
disperseversionoftheAHSmodel. Toourknowledgeno a worthwhile avenue for further study.
simulationestimatesoftheliquid-gastransitioncurrently
exist for a polydisperse AHS model. Indeed, the matter
is complicated by the fact that there is no unique model VI. CONCLUSIONS AND DISCUSSION
for polydispersity in such a system. Very recently, how-
ever, a number of physically reasonable models for poly-
To summarize, we have determined the effect of intro-
dispersityinAHSsystemhavebeenproposedbyFantoni
ducing polydispersity on the fluid-fluid critical point pa-
et al [44], who investigated the corresponding phase be-
rametersof a model depletion potential for highly asym-
haviour using integral equation theory. From ref. [44],
metric additive hardsphere mixtures. For the particular
one can deduce that the presence of polydispersity sig-
realization of the polydispersity considered, the critical
nificantly decreases the magnitude of B at the critical
2 pointis foundto shift (with respectto the monodisperse
pointcomparedtothemonodisperselimit. Thistrendin
limit)tosmallervaluesofthereservoirvolumefractionof
B is of the correct sign and overall magnitude to push
2 the small particles η˜ and to larger values of the volume
thepredictionsforη˜crit forourdepletionpotentialscloser s
s fraction of the large particle η . It seems reasonable to
tothe simulationestimates. Unfortunatelysincenodata b
assumethatthedirectionoftheshiftsshouldbeageneral
were reported for exactly the same degree of polydisper-
trend, common to other functional forms and degree of
sity (δ = 14%) studied in the present work, no direct
polydispersity. Indeed the same trend has also recently
comparison of B values is possible.
2 beenobservedinastudyofcolloidalpolydispersityinthe
In an attempt to throw additional light (albeit indi-
AO model [11].
rectly) on the discrepancy between the measured and
Beyondthis,ourresultsshowthatinclusionofpolydis-
predicted critical point parameters, we have studied the
persity pushes the whole fluid-fluid binodal deeper into
effect of introducing polydispersity on the critical point
the percolating regime. Since colloidal fluids are known
B value for the Lennard-Jones fluid, which is a com-
2 to form a gel [46, 47] for sufficiently high η˜ , it would
putationally more tractable systemthan the AHS model s
seem that the presence of polydispersity increases the
[45]. The corresponding potential is
likelihoodthatdirectobservationsoffluid-fluidphaseco-
existence is complicated by dynamical arrest.
Additionally we demonstrated that excellent predic-
σ 12 σ 6
uij =ǫij ij ij (20) tions for the value of η˜scrit follow from matching the sec-
"(cid:18)rij(cid:19) −(cid:18)rij(cid:19) # ond virial coefficient B2 of depletion potentials to the
criticalB value ofthe adhesivehardspheremodel. The
2
with ǫij = σiσjǫ, σij = (σi +σj)/2 and rij = ri rj . quantitative accuracy of the predictions is undoubtedly
| − |
The potential was cutoff for rij >2.5σij and no tail cor- due in large part to the very short ranged nature of the
rections were applied. For the monodisperse limit, the depletion potentials; similar studies comparing critical
critical temperature occurs at Tc = 1.1876(3) [32] and point B2 values for a range of other potentials [39, 40]
one finds B = 6.621v , which lies within the range of didnotfindsuchahighdegreeofaccuracy. Nevertheless
2 0
−
‘typical’ critical point B values found in surveys of a our observation should prove generally useful in reduc-
2
wide range of model potentials [39, 40]. If, on the other ing the effort required to locate criticality in depletion
hand,σ isdistributedaccordingto aSchulz form(eq.12) potentials. It is intriguing however, that the accuracy of
with z =50, as used elsewhere in this work, simulations thepredictionswasreducedonincorporatingpolydisper-
yield a critical temperature of Tc = 1.384(1), for which sity, suggesting that (perhaps due to changes in packing
(using eq. 19), one finds B2 = 5.759v0, which is sig- abilityduetolocalsizesegregationeffects),theinclusion
−
nificantly smaller in magnitude that the monodisperse of polydispersity in a model does not leave B invari-
2
value. ant at the critical point. Comparisons of the critical B
2
If one now makes the (not unreasonable) assumption value for a monodisperse and polydisperse LJ fluid con-
that given the same form and degree of polydispersity, a firmed a significant difference in this regard. Moreover,
comparablefractionalchangeofB willensueinotherin- the magnitude of the effect was sufficient to account for
2
teraction potentials, one can estimate the expected crit- the discrepancy in the observed and predicted η˜crit for
s
ical point value of B for the polydisperse AHS model the polydisperse depletion potential. Clearly, however,
2
10
further work is called for in order to elucidate this mat- text [48, 49].
ter more fully. Finally,weremarkthatwhilethepresentworkhascon-
Obviously knowledge of the shift in the critical point sidered solely the case of polydispersity of the large par-
parameters is in itself insufficient to determine whether ticles, the converse situation of small particles polydis-
polydispersity renders the fluid-fluid transition stable persity is clearly of interest and practical relevance too.
with respect to fluid-solid coexistence. Although we did While there have been several theoretical studies which
observeaspontaneousrelaxationofthenearcriticalpoly- have considered this scenario (see eg. refs. [17, 50]), we
disperse system to a high density state showing some know of no simulation studies to date. An extension
crystallineorder,thisfindingshouldbetreatedwithcau- of the methods utilised here to address this case would
tion because the system departs from the dilution line doubtless be a worthwhile endeavor–one which we hope
during the formation of the new state. It would thus be to undertake in future work.
interesting in future work to try to study explicitly the
effects of polydispersity on the freezing transition. As
wellasprovidingassessmentoftheoverallstabilityofthe
fluid-fluid critical point, freezing in polydisperse fluids is
Acknowledgments
amatter ofconsiderableinterestinits ownright. Indeed
recenttheoreticalcalculationsforthe AOmodelindicate
an increasing richness of fluid-solid and solid-solid phase The authors are grateful to R. Evans and A. Louis
behaviourasthedegreeofpolydispersityisincreased[11]. forhelpfuldiscussionsandadvice,andtoM.Dijkstrafor
Totacklethiscomputationally,however,isaconsiderable supplyingdatafromref.[3]. Thisworkwassupportedby
challenge, but one which might be met by extending to EPSRC, grant number GR/S59208/01. JL acknowledge
polydispersesystemnovelcomputationalmethodswhich support from Spanish Ministerio de Educacion, Cultura
havehitherto only be deployedinthe monodisperse con- y Deporte.
[1] For a general overview of colloidal interactions, see L. [18] D.FrydelandS.A.Rice,Phys.Rev.E71,041403(2005).
Belloni, J. Phys.: Condens. Matter, 12, 549 (2000). [19] S. Asakura and F. Oosawa, J. Chem. Phys. 22 1255
[2] T. Biben and J.P. Hansen, Phys. Rev. Lett. 66, 2215 (1954).
(1991). [20] M. Dijkstra, J.M. Brader and R. Evans, J. Phys. Con-
[3] M. Dijkstra, R. van Roij, R. Evans, Phys. Rev. Lett., dens. Matter 11, 10079 (1999).
81,2268(1998); M.Dijkstra,R.vanRoij,andR.Evans, [21] R.L.C.Vink,andJ.Horbach,J.Phys.: Condens.Matter,
Phys. Rev. Lett., 82, 117 (1999); M. Dijkstra, R. van 16, 3807 (2004); R.L.C. Vink and M. Schmidt, Phys.
Roij, and R.Evans, Phys.Rev. E, 59, 5744 (1999). Rev. E71, 051406 (2005).
[4] N.G. Almarza and E. Enciso, Phys. Rev. E59, 4426 [22] B.G¨otzelmann,R.Evans,andS.Dietrich,Phys.Rev.E,
(1999). 57, 6785 (1998).
[5] W.C.K. Poon, J. Phys.: Condens. Matter, 14, 859 [23] B. G¨otzelmann, R. Roth, S. Dietrich, M. Dijkstra, and
(2002). R. Evans, Europhys.Lett., 47, 398 (1999).
[6] A. Buhot and W. Krauth, Phys. Rev. Let 80, 3787 [24] B.V. Derjaguin. Kolloid-Z. 69, 155 (1934).
(1998). [25] W.B. Russel, D.A. Saville, and W.R. Schowalter, Col-
[7] J. Liu and E. Luijten, Phys. Rev. E71, 066701 (2005). loidal Dispersions, Cambridge University Press (1989)
[8] R. Roth, R. Evans, and S. Dietrich, Phys. Rev. E, 62, [26] G.V. Schulz. Z. Physik. Chem. B43, 25 (1939); B.H.
5360 (2000). Zimm, J. Chem. Phys. 16, 1099 (1948).
[9] L. Bellier-Castella, H. Xu and M. Baus, J. Chem. Phys. [27] P. Sollich. J. Phys. Condens. Matt. 14, R79 (2002).
113, 8337 (2000). [28] N.B Wilding and P. Sollich, Europhys. Lett. 67 219
[10] N.B. Wilding, M. Fasolo and P. Sollich, J. Chem. Phys. (2004).
121, 6887 (2004). [29] N.B. Wilding, J. Chem. Phys. 119, 12163 (2003).
[11] M. Fasolo and P. Sollich, J. Chem. Phys. 122, 074904 [30] N.B. Wilding and P. Sollich, J. Chem. Phys. 116, 7116
(2005). (2002).
[12] D.A. Kofke and P.G. Bolhuis, Phys. Rev. E 54, 634 [31] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett.
(1996); ibid59, 618 (1999). 61, 2635 (1988); ibid63, 1195 (1989).
[13] P.B. Warren, Europhys.Lett., 46, 295 (1999). [32] N.B.Wilding,Phys.Rev.E52602(1995); N.B.Wilding
[14] P. Sollich and M.E. Cates, Phys. Rev. Lett. 80 1365 andA.D.Bruce.J.Phys.Condens.Matt.4,3087(1992).
(1998); P.B. Warren, Phys. Rev.Lett. 80 1369 (1998); [33] Ph. Germain, J.G. Malherbe and S. Amokrane, Phys.
[15] T. Boubl´ik, J. Chem. Phys. 53, 471 (1970); G.A. Man- Rev. E70, 041409 (2004).
soori. N.F. Carnahan, K.E. Starling, T.W. Leland, J. [34] J.Cl´ement-Cottuz,S.Amokrane,andC.Regnaut,Phys.
Chem. Phys.54, 1523 (1971). Rev. E61, 1692 (2000).
[16] Y.Mao, J. Phys.II, 5, 1761 (1995). [35] F. Lo Verso, D. Pini, and L. Reatto, J. Phys.: Condens.
[17] D. Goulding, and J.-P. Hansen, Mol. Phys., 99, 865 Matter, 17, 771 (2005).
(2001). [36] J. Hoshen, and R. Kopelman, Phys. Rev. B, 14, 3438
