Cooling-rate dependence of kinetic and mechanical stability of simulated glasses Hannah Staley,1 Elijah Flenner,2 and Grzegorz Szamel2 1Department of Physics, Colorado State University, Fort Collins, Colorado 80523, USA 2Department of Chemistry, Colorado State University, Fort Collins, Colorado 80523, USA (Dated: January 15, 2015) Recently,ultrastableglasseshavebeencreatedthroughvapordeposition. Subsequently,computer simulation algorithms have been proposed that mimic the vapor deposition process and result in simulated glasses with increased stability. In addition, random pinning has been used to generate very stable glassy configurations without the need for lengthy annealing or special algorithms in- spired by vapor deposition. Kinetic and mechanical stability of experimental ultrastable glasses is compared to those of experimental glasses formed by cooling. We provide the basis for a similar 5 comparison for simulated stable glasses: we analyze the kinetic and mechanical stability of simu- 1 latedglassesformedbycoolingataconstantratebyexaminingthetransformationtimetoaliquid 0 upon rapid re-heating, the inherent structure energies, and the shear modulus. The kinetic and 2 structuralstabilityincreasesslowlywithdecreasingcoolingrate. Themethodsoutlinedherecanbe n used to assess kinetic and mechanical stability of simulated glasses generated by using specialized a algorithms. J 4 PACSnumbers: 64.70.Q-,81.05.Kf 1 ] Vapordepositiononasubstrateheldataround85%of and de Pablo [6] and Lyubimov, Ediger, and de Pablo t the glass transition temperature is used to create glasses [7] modeled the vapor deposition of a binary mixture of f o that have a higher kinetic stability [1] and larger elastic Lennard-Jonesparticlestosimulatethecreationofanul- s moduli [2] than glasses created by cooling from a liq- trastableglass. Theyfoundanenhancedkineticstability . t a uid. Due to these and other desirable properties of these compared to the glass obtained by cooling from a liquid m so-called ultrastable glasses, it is of great interest to un- at a constant rate, but the enhanced stability was less - derstand the differences between glasses created through pronouncedthaninexperiments. Furthermore,thesimu- d vapor deposition and by cooling from a liquid. Since va- latedvapordepositionresultedinglasseswithstructural n por deposited glasses are created one layer at a time, it anisotropy, a higher density, and non-uniform composi- o issuspectedthatenhancedmobilityatthesurfaceallows tion compared to the simulated glasses created by cool- c [ the molecules to find more stable configurations. ing from a liquid. Very recently, Hocky et al. [8] created two dimensional binary Lennard-Jones mixture glasses 1 Discoveryofultrastableglassesinthelabinspiredsev- by randomly pinning particles. Glasses formed by ran- v eralcomputersimulationstudies. Thegoalofthesestud- 9 dom pinning are, by construction, isotropic and uniform ies has been two-fold. First, modeling the vapor depo- 4 (on average). Hocky et al. found that these glasses have sition process can provide insight into the origin of the 4 increasedkineticstabilitybutthepresenceofpinnedpar- 3 stability of vapor deposited glasses. Second, the experi- ticles prevented them from investigating whether these 0 mentalprocedurecanbeusedasaninspirationforthede- glasses also have a larger shear modulus. 1. velopment of computer simulation algorithms that could 0 generate very stable simulated glasses without the need The above mentioned simulational studies focused on 5 for lengthy annealing. L´eonard and Harrowell [3] mod- the properties of stable glasses created by using different 1 eledultrastableglassesusingathree-spinfacilitatedIsing specialized algorithms. However, it is difficult to assess : v model. Theirstudyprovidedevidencethatenhancedsur- the stability of these glasses without knowledge of the i facemobilityduringvapordepositionresultsinimproved stability of simulated glasses generated in a conventional X kinetic stability, but the simplicity of their model pre- way, i.e. by cooling at a constant rate. We emphasize r a ventedthemfromstudyingotherpropertiesofultrastable that simulated glasses generated by cooling are isotropic glasses, such as the shear modulus. On the other end of anduniform. Thus,studyingtheseglassesmayalsoshed the complexity spectrum, Singh et al. [4] modeled vapor somelightontheimportanceoftheanisotropyandcom- depositionoftheorganicmoleculetrehalosebyintroduc- positionalnon-uniformityfortheglassstability. Herewe ing1-5“hot”moleculesaboveasubstrateandthen(arti- createglassesbycoolingfromthesupercooledliquidand ficially)coolingthemtothesubstratetemperature. This weassesstheirkineticandmechanicalstabilitybyexam- procedure resulted in a simulated glass with a higher ki- ining their kinetic stability against reheating, properties neticstability,butalsopronouncedstructuralanisotropy oftheirpotentialenergylandscape,andtheirshearmod- in the direction perpendicular to the substrate. ulus. More recent studies used well-known glass-forming bi- We simulated the Kob-Andersen (KA) [5] 80:20 nary Lennard-Jones mixtures [5]. First, Singh, Ediger, binary Lennard-Jones mixture in three dimensions. 2 The particles interact via the potential V (r) = 1 αβ 4(cid:15) (cid:104)(cid:0)σαβ(cid:1)12−(cid:0)σαβ(cid:1)6(cid:105), and the parameters are (cid:15) = αβ r r BB 0.8 0.5(cid:15) , (cid:15) = 1.5(cid:15) , σ = 0.88σ , and σ = AA AB AA BB AA AB 0.8σ . The majority species are of type A and both AA 0) 0.6 masses are equal. We cut the potential at 2.5σαβ. The =w raapentasudrurtlei(cid:112)tc,sleamasnrAaedtσptAa2riemAns/uee(cid:15)n,mAtrbeAedesrpbineedcierntneigvdsieuttlhycyee.NduW/unVneitits=ssifmwo1r.ui2tl0hlaetnuσegsdAitAnhNg,,(cid:15)ptA=eeAmri8/op0kde0Bri0c-, q(t,ts00..24 2〈δ〉r(t)11110000-0211 10-2 boundary conditions. 100 102 104 t 0 WeranNVTsimulationsusingaNos´e-Hooverthermo- 0 1 2 3 4 10 10 10 10 10 stat. We started with an equilibrated supercooled liquid t at the temperature T = 0.5 (for this system the onset temperature for the slow dynamics is To ≈ 1.0 and the FIG. 1: Average overlap function qs(t,tw = 0) monitored during the heating trajectories (solid lines) compared to the mode-coupling temperature is T = 0.435). We ran 80 c cooling runs at cooling rates T˙ =∆T/∆t of 3.33×10−n averageoverlapfunctionfortheequilibriumliquidatT =0.5 (dashedline). Solidlinescorrespondtoglassesobtainedwith where n = 3, 4, 5, 6, and 7. Each run started at an cooling rates T˙ = 3.33 × 10−n, where n = 3, 4, 5, 6, 7, equilibrium configuration at T = 0.5 and was cooled to and 8 from left to right. The inset shows the mean square T = 0.3. We also ran 4 cooling runs at a cooling rate displacement versus time for the same cooling rates (solid of 3.33×10−8. For T˙ > 3.33×10−8, we started from lines)andfortheequilibriumliquidatT =0.5(dashedline). the final configuration of each of the 80 cooling runs and annealedthesystemforatleast100timesτ determined α from T = 0.5 (we define τ through the average over- itythataparticlemovedoverdistanceaduringthetime α lap function discussed later). For the cooling rate of between t and t+t , where t is the waiting time. w w w 3.33×10−8, we performed 5 different annealing runs for (cid:42) (cid:43) each of the 4 configurations, each beginning with a dif- 1 (cid:88) q (t,t )= q (t,t ) , (1) ferent initial random velocity. In addition, we heated s w N m w m the final configurations of the cooling runs by ramping up the temperature from T = 0.3 to T = 0.5 at a con- where q (t,t ) = Θ(a−|r (t+t )−r (t )|), Θ is m w m w m w stant rate over a time t = 10 (we tried increasing the Heaviside’s step function, and r (t) is the position of m temperature instantaneously from T = 0.3 to T = 0.5 a particle m at a time t. t , the waiting time, is mea- w but found that this resulted in large, non-physical, os- sured from the start of the heating trajectory. We chose cillations of the kinetic and potential energies). After a = 0.25 to be consistent with previous work [16]. We ramping up the temperature to T = 0.5, we continued define τ as when q (τ ,t ) = e−1 for an equilibrium α s α w constant temperature simulations (at T = 0.5) until the system (note that for an equilibrium system q (t,t ) s w mean-square displacement started growing linearly with does not depend on the waiting time, i.e. it is time- time. We refer to this ramping up of temperature and translationally invariant). the subsequent run at T = 0.5 as a heating trajectory. Fig. 1 shows average overlap function q (t,t =0) ob- For T˙ > 3.33×10−8, we ran 80 heating trajectories, in s w tainedwhileheatingglassespreparedatdifferentcooling eachcasestartingfromthefinalconfigurationofadiffer- rates (solid lines) and obtained from the equilibrium run entcoolingrun. ForT˙ =3.33×10−8,werantheheating at T = 0.5 (dashed line). The feature at t = 10 is a trajectories15timesusingdifferentinitialrandomveloc- consequence of changing the thermostat from a ramp-up ities for each of the configurations obtained from the 4 of the temperature at a constant rate from T = 0.3 to cooling runs. Finally, we also performed an equilibrium T = 0.5 to maintaining constant temperature T = 0.5. run at T = 0.5, for at least 100τ . We used HOOMD- α For t > 10, q (t,t = 0) exhibits prolonged plateaus for s w blue [9, 10] for equilibration at T = 0.5 and the cooling the three slowest cooling rates followed by a very rapid simulations. The heating trajectories, the T = 0.3 runs, decayatlatertimes(theglassgeneratedusingtheslowest andthesinglerunatT =0.5usedLAMMPS[11,12]run cooling rate,3.33×10−8, relaxes faster than for ballistic on a GPU [13–15]. motion). With decreasing cooling rate both the plateau We start with an assessment of the kinetic stability of height and its extent increase. This indicates that with the simulated glasses. First, we we examine the time it decreasing cooling rate the cage diameter decreases and takes for the particles to rearrange after the sudden in- ittakeslongerfortheparticlestomoveoutoftheircages. creaseintemperature,i.e. alongtheheatingtrajectories. This is the first indication of increasing kinetic stability We quantify these rearrangements through the average against re-heating. overlap function q (t,t ), which measures the probabil- The inset in Fig. 1 shows the mean square displace- s w 3 2 -7.78 10 -8.22 -7.8 S-8.26 EI -8.3 1 -7.82 S101 0.8 〉U 10-8 10-6 • 10-4 )w0.6 〈 T q(t,ts0.4 -7.84 0.2 0 100 101 102 103 104 -7.86 t 0 1010-8 10-7 10-6 10-5 10-4 10-3 10-2 10-8 10-7 10-6 10-5 10-4 10-3 • • T T FIG.2: DependenceofstabilityratioS =t /τ oncooling FIG. 3: Average potential per particle, (cid:104)U(cid:105), as a function of trans α rate T˙. The inset shows the average overlap function versus cooling rate, T˙. The inset shows average inherent structure time for the cooling rate T˙ =3.33×10−6 and waiting times, energy per particle, EIS, as a function of cooling rate. The from top to bottom, 0, 10, 500, 1000, 3000, 4250, 5000. calculationsweredoneatcoolingrates3.33×10−n,withn= 8,7,6,5,4,3, the lines are to guide the eye. ment (cid:10)δr2(t,t )(cid:11) = N−1(cid:68)(cid:80) [r (t+t )−r (t )]2(cid:69) w n n w n w bility ratios of the most stable simulated glasses, which for a waiting time t =0 obtained while heating glasses w approach approximately 400 [8], we would need to de- prepared at different cooling rates (solid lines) and ob- crease the cooling rate by about 3 decades (to achieve tainedfromtheequilibriumrunatT =0.5(dashedline). We find that (cid:10)δr2(t,t )(cid:11) encodes similar information as stabilityratiosofexperimentalultrastableglasses,which w approach approximately 103.5 [17] we would need to de- q (t,t = 0). In particular, we note that with decreas- s w crease the cooling rate by about 22 decades). We con- ing cooling rate the time required for the mean-square cludethatevensimulatedstableglasseswouldbeimpos- displacement to reach the equilibrium curve increases. sibletoobtainbycooling,atleastwiththecomputational This is another indication of increasing kinetic stability resources available to us at present. against re-heating. To quantify kinetic stability we follow earlier experi- Next, we turn to an assessment of the mechanical sta- mental [17] and simulational [8] studies and we evaluate bility of the simulated glasses. The simplest, albeit in- the transformation time t and the stability ratio S. direct, measure of the mechanical stability is provided trans The transformation time is defined as the time it takes by the inherent structure energy [18], since in the en- for the liquid obtained by heating the glass to return ergy landscape picture of the glass transition, more sta- to equilibrium. In the experimental study of Ref. [17] ble glasses are in a lower energy basin of the landscape the transformation time was obtained by monitoring the and have to overcome a larger energy barrier to deform. time evolution of the dielectric response. We follow the We show the dependence of the average inherent struc- simulationalstudyofRef. [8]andobtainthetransforma- ture energy (obtained using the Fire algorithm [19] in tion time from monitoring the waiting time dependence HOOMD-blue) at T = 0.3 on the cooling rate used to of the average overlap function. Specifically, we calcu- reachthistemperatureintheinsetinFig.3. Weseethat late q (t,t ) for a range of t (see the inset to Fig. 2 for the average inherent structure energy decreases with de- s w w q (t,t ) for various waiting times for T˙ = 3.33×10−6) creasing cooling rate. In Fig. 3 we show the dependence s w andwedefinethewaitingtime-dependentrelaxationtime oftheaveragepotentialenergyperparticleatT =0.3on τ (t ) through the equation q (τ ,t ) = e−1. We define the cooling rate. It also decreases with decreasing cool- s w s s w the transformation time, t , as the waiting time such ing rate. We note that for their vapor deposited glass at trans that τ (t ) = τ , where τ is the equilibrium relaxation T = 0.3, Ediger and de Pablo got an inherent structure s w α α time. Finally, we follow Refs. [8, 17] and we define a energy of -8.35 [7] and a potential energy of -7.8 (esti- stability ratio S =t /τ . mated from [7]). These values are close to our values of trans α In Fig. 2 we show the cooling rate dependence of the -8.30 and -7.86. We note, however, that the vapor de- stability ratio. At larger cooling rates S increases rather position inspired algorithm used in Ref. [7] resulted in quicklywithdecreasingT˙. However,forcoolingratesless glasses that were non-uniform in contrast to our simu- than 3.33×10−5 we get a slower increase of the stability lated glasses, and thus, a literal comparison of the two ratio with cooling rate. For the slowest cooling rates, an studies is impossible. order of magnitude decrease in the cooling rate results A more direct measure of mechanical stability is pro- in a doubling of the stability ratio. Thus, to achieve sta- vided by elastic constants. Ultrastable glasses prepared 4 in experiments [2] and stable glasses obtained in simula- 19 tions using vapor deposition-inspired algorithms [6] have largerelasticconstantsthanglassespreparedbyordinary 18 cooling. Weexamineherethecoolingratedependenceof twheeussheeaanr mapopdruolaucshorfecoeunrtlsyiminutlraotdeducgeldasbsyest.wTooofthuiss[e2n0d]. µ 17 -1q;t)] This method is based on the relationship between corre- 16 S(4 lations of transverse particle displacements in the glass 2q ρ[ 101 and the glass’s shear modulus. The final result is that T 15 kB the shear modulus of a glass can be measured by study- 2 100 q ingthesmallwave-vectorq behavioratlongtimesofthe 14 four-point structure factor -8 -7 -6 -5 -4 10 10 10 10 10 • T (cid:42) (cid:43) 1 (cid:88) S4(q;t)= N g[δrn(t)]g∗[δrm(t)]eiq·[rn(0)−rm(0)] , FIG. 4: Dependence of the shear modulus, µ, on the cool- n,m ing rate T˙ for glasses at T = 0.3. The inset shows (2) 2k Tρ(cid:2)q2S (q;t)(cid:3)−1 versus q plotted at a time t = 100 for B 4 where the weighting function g[δrn(t)] = δr⊥n(t) = T =0.3,forthefourslowestcoolingrates: 3.33×10−n where (r (t)−r (0))⊥ and δr⊥(t) is the component of n=5,6,7,8, from bottom to top. n n n the displacement of particle n perpendicular to the wave-vector q, δr⊥(t) · q = 0. As argued in n Ref. [20], the shear modulus for a glass can be ob- relative increases of the shear modulus achieved in our tained from the small q, long t behavior of S4(q;t), and Ashwin et al. simulational studies are comparable µ = lim lim 2k Tρ(cid:2)q2S (q;t)(cid:3)−1. In practice to the relative differences between the shear modulus of q→0 t→∞ B 4 one can examine the small wave-vector behavior of experimentalultrastableglassesandexperiemntalglasses 2k Tρ(cid:2)q2S (q;t)(cid:3)−1 for times longer than the initial generatedbycooling(19%forindomethacinand15%for B 4 transient dynamics but shorter than time scale for the tris-naphtylbenzene [2]). It is unclear to us whether this aging process of a glass. Also, in order to get good fact has some deeper meaning. statistics one either needs many more independent con- In summary, we showed here that by decreasing the figurationsattheendofeachcoolingprocessthanthe80 cooling rate in simulations, we generated simulated that we generated or one needs to continue running at glasses that have increasing kinetic and mechanical sta- T =0.3 and average over different time origins. In order bility. The highest kinetic stability that we achieved is to measure properties of the glass as prepared using a about 66. Decreasing the cooling rate by one decade, specific cooling rate, one can only average over different which would result in the kinetic stability of about 130, time origins before the start of the aging process. For would require approximately 1 year on a single GPU us- faster cooling rates the aging starts very quickly and the ing HOOMD-blue per one cooling trajectory. This re- range of available time origins is rather small. Fortu- sult establishes a baseline against which specialized al- nately, the present method to calculate the shear mod- gorithms that are inspired by vapor deposition or use ulus does not require very long trajectories [20]. We random pinning should be measured. show 2kBTρ(cid:2)q2S4(q;t)(cid:3)−1 versus q at a time t = 100 We believe that the kinetic stability, quantified by the for T =0.3 in the inset to Fig. 4. We averaged over time stability ratio S =t /τ is the best quantity to com- trans α origins for t≤100 for the 3 faster cooling rates, and for pare stability of glasses generated using different proto- t≤400 for the slowest cooling rate. cols and algorithms. We note that the average potential In Fig. 4 we show the cooling rate dependence of the energy is very sensitive to the average density, compo- shear modulus. We find that the modulus increases with sitional nonuniformities and possibly to anisotropy. In decreasing cooling rate, although the rate of increase fact, our glasses generated by cooling at constant rate seemstobedecreasingwithdecreasingT˙. Attheslowest have lower average potential energy than stable glasses cooling rates an order of magnitude decrease of the cool- generated by Singh et al. [6] using a vapor deposition ing rate results in an increase of the modulus by 4.6%. inspired algorithm. We note the average inherent struc- We note that Ashwin, Bouchbinder, and Procaccia [21] ture energy seems to correlate better with the stability alsoinvestigatedthecoolingratedependenceoftheshear but most likely it is also sensitive to the density and modulus. Theyusedadifferentsystem,intwospatialdi- compositional nonuniformities. Finally, we note that the mensions and thus their results cannot be literally com- relativechangeoftheshearmodulusofsimulatedglasses pared to ours. However, we do note that in a similar formed by cooling at different cooling rates is compara- range of cooling rates our shear modulus increases by ble to the relative difference between the shear modulus 15% and theirs increases by 12%. Interestingly, these ofexperimentalglassescreatedthroughvapordeposition 5 and by cooling. This surprising finding, which deserves [7] I.Lyubimov,M.D.Ediger,andJ.J.dePablo,J.Chem. further investigation, suggests that that the kinetic and Phys. 139, 144505 (2013). mechanical stability are not necessarily correlated. [8] G.M.Hocky,L.Berthier,andD.R.Reichman,J.Chem. Phys. 141, 224503 (2014). The methods presented here can be used to assess ki- [9] http://codeblue.umich.edu/hoomd-blue. netic and mechanical stability of simulated glasses gen- [10] J.A.Anderson,C.D.Lorenz,andA.Travesset,J.Com- erated by using specialized algorithms. put. Phys. 227, 5342 (2008). We gratefully acknowledge the support of NSF grant [11] http://lammps.sandia.gov. CHE 1213401. [12] S. Plimpton, J. Comput. Phys. 117, 1 (1995). [13] The GPU package was created by Mike Brown. [14] W.M.Brown,P.Wang,S.J.Plimpton,andA.N.Thar- rington, Comput. Phys. Commun. 182, 898 (2011). [15] W. M. Brown, A. Kohlmeyer, S. J. Plimpton, and A. N. [1] S.F.Swallen,K.L.Kearns,M.K.Mapes,Y.S.Kim,R. Tharrington,Comput.Phys.Commun.183,449(2012). J.McMahon,M.D.Ediger,T.Wu,L.Yu,andS.Satija, [16] E. Flenner and G. Szamel, J. Phys.-Condens. Mat. 19, Science 315, 353 (2007). 205125 (2007). [2] K.L.Kearns,T.Still,G.Fytas,andM.D.Ediger,Adv. [17] A.Sepu´lveda,M.Tylinski,A.Guiseppi-Elie,R.Richert, Mater. 22, 39 (2010). and M.D. Ediger, Phys. Rev. Lett. 113, 045901 (2014) [3] S.L´eonardandP.Harrowell,J.Chem.Phys.133,244502 [18] F. Sciortino, J. Stat. Mech. P05015 (2005). (2010). [19] E. Bitzek, P. Koskinen, F. Gahler, M. Moseler, and P. [4] S.SinghandJ.J.dePablo,J.Chem.Phys.134,194903 Gumbsch, Phys. Rev. Lett. 97, 170201 (2006). (2011). [20] E.FlennerandG.Szamel,Phys.Rev.Lett.114,025501 [5] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (2015). (1995). [21] AshwinJ.,E.Bouchbinder,andI.Procaccia,Phys.Rev. [6] S. Singh, M. D. Ediger, and J. J. de Pablo, Nat. Mater. E 87, 042310 (2013). 12, 139 (2013).