Model Energy Landscapes of Low-Temperature Fluids: Dipolar Hard Spheres

0.16 MB·English

by  Dmitry V. Matyushov

Model Energy Landscapes of Low-Temperature Fluids: Dipolar Hard
Spheres
∗
Dmitry V. Matyushov
7
Department of Physics and Astronomy and Center for Biophysics,
0
Arizona State University, PO Box 871604, Tempe, AZ 85287-1604
0
(Dated: February 6, 2008)
2
An analytical model of non-Gaussian energy landscape of low-temperature fluids is developed
n
a based on the thermodynamics of the fluid of dipolar hard spheres. The entire excitation profile of
J theliquid,fromthehightemperaturestothepointofideal-glasstransition,hasbeenobtainedfrom
the Monte Carlo simulations. The fluid of dipolar hard spheres loses stability when reaching the
3
point of ideal-glass transition transforming via afirst-order transition into a columnar liquid phase

] of dipolar chains locally arranged in a body-centeredtetragonal order.
t
f
o PACSnumbers: 64.70.Pf,65.40.Gr,66.20.+d
s
.
t
a Thermodynamic and dynamic properties of super- connection between S (T) and Sφ(T) is, however, still
c c
m cooled liquids are often related to the excess of config- unclear.
- urationstheypossessrelativetothecrystallinephase[1].
Analytical modeling of the landscape of supercooled
d
The number of states of the liquid state is reflected by
n liquids is to a large extent based on the ideas advanced
the configurationalentropyS (T)conventionallydefined
o c for spin glasses with quenched disorder [13, 14]. The
[2, 3, 4] as the logarithm of the density of states Ω(E)
c most prominent role in application to structural glasses
[ takenatthe energyequalto the internalenergyE(T) at
is played by Derrida’s random-energy,or Gaussian land-
temperature T:
1 scape,model(REM)[15]. Computersimulations,mostly
v limited to temperatures above the mode-coupling criti-
S (T)=ln[Ω(E(T))]. (1)
6 c cal temperature, basically support the REM, in particu-
5 lar the prediction of the 1/T falloff of the average basin
0 TheenergyintegralofthedensityofstatesandtheBoltz-
energy φ(T) from its high-temperature plateau [7, 8].
1 mann factor provide the canonical configurationintegral
However, general theoretical arguments [6, 16] and re-
0
7 centsimulations[12]suggestthatthe low-energyportion
0 Z(β)= Ω(E)e−βEdE, (2) of the Gaussian landscape might be inaccurate. In par-
Z
/ ticular,combinatorialargumentssuggestthatthederiva-
at tiveoftheenumerationfunction,σ (φ)=N−1ln[Ω (φ)],
m where β is the inverse temperature. φ φ
approaches infinity at the point of ideal glass transi-
Stillinger [5, 6] gave an alternative definition of the
- tion when the system runs out of configurations and
d configurational entropy in terms of the density of states
σ (φ )= 0 (N is the number of liquid particles). This
n Ω (E)ofinherentstructures,i.e.minimaofthetotalen- φ IG
φ infinite derivative eliminates the ideal glass transition at
o ergyofthesystemasafunctionofitsalltranslationalan
a positive temperature [6]. A phenomenological descrip-
c rotational degrees of freedom (energy landscape). The
: tion, patching together logarithmic and Gaussian enu-
v canonicalconfigurationintegralisnowgivenbyintegrat-
meration functions, was suggested to provide a correct
i ing over the basin depths φ and the free energy of inter-
X low-energy portion of σ (φ) [17].
basin vibrations F (φ), which is a weak, approximately φ
v
ar linear, function of φ [7, 8] Generally, unlike the case of spin glasses, there has
been a lack of simple, solvable models of structural
glasses. In this Letter we propose a new model of non-
Z(β)= dφΩ (φ)e−βφ−βFv(φ). (3)
φ Gaussian landscape of low-temperature fluids with no
Z
quencheddisorder. Ourderivationis basedonthe estab-
Because of the central role of the configurational en- lishedthermodynamicsofthemodelfluidofdipolarhard
tropyinearlytheoriesofviscousliquidsbyAdam,Gibbs, spheres (DHS) [18]. This monoatomic glass former (in
andDiMarzio[9,10],andmorerecentrandomfirst-order contrastto binary mixtures used in many recent simula-
transition models by Wolynes and co-workers[11], much tions[7,8])islongknowntoresistphasetransformations.
effort has been invested in calculating the Stillinger con- This fluid lacks the liquid-vapor transition decomposing
figurational entropy, Sφ(T) = ln[Ω (φ(T))], from com- instead into low-density dipolar chains [19, 20]. Upon
c φ
puter simulations of model fluids [7, 8, 12]. A general cooling, dipolar hard spheres transform into a ferroelec-
tric fluid [21], but the ferroelectric phase is stabilized by
tin-coilboundaryconditions usedinsimulationsandcan
be prevented by using a low dielectric constant for the
∗E-mail:[email protected]. reaction-fieldor Ewaldcorrectionsfor the cutoff ofdipo-
2
lar interactions [22]. This strategy has been employedin 5
(a)
this study. We have carried out the Monte Carlo (MC) 4
simulations of the DHS fluid within the standard NVT e)3
Metropolis protocol and the reaction-field cutoff correc- σ(2
tions. The reaction-field dielectric constant of ǫ =10,
RF 1
below the lowest value 18 permitting the ferroelectric
≃ 0
phase [22], has allowed us to eliminate the transition to
(b)
15
liquid ferroelectric. The initial configuration was set up
as a face-centered cubic lattice of 108 dipolar spheres. e)10
The cubic simulation box also helps to suppress crystal- P(
lization of dipoles which do not favor highly symmetric 5
lattice structures [23, 24]. The data were collected for
0
(2 10) 107 MC cycles. -2 -1.5 -1 -0.5 0
− × e
Apart from avoiding crystallization, a significant ad-
vantage of the DHS fluid is the existence of a simple
FIG.1: Enumerationfunctionσ(e)(a)anddistributionfunc-
analytical form for the free energy βF(β) = ln[Z(β)]. tion P(e) (b) of the fluid of dipolar hard spheres at ρ∗ =0.8
−
The Pad´etruncation of perturbation series suggested by and the temperatures T (from left to right in (b)): 0.125,
Stell and co-workers [25] turned out to be very success- 0.167,0.2,0.25,0.33,0.5,1.0,2.0,10.0. Thesimulationpoints
fulwhentestedagainstsimulations[18]. The free energy in (a) are obtained from P(e) within 95% of the maximum
of dipolar hard spheres depends on two parameters, the probability at each temperature.
reduced density ρ∗ = ρσ3 and the reduced temperature
T = k Tσ3/m2. Here, ρ is the number density, σ is the
B
hard-sphere diameter, and m is the dipole moment. All σ(e)=N−1ln[Ω(e)]:
calculationsandsimulationsherehavebeendoneatcon-
stant volume with ρ∗ =0.8 thus reducing the number of σ(e)=f(η) b−1 √c+e √c 2. (6)
variables to one. − −
(cid:0) (cid:1)
Stell’s Pad´e solution for F(β) is
The enumeration function in Eq. (6) is Gaussian near
itstop,e 0,wherethedipole-dipoleinteractionsarein-
Aβ2 ≃
βF(β)=Nf(η) . (4) significantandthelimitofahard-spherefluidisreached.
− 1+bβ It deviates from the parabolic shape in its low-energy
wing, in particular close to the low-energy cutoff at
Here f(η) is the reduced free energy of the fluid of hard
e = c (Fig. 1a). Depending on the parameters, Eq.
0
spheres with zero dipole moment as a function of the −
(6) gives two possible resolutions of the entropy crisis
∗
packing density η =(π/6)ρ . For the Carnahan-Starling
of low-temperature fluids [1]. When f(η) c/b 0, the
equationofstate,onhasf(η)=(4η 3η2)/(1 η)2. β = − ≥
− − enumerationfunctionhasaninfinitederivativeate0and,
1/T inEq.(4)isgiveninreducedunitsandthusthecoef-
according to Stllinger’s arguments [6], there is no ideal
ficients A=Na and b depend only on the liquid density
glass transition. When f(η) c/b is strictly positive,
through two-particle and three-particle perturbation in- −
the entropy is non-zero at the cutoff energy, a situation
tegrals, a = (ρ∗/6)I(2)(ρ∗), b = (ρ∗/9)I(3)(ρ∗)/I(2)(ρ∗),
observed for network glass formers [26]. Finally, when
tabulated by Larsen et al [25].
f(η) c/b<0,the idealglassarrestσ(e )=0happens
IG
Following Freed [4], the density of states can be ob- −
beforethecutoffisreached,andthetemperatureofideal-
tained by inverse Laplace transformation of Eq. (2) in
glasstransitionT ,e =e(T )ispositive. Theuseof
IG IG IG
which we use βF(β) in the Pad´e form given by Eq. (4).
Carnahan-Starling hard-sphere free energy, and pertur-
The inverse Laplace transform is calculated by expand-
bationintegralsfromRef. 25givesf(η) c/b= 1.19at
ing exp( βF(β)) in powers of the second summand in ∗ − −
ρ =0.8 and T 0.073.
− IG
Eq. (4) and performing pole integration. The result can ≃
The enumeration function from Eq. (6) is compared
be converted to a closed-form equation:
to the results of MC simulations in Fig. 1a. The simu-
−1 lation points were obtained by patching together central
Ω(e)= b 1+e/c exp[N(f(η) e/b 2c/b)] parts (within 95% of the maximum) of the distribution
(cid:16) p (cid:17) − − (5) functions P(e) exp(N(σ(e) βe)) at different temper-
I1 2(c/b)N 1+e/c , atures. Since P∝(e) defines σ−(e) up to a constant, this
(cid:16) p (cid:17) comparison only serves to show the non-Gaussian shape
where c = a/b, I (x) is the first-order modified Bessel of the curve. This non-Gaussian form of σ(e) forces the
1
function, and e = Eσ3/Nm2 is the reduced internal width of P(e) to decrease with cooling (Fig. 1b). This
energy per particle. Since the argument of the Bessel behaviorisquitedistinctfromthepredictionoftheGaus-
function in Eq. (5) is proportional to the number of sianlandscapemodelinwhichthe widthistemperature-
particles N, its thermodynamic-limit expansion gives independent [15].
the following expression for the enumeration function The temperature dependence of the width is directly
3
0 that the long-range dipolar forces produce essentially a
(a)
-1 mean-field potential for the local translations and rota-
-0.8 e
-2 tions and Eq. (6) for the density of internal energies can
e-1.6 e 0 0.5 1T 1.5 2 be used for the density of inherent structures as well,
0
Ω(e) Ω (e).
-2.4 ≃ φ
e The density of states from simulations approaches the
60 AC
(b) ideal-glassstateevensteeperthanEq.(6)(Fig.1a). This
0.8
trend results in a slightly steeper temperature drop of
40 0.4
0 e(T)fromsimulationscomparedtotheresultofthePad´e
V
c -0.4 approximation(Fig. 2a):
20 0 5 10 15
e(T)= c+c(1+b/T)−2. (8)
0 −
(c)
In particular, instead of leveling off while approaching
the cutoff energy e = 2.07, the system undergoes a
sc2 first-order liquid-liq0uid t−ransition (see below) accompa-
T-1 nied by a discontinuous drop of the internal energy and
IG
asharppeakinthe heatcapacity(Fig.2a,b). Theorien-
0
0 5 10 15 tational structure of the fluid also changes as indicated
1/T
by a step-wise increase in the nematic order parameter
P calculated as the largest eigenvalue of the Q-tensor
FIG. 2: Average energy e(T) (a), heat capacity cV (b), and 2
the configurational entropy sc (c) vs 1/T. The points in (a) Q=(2N)−1 (3ˆe ˆe I), (9)
aretheresultsofMCsimulations: openpointsfortheinternal j j −
energies and close points for the energies of inherent struc- Xj
tures. The solid line in (a) shows the Pad´e approximation
where ˆe is a unit vector along the dipole of particle j.
of Eq. (8). The dashed horizontal lines show the cutoff en- j
ergye0=−2.07, theenergyofnon-interactingdipolarchains The ferroelectric order parameter P1 (normalized total
e = −2.4 (middle line), and the energy of the closed-packed dipole momentofthe liquidM)remainscloseto zeroin-
fcc antiferroelectric crystal eAC = −2.56. The inset in (a) dicatingthattheliquidremainsunpolarized. TheBinder
shows e(T) from simulations, the inset in (b) shows the ne- parameter [28], M = 1 M4 /3 M2 2, also shows
B
− h i h i
maticorderparameterP2(squares)andtheBinderparameter downwardspikespointingtofirst-orderphasetransitions
MB (diamonds). The solid line in (c) is the configurational (inset, Fig. 2b). The internal energy at the first spike of
entropy sc(T) = σ(e(T)) in the Pad´e approximation [Eqs. c drops from the level e to the energy approximately
V 0
(6) and (8)]. The dots show the configurational entropy [Eq.
equaltothe energyofnon-interactingdipolarchains[29]
(6)] with e(T) from simulations. The dashed line in (c) is
e= 2.40(middledashedlineinFig.2a). Thisenergyis
the thermodynamic excess entropy over that of the ideal gas −
slightlyabovetheenergyofclose-packedantiferroelectric
calculated bytemperature integration of simulated cV(T).
fcc crystal e = 2.56 [30].
AC
−
The configurational entropy per particle s (T) =
c
related to the heat capacity since σ(e(T)) from Pad´e thermodynamics [Eqs. (6) and (8)]
deviates from the hyperbolic functions at low tempera-
(δe)2 tures approachingthe ideal glass transition at a positive
P(e) exp , (7)
∝ (cid:20)−2c T2(cid:21) temperature (Fig. 2c). In contrast, when e(T)from sim-
V
ulations is substituted into Eq. (6), the decay of s (T)is
c
where δe is the energy fluctuation and c is the excess faster,reachingtheideal-glassstatenearthepointofthe
V
constant-volumeheatcapacityoverthatofthe idealgas. first-order transition. The DHS fluids thus loses stabil-
The high-temperature portion of the heat capacity, cor- ity when it runs out of configurations when approaching
responding to the distributions shown in Fig. 1b, scales the ideal-glass transition and transforms into a different
linearly with the inverse temperature, c 1/T (Fig. thermodynamic state via a first-order transition.
V
∝
2b). This hyperbolic temperature scaling, often docu- The properties of the low-temperature fluid phase are
mented for structural glass formers at constant pressure peculiar. Both the density ( ρ 2 ) and polarization
k
[27],resultsinalineartemperaturescalingofthesquared ( M 2 )structurefactors∝dohn|ot|shiowcrystallinespa-
k
∝h| | i
width when the distribution P(e) is fitted to a Gaussian tial or orientational order. However, snapshots of simu-
function. This behavior was predicted by models of su- lated configurations (not shown here) indicate the ex-
percooled liquids in terms of configurational excitations istence of chains of dipoles aligned head-to-tail. These
[16]. chains are persistentfor (5 20) 103 MC cycles result-
− ×
The calculation of the energies of inherent structures ing in the overall slow convergence of simulations. The
by conjugate gradient minimization of configurations chainsofparalleldipolesarearrangedinbundleswiththe
along simulated trajectories results in energies φ(T) just locallybody-centeredtetragonal(bct)structureinwhich
slightlybelowe(T)(closedpointsinFig.2a). Thismeans thetwochainsaredisplacedrelativetoeachotherbythe
4
1.6 0.4 transition,a peak at zero transverseprojectionindicates
1.2 0.105 g(r)||0.2 0.095 0.069 tahsemeamlleerrgpeenacke,opfrtehseenctoalutmr/nσar=p1hainset.hAetistohtreospaimcdeiptiomlaer,
)g(r⊥0.8 0 -1 r/0σ 1 fluid, shifts to r/σ = √3/2 characteristic of the closest
0.25
0.167 bctdistance[29]. Theone-dimensionaldipolarchainsex-
0.4
perience strong Landau-Peierls fluctuations with the or-
0 0.071 thogonalmean-squaredisplacementscaling as r⊥2 lT
0 0.5 1 1.5 2 2.5 h i∝
r/σ [33], where l is the average length of the chain. The fits
of the initial portion of g(r⊥) result is almost invariant
FIG. 3: Pair distribution function of transverse interparti- r2 0.07σ2 suggesting that the length of the chains
⊥
cle separations (relativeto thenematic director) at tempera- hincriea≃ses as 1/T.
tures indicated in the plot. The high-temperature results for
The nematic columnar phase transforms into a smec-
T =0.167 (dashed line) and T =0.25 almost coincide on the
tic phase with further coolingas indicatedby the second
scaleoftheplot. Thedistributionofdistancesparallel tothe
peakoftheheatcapacity,adownwardspikeoftheBinder
directorintheinsetshowsthetransitiontothesmecticphase
parameter(Fig. 2b), andthe appearanceofcleardensity
of parallel chains with thelocal bct order.
modulationinthedistributionfunctionofmolecularsep-
arations parallel to the director (inset in Fig. 3).
hard-sphereradius. This arrangementhasthe lowesten- Inconclusion,anew modelofthermodynamicsoflow-
ergy for fluids forming columnar phases of parallel dipo- temperature fluids offers a description of non-Gaussian
lar chains [29, 31]. A body-centered orthorombic lattice landscape in a good agreement with simulations. For
was suggested as the ground state for Stockmayer ferro- the parametersofthe DHS fluids studied here the model
electrics [23]. In the present case, all the chain bundles predicts the excitation profile e(T) to terminate at the
areorientedalongthesamedirector,butthenetmoment ideal-glass transition. The DHS fluid nearly avoids this
is compensated among the oppositely oriented bundles. pointthroughafirst-ordertransitiontoafluidphasewith
The columnar structure is well seen in the pair corre- columnar order.
lation function of transverse separations, which is very This research was supported by the the Air Force
sensitive to columnar order [32] (Fig. 3). At the phase (FA9550-06-C-0084)and NSF (CHE-0616646).
[1] C. A. Angell, K. L. Ngai, G. B. McKenna, and S. W. [18] C.G.GrayandK.E.Gubbins,TheoryofMolecularLiq-
Martin, J. Appl.Phys.88, 3113 (2000). uids (Clarendon Press, Oxford,1984).
[2] I.M.Lifshitz,A.Y.Grosberg,andA.R.Khokhlov,Rev. [19] J. J. Weis and D. Levesque, Phys. Rev. Lett. 71, 2729
Mod. Phys. 50, 683 (1978). (1993).
[3] J. D. Bryngelson and P. G. Wolynes, Proc. Natl. Acad. [20] M. E. van Leeuwen and B. Smit, Phys. Rev. Lett. 71,
Sci. 84, 7524 (1987). 3991 (1993).
[4] K.F. Freed, J. Chem. Phys. 119, 5730 (2003). [21] D. Wei and G. N. Patey, Phys. Rev. Lett. 68, 2043
[5] F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1992).
(1982). [22] D. Wei, G. N. Patey, and A. Perera, Phys. Rev. E 47,
[6] F. H.Stillinger, J. Chem. Phys. 88, 7818 (1988). 506 (1993).
[7] S.Bu¨chnerandA.Heuer,Phys.Rev.E60,6507(1999). [23] G. T. Gao and X. C. Zeng, Phys. Rev. E 61, R2188
[8] S.Sastry, Nature409, 164 (2001). (2000).
[9] J.H.GibbsandE.A.D.Marzio,J.Chem.Phys.28,373 [24] B.GrohandS.Dietrich,Phys.Rev.E63,021203(2001).
(1958). [25] B. Larsen, J. C. Rasaiah, and G. Stell, Mol. Phys. 33,
[10] G.AdamandJ.H.Gibbs,J.Chem.Phys.43,139(1965). 987 (1977).
[11] X.XiaandP.G.Wolynes,Proc.Nat.Acad.Sci.97,2990 [26] A. Saksaengwijit, J. Reinisch, and A. Heuer, Phys. Rev.
(2000). Lett. 93, 235701 (2004).
[12] A.J. Moreno, I.Saika-Voivod,E. Zaccarelli, E. L.Nave, [27] R. Richert and A. C. Angell, J. Chem. Phys. 108, 9016
S.V. Buldyrev,P. Tartaglia, and F. Sciortino, J. Chem. (1998).
Phys.124, 204509 (2006). [28] M.S.S.Challa,D.P.Landau,andK.Binder,Phys.Rev.
[13] M.Mezard, G.Parisi, andM.Virasoro, Spin Glass The- B 34, 1841 (1986).
ory and Beyond (World Scientific,Singapore, 1987). [29] R. Tao and J. M. Sun,Phys. Rev.Lett. 67, 398 (1991).
[14] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, [30] J.M.LuttingerandL.Tisza,Phys.Rev.70,954(1946).
Phys.Rev.A 40, 1045 (1989). [31] A.-P. Hynninen and M. Dijkstra, Phys. Rev. Lett. 94,
[15] B. Derrida, Phys. Rev.B 24, 2613 (1981). 138303 (2005).
[16] D.V.MatyushovandC.A.Angell,J.Chem.Phys.123, [32] D. Wei, Phys. Rev.E 49, 2454 (1994).
034506 (2005); cond-mat/0612627. [33] P. G. de Gennes and P. A. Pinkus, Phys. Kondens. Ma-
[17] P. G. Debenedetti, F. H. Stillinger, and M. S. Shell, J. terie 11, 189 (1970).
Phys.Chem. B 107, 14434 (2003).