Table Of ContentPhase diagram of a polydisperse soft-spheres model for liquids and colloids
L. A. Ferna´ndez,1,2 V. Mart´ın-Mayor,1,2 and P. Verrocchio2,3
1Departamento de F´ısica Te´orica I, Universidad Complutense, 28040 Madrid, Spain.
2Instituto de Biocomputaci´on y F´ısica de Sistemas Complejos (BIFI), Spain.
3Dipartimento di Fisica, Universit`a di Trento and CRS, SOFT, INFM-CNR, 38050 Povo, Trento, Italy.
(Dated: February 6, 2008)
7 Thephasediagramofsoftsphereswithsizedispersionhasbeenstudiedbymeansofanoptimized
0 Monte Carlo algorithm which allows to equilibrate below the kinetic glass transition for all sizes
0
distribution. The system ubiquitously undergoes a first order freezing transition. While for small
2
sizedispersion thefrozen phasehasacrystalline structure,large densityinhomogeneities appear in
n thehighlydispersesystems. Studyingtheinterplaybetweentheequilibriumphasediagramandthe
a kinetic glass transition, we argue that the experimentally found terminal polydispersity of colloids
J is a purely kineticphenomenon.
6
2 PACSnumbers: 64.60.Fr,64.60.My,66.20.+d
]
t Theequilibriumphasediagramofdenseclassicalfluids needed to keep the system in thermodynamic equilib-
f
o is far from being fully understood, especially as regards rium often become exceedingly slow. Such processes are
s
theinfluenceoverthefreezingtransitionofsomedisorder the nucleation of the solid within the fluid and the α-
.
t
a intheparametersoftheinteraction. Whilefluidsmadeof relaxation in the super-cooled fluid [9].
m identical atoms crystallize easily upon lowering the tem- Herewe obtainthe equilibrium phasediagramofpoly-
- peratureorincreasingthedensity,thefateofpolydisperse disperse soft-spheres in the (N,V,T) ensemble, aiming
d systems (colloids, for instance), where the particle size σ to describe both liquids and colloids. This is made pos-
n
is sampled from a probability distribution, P(σ), is still sible by the combination of an optimized Monte Carlo
o
a matter of debate. (MC) method (which, unlike standard MC, equilibrates
c
[ On the theoretical side, the effect of size dispersion, below the kinetic glass temperature) and a novel Finite-
measured by the quantity δ (the ratio among the stan- Size Scaling study of the fluid-solid coexistence line. At
2
v dard deviation and the mean of P(σ)) over the phase all δ, a first-order freezing transition separates the fluid
4 diagram has been addressed mostly in the case of the phase from the low temperature solid. This rules out
0 hard-spheresmodel, which is meant to describe colloidal the thermodynamic originof the terminal polydispersity
2 systems [1]. At least for small polydispersity, δ seems scenario. However, we show that a Brownian dynamics
9
to be the only feature of P(σ) that controls the physical will not find crystallization for δ > 0.12, in quantita-
0
6 results. Differenttheories predictconflicting scenariosin tive agreement with colloids experiments [7]. Further-
0 theregionoflargeδ. Themomentfree-energymethod[2] more, depending on polydispersity the solid phase can
/ suggests phase separation between many different crys- be eitheracrystaloraninhomogeneousphase(hereafter
t
a talphases,eachonewithamuchnarrowersizedispersion I-phase). The freezing temperature shows a reentrant
m
than δ (a phenomenon termed fractionation). However, behavior when increasing δ. Interestingly, in the range
- a density functional analysis [3] predicts the existence of δ [0.12,0.38]thekineticglasstransitionoccursfortem-
d ∈
n a terminal polydispersity δt beyondwhich the formation peratures, T, and densities, ρ, in the stable fluid phase.
o ofthecrystalisthermodynamicallysuppressed. Numeric Specifically,inoursimulationssoft-spheresofradiusσi
c simulations of the hard sphere model find some agree- and σj at distance r interact via the pair potential [30]
:
v ment with both the fractionation [4, 5] and the terminal
12
σ +σ
Xi polydispersity [6] scenarios. As regards models with a Vij(r)= i j . (1)
soft potential (e.g. Lennard-Jones), which are more ap- (cid:18) r (cid:19)
r
a propriate to describe liquids, no extensive study on the The effect of T and ρ is encoded in the single thermo-
effects of polydispersity has been performed so far. dynamic parameter Γ ρT−1/4. Although (1) general-
On the experimental side, crystallization of very vis- izes well known models≡for simple liquids [10], its scale-
couscolloidalsampleswithδ >δt 0.12doesnotoccur, invariantformsuggeststhatitcoulddescribeaswellcol-
≈
evenaftermonthsspentfromthepreparation[7]. Further loids,whosesizeisinthemicrometerrange. Weconsider
evidence supporting the terminalpolydispersity scenario aflatsizedistribution,constantintherange[σ ,σ ].
min max
comes from the finding of reentrant melting (crystal to In order to eliminate sample-to-sample fluctuations we
fluid transition when increasing the density) on polydis- follow [11], which for a flat distribution amounts to pick
perse colloids in confined geometry [8]. σ =σ +∆(i 1), with ∆=(σ σ )/(N 1).
i min max min
− − −
Yet, these experimental results do not necessarily re- Note that σmax/σmin at δ∞=1/√3.
→∞
vealfeaturesofthephasediagram. Indeed,theprocesses The numerical investigation of equilibrium properties
2
0.006
0.3 N=500 N=500
0.4 110067 I-phase CV 0.2 N=256 /N Fi0.004 N=256
0.3 τ105 0.1 h0.002
4
10
3 0 0
10
δ 1.52 1.54 1.56 1.52 1.54 1.56
0.2 1.4 1.45 1.5 Γ Γ
Γ 0.8
Fluid 20 0.6
0.1
f f 0.4
Crystal
10
0.2
0 0 0
1.2 1.3 1.4 1.5
0.68 0.72 0.76 0 1 2 3 4 5
Γ e
FIG.1: TheΓ−δ phasediagram showsthreephases: afluid is F
FIG. 2: Equilibrium quantities for δ=0.24 (top-left) Spe-
at high temperatures, a crystal, and a inhomogeneous solid
cificheat vs. Γ,from two simulations, oneat thesize depen-
(I-phase). TheboundarybetweenthecrystalandtheI-phase
dentcriticalpoint(fullsymbolsdenotetheactualsimulations,
lies at 0.18<δc<0.21. The vertical line in the fluid phase
full lines coming from histogram reweighting), the other at a
is the kinetic glass transition. (Inset) Integrated relaxation
higherΓ (opensymbols,anddashedlines). Themaximumof
time [12], τ, for e (full) and F (open), both for standard
the specific heat scales with N. (Top-right) hFi/N vs. Γ for
(squares) and local swap (circles) MC, for N = 256. The
kineticglasstransitionariseswhenτ∼106standardMCsteps. thesamesimulationsoftop–leftpanel. (Bottom-left)Normal-
izedhistogram,f,oftheenergyoftheinherentstructures,for
N=256 (empty) and N=500 (filled), at the Γ value where
of fluids is limited by the practical impossibility to equi- thepeaksof theinstantaneous energyhaveequalheight[18].
librate when either the relaxation time τ within the Inherentstructureswerefoundevery105 MCsteps. (Bottom-
fluid phase, or the freezing time t [31] become com- right) AsBottom-left, forF. InthefluidphaseF isO(1)(in
fr
parable with the time scale of the simulation. In or- the solid, F is O(N)). We also show (dashed line) eis and
F histograms for N = 864, where only two back and forth
der to achieveequilibrium we straightforwardlyadapt to
tunnelingeventsbetween theliquid and solid phase wereob-
model (1) the local swap MC algorithm [12, 13], which
served. Yet, data nicely agree with the predicted N-scaling.
outperforms [13, 14] other algorithms, such as the non-
local swap MC [15] and the density of states MC [16].
In fact, the standard method of studying a first-order forδ 0.21. However,for localswapMCtfr remainsbe-
≥
transitionin a systemof finite size, L, looksfor a double tween106 and107 MC steps for allδ, while for standard
peak in the histogramof the internal energy per particle MC it grows beyond 109 steps for δ >0.12.
e,accompaniedbyapeakofthespecificheatC [17,18]. Todetectthepossibleexistenceoflargedensityfluctu-
V
Thisprocedureisextremelydemandingonthe qualityof ationswefocus onS(q)≡ N1 i,jexp[iq·(ri−rj)](ri is
the statistical sampling, and has been scarcely used (if the positionof thei-thparticPle). Note that (q) is the
hS i
at all) for glass-forming liquids or colloids. Fortunately, static structure factor. The longest wavelengths fluctua-
the local swap MC allows us to employ it. tionsarestudiedthrough [ (2π,0,0)+ (0,2π,0)+
Thethermalizationissueneedstobeaddressedinthree (0,0,2π)]/3. Infact,inaFn≡homSogLeneousliqSuidoLrcrys-
S L
different regimes, see Fig. 1: the liquid phase, at the tal phase we expect to be of order 1, but for a macro-
F
phase coexistence line, and the solid phase. The swap scopically inhomogeneous phase it becomes of order N.
algorithm avoids the cage effect, thus equilibrating eas- In Fig. 2 we show as an example the results obtained
ily the whole liquid phase (it reduces τ by two orders of for δ=0.24 where an evidence for a freezing transition
magnitude as compared with standard MC, see Fig. 1– is presented. The specific-heat displays a peak of height
inset). Thermalization in the deep solid phase, not at- proportional to N while, as usual for small systems, its
temptedinthiswork,wouldrequireadifferentapproach. position suffers a strong finite size shift. We found con-
At the phase coexistence line, our criterion for thermal- venienttousethehistogramoftheenergyperparticleof
ization required the observation of tenths of back and the inherentstructureseIS,i.e. the minima ofthe poten-
forthtunnelingeventsbetweentheliquidandsolidphase. tial energy [19]. The advantages are twofold (Fig. 2): it
The difficulty for meeting the criterion grows dramat- largelyabsorbstheeffectsofthefinite-systemshiftofthe
ically with the number of particles ( exp[ΣN2/3], Σ critical temperature (so that the position of the peaks is
∼
being the liquid-solid surface tension). A stronger, more almostN independent)anditmakeseachpeaknarrower.
quantitative check is the consistency of the histogram Note that δ =0.24 is much higher than the terminal
reweighting [32]. We equilibrated N =256 particles for polydispersity δ reported in experiments and simula-
t
δ 0.1 (and for the monodisperse system) and N=500 tions [7, 11]. Actually, a freezing transition arises in the
≥
3
full range of δ studied. The line of phase-coexistence,
as located by the arising of a double peak for N =256, 4.0 δ=0.24, N=256 δ=0.24, N=500
between the fluid and the solid phase is shown in Fig. 1. 3.0
For δ=0, we find a body-centered cubic (bcc) crystal,
2.0
as expected for modest N [20]. Indeed, (q) displays
peaks of order N atq 2π(1,1,0),2π(0,1,S1),2π(1,0,1) 1.0
with lattice spacing a∼ a0.78. Theasame Braagg peaks 0.0
∼ 0.7 0.75 0.8 e 0.7 0.75 0.8
are found at δ=0.12, broadened due to disorder [33] (as is
F
compared with δ =0, the intensity reduces by a factor 1.0
δ=0.18, N=256 δ=0.18, N=500
1/4, but it still increases linearly with N).
Interestingly, for δ=0.24 the histogram of (Fig. 2)
developsadoublepeakatthe freezingpoint. LFeft-peak’s 0.5
position is N-independent, as expected for the liquid
phase (see also the scatter plot in Fig. 3). On the other
0.0
hand,thesecondpeak(correspondingtothe solid)shifts 0.55 0.6 0.65 0.7 0.55 0.6 0.65 0.7
e
right proportionally to N, revealing that the solid phase is
isnotastandardcrystallineone. ForN=500andδ .0.2 FIG. 3: Scatter plot of F vs. eIS below (bottom) and above
(top) the I-phase-crystal transition line at the liquid-solid
it is risky to use histogram techniques, as the tunneling
phasecoexistence,forN=256(left)andN=500(right). The
from solid to liquid is harder to observethan the reverse numberofpointsis∼45000 (N=256)or∼15000 (N=500).
liquidtosolidtunneling,thusweshowinFig.3ascatter Thedarkertheshadeofgray,thehigherthedensityofpoints.
plot. It is clear that in the crystals we found at δ=0.18
(corresponding to the lower values of eIS in Fig. 3) 0.1
F
takes a N-independent value [34]. Summarizing, at high
polydispersities we havefound a freezing transitionfrom f
the liquid to a solid inhomogeneous phase, while at low 0.05 δ=0.24
polydispersities the low temperature high density phase
(largeΓ’s)isastandardhomogeneouscrystal. Thestudy
of the transition between the crystal and the I-phase is 0
left for future work. To gain some insight about the I- 0 0.1 0.2
phase we measured the intensity of the Bragg peaks for 0.1
0.15
theinherentstructuresinthesolidphase. Moreprecisely,
f f
we define as the maximum [35] of (q). In a crystal 0.1
is of ordBer N (the corresponding q iSs in the reciprocal 0.05 δ=0.30 δ=0.45
B 0.05
lattice), while in a fluid is always of order 1. In Fig. 4
B
we show the histogram of along the phase coexistence
B 0 0
line both for N=256 andN=500. For everyδ we find a 0 0.1 0.2 0 0.05 0.1
double peak structure. The solid’s peak shifts right pro- B
FIG. 4: (Top-left) A δ = 0.3 solid configuration of 500
portionallyN. Thus and provideaclassification
particles, at phase coexistence, as projected into the XY
hBi hFi
of the Γ δ plane: distinguishes the solid from the plane (circles are centered at particles positions, their diam-
− hBi
fluid phase, while is of order N only in the I-phase. eter being proportional to particle sizes). Larger particles
hFi
The solid peak in Fig. 4 shifts left with increasing δ form thecrystalline central band(theBragg peaks’s analysis
(at δ = 0.45 the two peaks are hardly resolved). The suggests a FCC ordering). Normalized histograms of B for
δ=0.24,0.3,0.45 with 256 (empty) and 500 (filled) particles.
I-phase might signal either fractionation [4, 5] or phase
coexistence of a crystalwith anamorphoussolid (Fig. 4,
top-left). Further work is needed to clarify this point. case and 10−2 secs. in the latter one [21].
∼
Theexistenceofafreezingtransitionforallδ rulesout WehavelocatedthekineticglasstransitioninFig.1as
the terminal polydispersity scenario. But even if a sta- the point where τ reaches 106 (standard) MC steps. In
ble solidphaseexistsforlargeδ,itmightbe dynamically the case of colloids, this roughly corresponds to a relax-
inaccessible on experimental time-scales. In order to de- ation time of three hours, hence our kinetic glass transi-
tect the occurrence of the kinetic glass transition on the tion corresponds to the experimental one. On the other
δ Γ plane we adopted standardMC simulations which hand, since in the case of liquids 106 MC steps 10−7
− ∼
(atvariancewithlocalswapMC),mimictherealdynam- secs., our kinetic glass transition is rather close to the
ics of the system by modeling it as a Brownian motion. mode-coupling transition [22]. Indeed, for most molecu-
Note that the correspondence between the single step of lar and polymeric glass-formers τ at the mode-coupling
standardMCandthephysicaltimeunitswidelydifferfor transition lies between 10−7.5 and 10−6.5 secs. [23]. By
liquids and colloids, being 10−13 secs. in the former the way, since τ grows by many order of magnitudes in
∼
4
a narrow range, for liquids the difference in Γ between versity Press, 1997).
the experimental and the mode coupling transition will [10] J.P.HansenandI.R.McDonald,Theory of SimpleLiq-
not be larger than 5% [24]. In the inset of Fig. 1 we uids (Academic Press, San Diego, 1986).
[11] L. Santen and W.Krauth, condmat/0107459 (2001).
observe that the kinetic glass transition is at Γ 1.45
c
≃ [12] L. A. Fern´andez, V. Mart´ın-Mayor, and P. Verrocchio,
for δ=0.24. This value coincides with the one found for
Phys. Rev.E 73, 020501(R) (2006).
the soft-spheres binary mixtures (i.e. a bimodal P(σ))
[13] L. A. Fern´andez, V. Mart´ın-Mayor, and P. Verrocchio,
with δ=0.09 [25] and δ=0.16 [26]. This suggests that Philosophical Magazine 87, 581 (2006).
Γc is independent on δ. The vertical line at Γc =1.45 [14] L. A. Fern´andez, V. Mart´ın-Mayor, and P. Verrocchio
intersects the freezing transition line at two intersection (AIP,Melville, NewYork,2005), no.832in AIPConfer-
points,creatingasortofnoman’slandwhereequilibrium ence Proceedings Series, pp.128–133.
[15] T.S.GrigeraandG.Parisi,Phys.Rev.E63,045102(R)
is experimentally unreachable on human time-scales. In
(2001).
fact, the freezing time t depends on δ, growingtremen-
fr [16] Q. Yan, T. S. Jain, and J. J. de Pablo, Phys. Rev. Lett.
douslyalongthephasecoexistenceline. Whenapproach- 92, 235701 (pages 4) (2004).
ing the first intersection point (δ 0.12, see Fig. 1), for [17] M.S.S.Challa,D.P.Landau,andK.Binder,Phys.Rev.
≈
standard MC we only know that tfr becomes largerthan B 34, 1841 (1986).
109 steps. For δ > 0.12 our standard MC simulations [18] J. Lee and J. M. Kosterlitz, Phys. Rev. Lett. 65, 137
did not find the solid phase. We thus argue that the (1990).
[19] F.H.StillingerandT.A.Weber,Phys.Rev.A28,2408
tremendous slowing down both of τ and t makes the
fr
(1983).
solid phase for δ &0.12 experimentally unattainable.
[20] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel,
Westudiedtheequilibriumphasediagramintheδ Γ Phys. Rev.Lett. 75, 2714 (1995).
−
plane for a polydisperse soft sphere model with empha- [21] N. B. Simeonova and W. K. Kegel, Phys.Rev. Lett. 93,
sis on the first order freezing transition. Optimized MC 035701 (pages 4) (2004).
algorithms and Finite Size Scaling analysis were crucial [22] W. G¨otze and L. Sj¨ogren, Rep. Prog. Phys. 55, 241
to obtainequilibrium properties. We found twodifferent (1992).
[23] V. N. Novikov and A. P. Sokolov, Phys. Rev. E 67,
solid phases, a standard crystal at small δ and a highly
031507 (pages 6) (2003).
inhomogeneousphaseatlargeδ. Moreover,forδlyingbe-
[24] E. R¨ossler, A.P.Sokolov, A.Kisliuk,and D.Quitmann,
tween0.12 and 0.38the glasstransitionis no longerpre- Phys. Rev.B 49, 14967 (1994).
emptedbythefreezingtransition,makingthisrangevery [25] B. Bernu, J. P. Hansen, Y. Hiwatari, and G. Pastore,
promising for the theoretical study of glasses. In experi- Phys. Rev.A 36, 4891 (1987).
mentsnon-equilibriumeffectsshoulddominatethewhole [26] C. C. Yu and H. M. Carruzzo, Phys. Rev. E 69, 051201
regionδ &0.12,whereeitherthekineticglasstransitions (2004).
[27] M. Falcioni, E. Marinari, M. L. Paciello, G. Parisi, and
or the growth of the freezing time-scale are expected to
B. Taglienti, Phys. Lett.B 108, 331 (1982).
hidethefreezingtransition. Thisisinquantitativeagree-
[28] A. M. Ferrenberg and R. H.Swendsen, Phys. Rev. Lett.
mentwiththefindingsforcolloids[7]suggestingthatthe 63, 1195 (1989).
polydisperse soft-spheres model (1) could catch the fea- [29] H.J.Sch¨ope,G.Bryant,andW.vanMegen, Phys.Rev.
tures both of molecular glass-formers and of colloids. Lett. 96, 175701 (pages 4) (2006).
We thank V. Erokhin and F. Zamponi for discus- [30] We use the long distance cut-off of [13, 16]. Our length
sions. WeweresupportedbyMEC(Spain),throughcon- unit σ0 is fixed by σ03=R dσidσjP(σi)P(σj)(σi+σj)3.
Wesimulated N (initially fullydisordered)particles ina
tracts BFM2003-08532, FIS2004-05073, FPA2004-02602
box with periodic boundary conditions at ρ=σ−3.
and by BSCH—UCM. The CPU time utilized (at BIFI 0
[31] We define tfr as the characteristic time for the first tun-
andCINECA)amountsto10yearsof3GHzPentiumIV.
neling from theliquid to thesolid phase.
[32] The average for a N particles system at temperature
1/(β+δβ),ofanyfunctionA[{x }],canbeobtained[27,
i
28] from mean-values at temperature 1/β through the
identity hAiβ+δβ=hAexp[−Nδβe]iβ/hexp[−Nδβe]iβ, (e
[1] F.SciortinoandP.Tartaglia,Adv.Phys.54,471(2005). is the energy per particle), that, in the δβ → 0 limit,
[2] M.FasoloandP.Sollich,Phys.Rev.E70,041410(2004). yieldsthe(differential)Fluctuation-DissipationTheorem
[3] P. Chaudhuri, S. Karmakar, C. Dasgupta, H. R. Krish- ∂ hAi = N[hAei − hAi hei ]. Thus, the reweight-
β β β β β
namurthy,and A. K. Sood, Phys. Rev. Lett. 95, 248301 ingmethodexploitsan(integral)Fluctuation-Dissipation
(pages 4) (2005). Theorem, and consistency in the extrapolation (in prac-
[4] P.Bartlett, J. Chem. Phys. 109, 10970 (1998). tice, δβ.[NC ]−1/2) is a strong thermalization test.
V
[5] D. A. Kofke and P. G. Bolhuis, Phys. Rev. E 59, 618 [33] See [29] for recent experiments at low polydispersity.
(1999). [34] The eIS dispersion at δ=0.18 crystals in Fig. 3 is due
[6] S.Auerand D.Frenkel, Nature 413, 711 (2001). totheinterplaybetweenspatialorientation andperiodic
[7] P.N.PuseyandW.vanMegen,Nature320,340(1986). boundary conditions.
[8] R. P. A. Dullens and W. K. Kegel, Phys. Rev. Lett. 92, [35] Maximum over q=2Lπ(n1,n2,n3),|ni|≤20 and q6=0.
195702 (pages 4) (2004).
[9] P. G. Debenedetti, Metastable Liquids (Princeton Uni-
