ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/chuckDissertation/metallicglass.tex
Revision: 3483
Committed: Tue Jan 13 14:39:50 2009 UTC (16 years, 6 months ago) by chuckv
Content type: application/x-tex
File size: 34535 byte(s)
Log Message:
Chuck's dissertation for PhD Jan 2009

File Contents

# Content
1 %!TEX root = /Users/charles/Documents/chuckDissertation/dissertation.tex
2 \chapter{\label{chap:metalglass}COMPARING MODELS FOR DIFFUSION IN SUPERCOOLED LIQUIDS: THE EUTECTIC COMPOSITION OF THE AG-CU ALLOY}
3
4 Solid solutions of silver and copper near the eutectic point have been
5 historical curiosities since Roman emperors used them to cut the
6 silver content in the ubiquitous {\em denarius} coin by plating copper
7 blanks with the eutectic mixture.\cite{Pense92} Curiosity in this
8 alloy continues to the current day since the discovery that it forms a
9 glassy material when rapidly quenched from the melt. The eutectic
10 composition (60.1\% Ag, 39.9\% Cu) was the first glassy metal reported
11 by Duwez in 1960.\cite{duwez:1136} The Duwez experiments give us an upper
12 bound on the cooling rate required ($\approx 10^{6} K/s$) in this
13 alloy, while other metallic glasses (based on Ti-Zr alloys) can be
14 formed with cooling rates of only $1 K/s$.\cite{Peker93} There is
15 clearly a great wealth of information waiting to be discovered in the
16 vitrification and crystallization in liquid alloys.
17
18 The Ag-Cu alloy is also of great theoretical interest because it
19 resembles a model glass-forming Lennard-Jones system which has
20 correlation functions that decay according to the the famous
21 Kohlrausch-Williams-Watts ({\sc kww}) law,
22 \begin{equation}
23 C(t) \approx A \exp\left[-(t/\tau)^{\beta}\right].
24 \label{eq:kww2}
25 \end{equation}
26 Kob and Andersen observed stretched exponential decay of the van Hove
27 correlation function in a system comprising an 80-20 mixture of
28 particles with different well depths $(\epsilon_{AA} \neq
29 \epsilon_{BB})$ and different range parameters $(\sigma_{AA} \neq
30 \sigma_{BB})$.\cite{Kob95a,Kob95b} They found that at low
31 temperatures, the best fitting value for the stretching parameter
32 $\beta$ was approximately $0.8$. The same system was later
33 investigated by Sastry, Debenedetti and Stillinger.\cite{Stillinger98}
34 They observed nearly identical stretching parameters in the time
35 dependence of the self intermediate scattering function
36 $F_{s}(k,t)$.\cite{Hansen86}
37
38 In single component Lennard-Jones systems, the stretching parameters
39 appear to be somewhat lower. Angelani {\em et al.} have reported
40 $\beta \approx 1/2$ for relatively low temperature Lennard-Jones
41 clusters, and Rabani {\em et al.} have reported similar values for the
42 decay of correlation functions in defective Lennard-Jones crystals.
43
44 There have been a few recent studies of amorphous metals using fairly
45 realistic potentials and molecular dynamics methodologies. Gaukel and
46 Schober studied diffusion in
47 $\mbox{Zr}_{67}\mbox{Cu}_{33}$,\cite{Gaukel98} and Qi {\em et al.}
48 have looked at the static properties ($g(r)$, $s({\bf k})$, etc.) for
49 alloys of Cu with both Ag and Ni.\cite{Qi99} None of the metallic
50 liquid simulation studies have observed time correlation functions in
51 the super-cooled liquid near $T_{g}$. It would therefore be
52 interesting to know whether these realistic models have the same
53 anomalous time dependence as the Lennard-Jones systems, and if so,
54 whether their stretching behavior is similar to that observed in the
55 two-component Lennard-Jones systems.
56
57 Additionally, bi-metallic alloys present an ideal opportunity for us
58 to apply the cage-correlation function methodology that one of us
59 developed to study the hopping rate in supercooled
60 liquids.\cite{Rabani97,Gezelter99,Rabani99,Rabani2000} In particular,
61 we want to use it to test two models for diffusion, both of which use
62 hopping times which are easily calculated by observing the long-time
63 decay of the cage correlation function. The two models are Zwanzig's
64 model,\cite{Zwanzig83} which is based on the periodic interruption of
65 harmonic motions around inherent structures, and the Continuous Time
66 Random Walk ({\sc ctrw})
67 model,\cite{Blumen83,Klafter94,Klafter96,Shlesinger99} which can be
68 used to derive transport properties from a random walk on a regular
69 lattice where the time between jumps is non-uniform.
70
71 \section{THEORY}
72
73 In this section we give a brief introduction to the models for
74 diffusion that we will be comparing, as well as a brief description of
75 how one can use the the cage-correlation function to obtain hopping
76 rates in liquids.
77
78 \subsection{ZWANZIG'S MODEL}
79 In his 1983 paper~\cite{Zwanzig83} on self-diffusion in liquids,
80 Zwanzig proposed a model for diffusion which consisted of ``cells'' or
81 basins in which the liquid's configuration oscillates until it
82 suddenly finds a saddle point on the potential energy surface and
83 jumps to another basin. This model was based on and supported by
84 simulations done by Stillinger and Weber
85 \cite{Stillinger82,Stillinger83,Weber84,Stillinger85} in which the
86 liquid configurations generated by molecular dynamics were quenched
87 periodically by following the steepest descent path to the nearest
88 local minima on the potential energy surface. Stillinger and Weber
89 found that as their simulations progressed, the quenched
90 configurations were stable for short periods of time and then suddenly
91 jumped (with some re-crossing) to other
92 configurations.\cite{Stillinger83}
93
94 The starting point in Zwanzig's theory is the Green-Kubo
95 formula~\cite{Berne90,Hansen86} for the self-diffusion constant,
96 \begin{equation}
97 D=\frac{1}{3} \int_{0}^{\infty}dt \langle {\bf v}(t) \cdot {\bf v}(0)
98 \rangle,
99 \label{eq:kubo}
100 \end{equation}
101 where $\langle {\bf v}(t) \cdot {\bf v}(0) \rangle$ is the
102 time-dependent velocity autocorrelation function. Zwanzig's model
103 predicts the diffusion constant using $\tau$, the lifetime which
104 characterizes the distribution ($\exp(-t/\tau)$) of residence times in
105 the cells.
106
107 Following a jump, the coherence of the harmonic oscillations is
108 disrupted, so all correlations between velocities will be destroyed
109 after each jump. Zwanzig writes the velocity autocorrelation function
110 in terms of the velocities of the normal modes in the nearest cell.
111 The normal mode frequencies are characterized by the normalized
112 distribution function $\rho(\omega)$, and the time integral can be
113 solved assuming the time dependence of a damped harmonic oscillator
114 for each of the normal modes. In the continuum limit of normal mode
115 frequencies, one obtains
116 \begin{equation}
117 D = \frac{kT}{M}\int_{0}^{\infty} d\omega \rho(\omega)
118 \frac{\tau}{(1 + \omega^{2}\tau^{2})},
119 \label{eq:zwanz}
120 \end{equation}
121 where $M$ is the mass of the particles.
122
123 Zwanzig does not explicitly derive the inherent structure normal modes
124 from the potential energy surface (he used the Debye spectrum for
125 $\rho(\omega)$. Moreover, the theory avoids the
126 problem of how to estimate the lifetime $\tau$ for cell jumps that
127 destroy the coherent oscillations in the sub-volume. Nevertheless, the
128 model fits the experimental results quite well for the self-diffusion
129 of tetramethylsilane (TMS) and benzene over large ranges in
130 temperature.\cite{Parkhurst75a,Parkhurst75b}
131
132 \subsection{THE {\sc ctrw} MODEL}
133 In the {\sc ctrw}
134 model,\cite{Blumen83,Klafter94,Klafter96,Shlesinger99} random walks
135 take place on a regular lattice but with a distribution of waiting
136 times,
137 \begin{equation}
138 \psi(t) \approx \left\{ \begin{array}{ll}
139 (t/\tau)^{-1-\gamma}, & 0 < \gamma < 1 \\
140 \frac{1}{\tau} e^{-t/\tau}, & \gamma = 1
141 \end{array}
142 \right.,
143 \label{eq:ctrw}
144 \end{equation}
145 between jumps. $\gamma=1$ represents normal transport, while $\gamma
146 < 1$ represents transport in a ``fractal time'' regime. In this
147 regime, the number of jumps grows in time $t$ as $t^\gamma$ and not
148 linearly with time. This dependence implies that the jumps do not
149 possess a well-defined mean waiting time, and that there is no way to
150 define a hopping rate for systems that are operating in the fractal
151 time regime.
152
153 Klafter and Zumofen have derived probability distributions for
154 transport in these systems.\cite{Klafter94} In the $\gamma=1$ limit,
155 they obtain the standard diffusive behavior in which the diffusion
156 constant is inversely proportional to $\tau$ and proportional to the
157 square of the lattice spacing for the random walk ($\sigma_{0}$):
158 \begin{equation}
159 \label{eq:CTRW_diff}
160 D = \frac {\sigma_{0}^{2}} {6\tau}
161 \end{equation}
162 This behavior is also suggested by our estimates of $\tau$ in
163 $\mbox{CS}_{2}$ and in Lennard-Jones systems using the cage
164 correlation function.\cite{Gezelter99,Rabani99}
165
166 When $\gamma < 1$, the situation is somewhat more complex. The
167 mean-square displacement can be approximated at long times as
168 \begin{equation}
169 \label{eq:ctrw_msd}
170 \langle r^{2}(t) \rangle =
171 \frac{4 \sigma_{0}^{2}}{3 \sqrt{\pi}}
172 \left(\frac{t}{\tau}\right)^{\gamma}
173 \left(\frac{2-\gamma}{\gamma}\right)^{\gamma}
174 \frac{1}{\left(2-\gamma\right)^{2}} \Gamma(\frac{7}{2} - \gamma).
175 \end{equation}
176 Note that there is no well-defined diffusion constant for transport
177 that behaves according to this expression.
178
179 The waiting time distribution in Eq. (\ref{eq:ctrw}) can be used to
180 derive a sticking probability,
181 \begin{equation}
182 \label{eq:ctrw_phi}
183 \Phi(t) \approx 1 - \int_{0}^{t} dt' \psi(t'),
184 \end{equation}
185 which is the probability of not having made a jump until time t. The
186 long-time behavior of this function should be directly related to the
187 long-time behavior of the cage correlation function, which is
188 essentially a measurement of the fraction of atoms that are still
189 within their initial surroundings. For the $\gamma < 1$ case, the
190 distribution of waiting times in Eq. (\ref{eq:ctrw}) results in
191 sticking probability that decays as $t^{-\gamma}$. This is a very
192 different behavior than the {\sc kww} law (Eq. (\ref{eq:kww2}) that
193 has been observed in the defective Lennard-Jones
194 crystals.\cite{Rabani99}
195
196 Ngai and Liu have asserted that the {\sc kww} decay law
197 (Eq. (\ref{eq:kww2}) cannot be explained through the use of the
198 waiting time distribution in Eq. (\ref{eq:ctrw}).\cite{Ngai81}
199 Instead, they propose a distribution,
200 \begin{equation}
201 \psi(t) = c \alpha t^{\alpha -1} e^{-c t^\alpha}
202 \label{eq:ngai}
203 \end{equation}
204 which can explain the experimentally-observed decay. The form of the
205 waiting time distribution in Eq. (\ref{eq:ctrw}) is more amenable to
206 analytical approaches to calculating transport behavior. Blumen, {\it
207 et al.} have investigated both proposed functional forms for $\psi(t)$
208 for recombination phenomena,\cite{Blumen83} and found that
209 Eq. (\ref{eq:ngai}) leads to time-independent rates at longer times.
210 They conclude that in order to observe anomalous transport, waiting
211 time distributions with pathological long-time tails are required.
212
213 \subsection{THE CAGE CORRELATION FUNCTION}
214
215 In a recent series of
216 papers,\cite{Gezelter97,Rabani97,Gezelter98a,Gezelter99,Rabani99,Rabani2000}
217 one of us has been investigating approaches to calculating the hopping
218 rate ($k_{h} = 1/\tau$) in liquids, supercooled liquids and defective
219 crystals. To obtain this estimate, we introduced the {\it cage
220 correlation function} which measures the rate of change of atomic
221 surroundings, and associate the long-time decay of this function with
222 the basin hopping rate for diffusion. The details on calculating the
223 cage correlation function can be found in Refs. \citen{Rabani97} and
224 \citen{Gezelter99}, but we will briefly review the concept here.
225
226 An atom's immediate surroundings are best described by the list of
227 other atoms in the liquid that make up the first solvation shell.
228 When a diffusive barrier crossing involving the atom has occurred, the
229 atom has left its immediate surroundings. Following the barrier
230 crossing, it will have a slightly different group of atoms surrounding
231 it. If one were able to paint identifying numbers on each of the
232 atoms in a simulation, and kept track of the list of numbers that each
233 atom could see at any time, then the barrier crossing event would be
234 evident as a substantial change in this list of neighbors.
235
236 The cage correlation function uses a generalized neighbor list to keep
237 track of each atom's neighbors. If the list of an atom's neighbors at
238 time $t$ is identical to the list of neighbors at time $0$, the cage
239 correlation function has a value of $1$ for that atom. If any of the
240 original neighbors are {\em missing} at time $t$, it is assumed that
241 the atom participated in a hopping event, and the cage correlation
242 function is $0$. The atom's surroundings can also change due to
243 vibrational motion, but at longer times, the cage will reconstitute
244 itself to include the original members. Only those events which
245 result in irreversible changes to the surroundings will cause the cage
246 to de-correlate at long times. The mathematical formulation of the
247 cage correlation function is given in Refs. \citen{Rabani97} and
248 \citen{Gezelter99}.
249
250 Averaging over all atoms in the simulation, and studying the decay of
251 the cage correlation function gives us a way to measure the hopping
252 rates directly from relatively short simulations. We have used the
253 cage correlation function to predict the hopping rates in
254 atomic~\cite{Rabani97} and molecular~\cite{Gezelter99} liquids, as
255 well as in defective crystals.\cite{Rabani99,Rabani2000} In the
256 defective crystals, we found that the cage correlation function, after
257 being corrected for the initial vibrational behavior at short times,
258 decayed according to the {\sc kww} law (Eq. (\ref{eq:kww2})) with a
259 stretching parameter $\beta \approx 1/2$. Angelani {\em et al.} have
260 also reported $\beta \approx 1/2$ for relatively low temperature
261 Lennard-Jones clusters. This is notably different behavior from the
262 correlation functions calculated for the Lennard-Jones mixtures that
263 have been studied by Kob and Andersen,\cite{Kob95a,Kob95b} and Sastry
264 {\em et al.}\cite{Stillinger98} In this system, the stretching
265 parameter was closer to 0.8.
266
267 This paper will concentrate on the use of the cage correlation
268 function to obtain hopping times for use by the Zwanzig and {\sc ctrw}
269 models for diffusion in the $\mbox{Ag}_{6}\mbox{Cu}_{4}$ melt. We are
270 also looking for the presence of anomalous dynamical behavior in the
271 supercooled liquid at temperatures just above the glass transition to
272 confirm the behavior observed in the simpler Lennard-Jones system.
273 Section \ref{metglass:sec:details} will outline the computational methods used
274 to perform the simulations. Section \ref{metglass:sec:results} contains our
275 results, and section \ref{metglass:sec:discuss} concludes.
276
277 \section{COMPUTATIONAL DETAILS}
278 \label{metglass:sec:details}
279
280 We have chosen the Sutton-Chen potential with the same
281 parameterization as Qi {\it et al.}\cite{Qi99}. Details of this potential can be found in section \ref{introSec:tightbind}. This model was chosen over {\sc EAM} because of better experimental agreement with the heat of solution for Ag and Cu alloys and for better prediction of melting point and liquid state properties. This particular parametarization has been shown to reproduce the experimentally available heat of mixing data for both FCC solid solutions of Ag-Cu and the high-temperature liquid.\cite{sheng:184203} In contrast, the {\sc eam} potential does not reproduce the experimentally observed heat of mixing for the liquid alloy.\cite{MURRAY:1984lr}
282
283 In order to study the long-time portions of the correlation functions
284 in this system without interference from the simulation methodology,
285 we carried out molecular dynamics simulations in the constant-NVE
286 ensemble. The density of the system was taken to be $8.742
287 \mbox{g}/\mbox{cm}^3$. This density was chosen immediately to the
288 liquid side of the melting transition from the constant thermodynamic
289 tension simulations of Qi {\it et al.}. Their simulations gave
290 excellent estimates of phase and structural behavior, and should be
291 seen as a starting point for investigations of these materials.
292
293 The equations of motion were integrated using the velocity Verlet
294 integrator with a timestep of 1 fs. A cutoff radius was used in
295 which $r^{cut}_{ij}=2\alpha_{ij}$. Nine independent configurations of
296 500 atoms were generated by starting from an FCC lattice and randomly
297 choosing the identities of the particles to match the
298 $\mbox{Ag}_{6}\mbox{Cu}_{4}$ composition which is near the eutectic
299 composition for this alloy.\cite{Banhart:1992sv}
300
301 Each configuration was started at a temperature of 1350 K (well in
302 excess of the melting temperature at this density) and then cooled in
303 50 K increments to 400 K. At each temperature increment, the systems
304 re-sampled velocities from a Maxwell-Boltzmann distribution every ps
305 for 20 ps, then were allowed to equilibrate for 50 ps. Following the
306 equilibration period, we collected particle positions and velocities
307 every ps for 250 ps. The lower temperature runs (375 K, 350 K, 325 K,
308 and 300 K) were sampled for 500 ps, 1 ns, 1 ns and 7 ns respectively
309 to accumulate more accurate long-time statistics. Cooling in this
310 manner leads to a effective quenching rate of approximately 10$^{11}$
311 K/s. This quenching rate is much larger then those that are achieved
312 through experimental methods which typically are on the order of
313 10$^{6}$ or 10$^{7}$ K/s depending on the quenching method.\cite{duwez:1136}
314
315 \section{RESULTS}
316 \label{metglass:sec:results}
317
318 The simulation results were analyzed for several different structural
319 and dynamical properties. Fig \ref{fig:gofr}. shows the pair
320 correlation function,
321 \begin{equation}
322 \label{eq:gofr}
323 g(r) = \frac {V} {N^2} \left \langle \sum _i \sum _{j\ne i}\delta(\mathbf
324 {r-r}_{ij}) \right \rangle,
325 \end{equation}
326 at seven temperatures ranging from liquid to supercooled liquid. The
327 appearance of a split second peak in the radial distribution function
328 has been proposed as a signature of the glass transition.\cite{Nagel96} This
329 behavior is evident in Fig. \ref{fig:gofr} for temperatures below 500
330 K.
331
332 \begin{figure}
333
334 \centering
335 \includegraphics[height=3in]{images/gofr_all.pdf}
336 \caption{Radial distribution functions for Ag$_6$Cu$_4$. The appearance
337 of the split second peak at 500 K indicates the onset of a structural
338 change in the super-cooled liquid at temperatures above $T_{g}$.}
339 \label{fig:gofr}
340 \end{figure}
341
342 Wendt and Abraham have proposed another structural feature to estimate
343 the location of the glass transition. Their measure,
344 \begin{equation}
345 \label{eq:wendt}
346 R_{WA} = \frac {g_{min}} {g_{max}},
347 \end{equation}
348 is the ratio of the magnitude of the first minimum, $g_{min}$ to
349 the first maximum, $g_{max}$ in the radial distribution function.
350 According to their estimates, when the value of $R_{WA}$ reaches 0.14,
351 the system has passed through the glass transition.\cite{Wendt78} We
352 observed a $T_{g}^{WA}$ of 547 K given a cooling rate of $1.56 \times
353 10^{11}$ K/s. Goddard, {\it et al.} observed a $T_g^{WA}\approx 500
354 K$ for a cooling rate of $2\times 10^12$ K/s in constant temperature,
355 constant thermodynamic tension (TtN) simulations\cite{Qi99}.
356
357 We note that the split second peak in $g(r)$ appears at $T_g^{WA}$, a
358 temperature for which the diffusion constant is still an appreciable
359 fraction of it's value in the liquid phase. We also note that
360 diffusive behavior continues at the lowest simulated temperature of
361 300 K indicating that glass transition for our system lies below 300
362 K. Taking the operational definition of the glass transition to be the
363 temperature at which the viscosity exceeds $10^{13}$
364 poise,\cite{Nagel96} it is obvious that neither of these structural
365 factors should be considered a definitive marker of $T_{g}$. It has
366 also been noted in previous papers that changes in these structural
367 factors can occur in certain supercooled metallic liquids
368 independently of glassy behavior.\cite{Lewis91,Liu92} However, both
369 structural features seem to mark the onset of anomalous thermodynamic
370 and dynamical behavior, and should at best be taken as evidence of the
371 location of the mode coupling theory ({\sc mct}) critical temperature,
372 $T_{c}$, and not $T_{g}$ itself.
373 \begin{figure}[htbp]
374 \centering
375 \includegraphics[height=3in]{images/w_a.pdf}
376 \caption{Wendt-Abraham parameter ($R_{WA}$) as a
377 function of temperature. $T^{WA}_{g} \approx$ 540 K for a cooling rate
378 of $1.56 \times 10^{11}$ K/s.}
379 \label{fig:Wendt_abraham}
380 \end{figure}
381
382 In order to compute diffusion constants using Zwanzig's model
383 (Eq. (\ref{eq:zwanz}), one must obtain an estimate of the density of
384 vibrational states ($\rho(\omega)$) of the liquid. We obtained the
385 $\rho(\omega)$ in two different ways, first by quenching twenty
386 high-temperature (1350 K) configurations to the nearest local minimum
387 on the potential energy surface. These structures correspond to
388 inherent structures on the liquid state potential energy surface.
389 Normal modes were then calculated for each of these inherent
390 structures by diagonalizing the mass-weighted Hessian to give the
391 squares of the normal mode frequencies. The second method we used for
392 determining the vibrational density of states was to compute the
393 normalized power spectrum of the velocity autocorrelation function
394 \begin{equation}
395 \label{eq:vcorr}
396 \rho_0(\omega) = \int_{-\infty}^\infty \left \langle \mathbf{v}(t) \cdot
397 \mathbf{v}(0) \right \rangle e^{-i \omega t} dt
398 \end{equation}
399 for trajectories which hover just above the inherent structures. To
400 calculate the power spectrum, a small amount of kinetic energy (8 K)
401 was given to each of the twenty inherent structures and the system was
402 allowed to equilibrate for 30 ps. After equilibration, the velocity
403 autocorrelation functions were calculated from relatively short (30
404 ps) trajectories and the density of states was obtained from a simple
405 discrete Fourier transform. At this low temperature, the particles do
406 not diffuse, so the velocity correlation function doesn't decay due to
407 hopping. This implies that the Fourier transform of $\langle
408 \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle$ is $0$ at $\omega=0$
409
410 \begin{figure}[htbp]
411 \centering
412 \includegraphics[height=3in]{images/rho.pdf}
413 \caption{Density of states calculated from quenched normal modes,
414 $\rho_q(\omega)$, and from the Fourier transform of the
415 velocity autocorrelation function $\rho_o(\omega)$.}
416 \label{fig:rho}
417 \end{figure}
418
419
420
421
422 Both of these methods result in similar gross features for the density
423 of states, but there are some important differences in their estimates
424 at low frequencies. Notably the normal mode density of states
425 ($\rho_{q}(\omega)$) is missing some of the low frequency modes at
426 frequencies below $10 \mbox{cm}^{-1}$. The power spectrum
427 ($\rho_{0}(\omega)$) recovers these modes, but gives a much noisier
428 estimate for the density of states.
429
430 We determined that with a radial cutoff of $2 \alpha_{ij}$ the
431 potential (and forces) exhibited a large number of discontinuities
432 which made the minimization into the inherent structures somewhat
433 challenging. In this type of potential, the discontinuities at the
434 cutoff radius {\em cannot} be fixed by shifting the potential as is
435 done with the Lennard-Jones potential. This limitation is due to the
436 non-additive nature of the density portion of the potential function.
437 In order to address this problem, the minimizations and frequencies
438 were obtained with a radial cutoff of $3 \alpha_{ij}$, which provided
439 a surface of adequate smoothness. The larger cutoff radius requires a
440 larger system size, and so for the minimizations and the calculation
441 of the densities of states, we used a total of 1372 atoms.
442
443 We note that one possible way to provide a surface without
444 discontinuities would be to use a shifted form of the density,
445 \begin{equation}
446 \rho_{i}=\sum_{j\neq
447 i}\left(\frac{\alpha_{ij}}{r_{ij}}\right)^{m_{ij}} -
448 \left(\frac{\alpha_{ij}}{r^{cut}_{ij}}\right)^{m_{ij}}.
449 \end{equation}
450 This modification substantially alters the attractive portion of the
451 potential, and has a deleterious effect on the forces. A similar
452 shifting in the Lennard-Jones or other pairwise-additive potential
453 would not exhibit this problem. Although this modification does
454 provide a much smoother potential, it would properly require a
455 re-parameterization of $c_{i}$ and $D_{ii}$, tasks which are outside
456 the purview of the current work.
457
458 \subsection{DIFFUSIVE TRANSPORT AND EXPONENTIAL DECAY}
459
460
461 Translational diffusion constants were calculated via the Einstein
462 relation\cite{Berne90}
463 \begin{equation}
464 \label{eq:einstein}
465 D = \lim _{t \rightarrow \infty} \frac {1} {6t} \left \langle \vert
466 \mathbf {r} _{i} (t) - \mathbf {r} _{i} (0)
467 \vert ^2 \right \rangle
468 \end{equation}
469 using the slope of the long-time portion of the mean square
470 displacement. We show an Arrhenius plot of the natural logarithm of
471 the diffusion constant vs. $1/T$ in Fig. \ref{fig:arrhenius}
472
473 The structural features that we noted previously (i.e. the
474 split-second peak in $g(r)$ and $T_{g}^{WA}$) appear at the same
475 temperature at which the log of the diffusion constant leaves the
476 traditional straight-line Arrhenius plot.
477
478
479 \begin{figure}[htbp]
480 \centering
481 \includegraphics[height=3in]{images/arrhenius.pdf}
482 \caption{Arrhenius plot of the self-diffusion constants indicating
483 significant deviation from Arrhenius behavior at temperatures below 450 K.}
484 \label{fig:arrhenius}
485 \end{figure}
486
487
488
489 Truhlar and Kohen have suggested an interpretation of this type of
490 plot based on Tolman's interpretation of the activation
491 energy.\cite{Truhlar00,Tolman20,Tolman27} Specifically, positive
492 convexity in Arrhenius plots requires the average energy of all
493 diffusing particles rise more slowly then the average energy of all
494 particles in the ensemble with increasing temperature. This only
495 occurs if the microcanonical-ensemble diffusion constant, D(E),
496 decreases with increasing energy. Truhlar and Kohen explain that one
497 possible explanation of this behavior is if configuration space can be
498 decomposed into two types of conformation, one which is ``reactive''
499 in which diffusive hopping occurs with a rate constant of $K_{R}(E)$,
500 and one in which there is no diffusion. If the probability of being
501 found in the reactive region is given by $P_{R}(E)$, then the overall
502 diffusion constant,
503 \begin{equation}
504 \label{eq:truhlar}
505 D(E) = P_{R}(E) K_{R}(E),
506 \end{equation}
507 would allow $P_{R}(E)$ to decrease with $E$ faster than $K_{R}(E)$ can
508 increase. In other words, higher energy systems explore a larger
509 segment of phase space and spend a smaller fraction of their time in
510 regions where there is a probability of diffusive hopping. As the
511 energy is decreased, more time is spent in regions of configuration
512 space where diffusive hopping is allowed so $D(E)$ should increase
513 with decreasing $E$.
514
515 In Fig. \ref{fig:diffplot} we show the diffusion constants calculated
516 via the Einstein relation Eq.(\ref{eq:einstein}) along with results
517 for the Zwanzig and {\sc ctrw} ($\gamma=1$) models using hopping times
518 obtained from simple linear fits to the long-time decay of the cage
519 correlation function. For Zwanzig's model we show results using the
520 two different estimates for the vibrational density of states. For
521 the {\sc ctrw} model, we have used a jump length ($\sigma_0$) of
522 1.016 {\AA} for {\em all} temperatures shown.
523
524 \begin{figure}[htbp]
525 \centering
526 \includegraphics[height=3in]{images/diffplot.pdf}
527 \caption{Self-Diffusion constant for the two models under
528 consideration compared to the values computed via standard techniques.}
529 \label{fig:diffplot}
530 \end{figure}
531
532
533 Zwanzig's model is extremely sensitive to the low-$\omega$ portion of
534 $\rho(\omega)$,\cite{Gezelter99} so it is no surprise that the choice
535 of the density of states gives a large variation in the predicted
536 results. The agreement with the diffusion constants is better at lower
537 temperatures, but we observe an obviously incorrect temperature
538 dependence in the higher temperature liquid regime.
539
540 The {\sc ctrw} model with $\gamma=1$ and an assumption of fixed jump
541 distances gives much better agreement in the liquid regime, and the
542 trend with changing temperature seems to be in excellent agreement
543 with the Einstein relation.
544
545 Delving a bit more deeply into the {\sc ctrw} predictions, we can
546 assume that the distribution of hopping times is well behaved
547 (i.e. $\gamma=1$) but that the jump {\it distance} is temperature
548 dependent. To obtain estimates of the jump distance as a function of
549 temperature, $\sigma_0(T)$, we can invert (Eq.(\ref{eq:CTRW_diff})) by
550 multiplying the cage-correlation hopping time by the self-diffusion
551 constant. We show the temperature-dependent hopping distances in
552 Fig. \ref{fig:r0_ctrw}. Note that these assumptions would lead us to
553 believe that the average jump distance is increasing sharply as one
554 approaches the glass transition, which could be an indicator of motion
555 dominated by Levy flights.\cite{Klafter96,Shlesinger99}
556
557 \begin{figure}[htbp]
558 \centering
559 \includegraphics[height=3in]{images/sig0_ctrw.pdf}
560 \caption{{\sc ctrw} hopping distance, $\sigma_0(T)$, as a function of
561 temperature assuming $\gamma=1$ and using the cage correlation hopping
562 times.}
563 \label{fig:r0_ctrw}
564 \end{figure}
565
566
567 \subsection{NON-DIFFUSIVE TRNASPORT AND NON-EXPONENTIAL DECAY}
568
569
570 A much more realistic scenario is that the distribution of hopping
571 {\it times} changes while the jump distance remains the same at all
572 temperatures. In Figs. \ref{fig:exponent} and
573 \ref{fig:tauhop}, we show the effects of relaxing the linear fits of
574 $\langle r^{2}(t) \rangle$ and of the long time portion of
575 $ln[C_{cage}(t)]$. To fit the mean square displacements, we performed
576 weighted non-linear least-squared fits using the {\sc ctrw} ($\gamma <
577 1$) expression for the mean square displacement
578 (Eq. (\ref{eq:ctrw_msd})). The only free parameter in the fit is the
579 jump distance, which we fixed at a value of 1.016 \AA, the optimal
580 jump distance for the $\gamma=1$ case. These fits allow us to obtain
581 estimates of $\gamma$ and the hopping times $\tau$ directly from the
582 mean square displacements.
583
584 \begin{figure}[htbp]
585 \centering
586 \includegraphics[height=3in]{images/exponent.pdf}
587 \caption{Comparison of exponential stretching coefficients, $\beta$
588 from the cage-correlation function and $\gamma$ from CTRW theory.}
589 \label{fig:exponent}
590 \end{figure}
591
592
593
594 The cage correlation functions were fit to the {\sc kww} law
595 (Eq. (\ref{eq:kww2})), after correcting for the short-time (i.e. $< 1$
596 ps) vibrational decay. The values of $\gamma$ from the $\langle
597 r^{2}(t) \rangle$ fits are shown with the {\sc kww} stretching
598 parameters in Fig. \ref{fig:exponent}. Note that fits for $\gamma$
599 and the {\sc kww} stretching parameters both begin to show deviation
600 from their more normal high-temperature behavior at approximately the
601 same temperature as the appearance of the structural features
602 ($T_{g}^{WA}$, and the split second peak in $g(r)$) noted above.
603
604 There is a small discrepancy between the stretching parameters
605 ($\beta$) for the {\sc kww} fits and the values of $\gamma$ for the
606 {\sc ctrw} fits to the time dependence of the mean-square
607 displacement. Although there is not enough data to make a firm
608 conclusion, it would appear that there is a small constant offset,
609 \begin{equation}
610 \beta \approx \gamma - a
611 \label{eq:exprel}
612 \end{equation}
613 with $a \approx 0.06$ for this system. This could be due in part to
614 the sensitivity of the cage correlation function to intermediate
615 timescales (i.e. 10-30 ps). Observations at these times will be
616 susceptible to the effects of short-lived inhomogeneities, while the
617 {\sc ctrw} fits (done over timescales from 10 ps - 1 ns) will show
618 only the effects of inhomogeneities which persist for longer times.
619 Since inhomogeneity can lead to the stretching behavior,\cite{Nagel96}
620 the values of the stretching parameter obtained from the longer time
621 fits are the most relevant for understanding the process of
622 vitrification.
623
624 \begin{figure}[htbp]
625 \centering
626 \includegraphics[height=3in]{images/tauhop.pdf}
627 \caption{The characteristic hopping time, $\tau_{hop}$, which
628 characterizes the waiting time distribution. The solid circles
629 represent hopping times predicted from {\sc kww} fits to the
630 cage-correlation function. The open triangles represent the
631 characteristic time calculated via the {\sc ctrw} model for $\langle
632 r^{2}(t) \rangle$.}
633 \label{fig:tauhop}
634 \end{figure}
635
636 In Fig. \ref{fig:tauhop}, we show the hopping times for both
637 models. As expected, the hopping times diverge as the temperature is
638 lowered closer to the glass transition. There is relatively good
639 agreement between the hopping times calculated via {\sc ctrw} approach
640 with those calculated via the cage correlation function, although the
641 error is as much as a factor of 3 discrepancy at the lowest
642 temperatures.
643
644 \section{DISCUSSION}
645 \label{metglass:sec:discuss}
646
647 It is relatively clear from using the cage correlation function to
648 obtain hopping times that the {\sc ctrw} approach to dispersive
649 transport gives substantially better agreement than Zwanzig's model
650 which is based on interrupted harmonic motion in the inherent
651 structures. The missing piece of the {\sc ctrw} theory is the average
652 hopping distance $\sigma_{0}$. A method for calculating $\sigma_{0}$
653 from short trajectories would make it much more useful for the
654 calculation of diffusion constants in the $\gamma=1$ limit.
655
656 In the low-temperature supercooled liquid, there are substantial
657 deviations from the $\gamma=1$ limit of the theory at temperatures
658 below 500 K. Below this temperature, the mean square displacement
659 cannot be fit well with a linear function in time, and the cage
660 correlation function no longer has a simple exponential decay. In
661 this paper, we have attempted to use the waiting time distribution due
662 to Klafter {\it et
663 al.}\cite{Blumen83,Klafter94,Klafter96,Shlesinger99} The Laplace
664 transform of Eq. (\ref{eq:ctrw}) lends itself to the derivation of an
665 analytical form for the propagator for dispersive
666 transport.\cite{Klafter94} However, it is somewhat troubling that the
667 cage correlation function (which has proven itself to be a good
668 indicator of the average decay of atoms from their initial sites)
669 cannot be fit with a sticking probability that is commensurate with
670 Eq. (\ref{eq:ctrw}).
671
672 Instead, we have been able to fit the long-time decay of the cage
673 correlation function with the more familiar Kohlrausch-Williams-Watts
674 law (Eq. (\ref{eq:kww2})), which appears to be a more accurate model
675 for the sticking probability. This point has been addressed by Ngai
676 and Liu.\cite{Ngai81} It is of some interest that there is a
677 relatively simple relationship between the stretching parameter in the
678 fit to the cage correlation function and the value of $\gamma$ from
679 the {\sc ctrw} model. One hopes that there is some fundamental
680 relationship waiting to be discovered between {\sc kww}-like decay and
681 a {\sc ctrw}-analysis of the propagator for dispersive transport.