Table Of Content

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 ferroelectric 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étruncation 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é 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é 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é 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ê ê 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é approximation where ê 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é 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é 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).