trunk/chuckDissertation/metallicglass.tex

%!TEX root = /Users/charles/Documents/chuckDissertation/dissertation.tex
\chapter{\label{chap:metalglass}COMPARING MODELS FOR DIFFUSION IN SUPERCOOLED LIQUIDS: THE EUTECTIC COMPOSITION OF THE Ag-Cu ALLOY}

Solid solutions of silver and copper near the eutectic point have been
historical curiosities since Roman emperors used them to cut the
silver content in the ubiquitous {\em denarius} coin by plating copper
blanks with the eutectic mixture.\cite{Pense92} Curiosity in this
alloy continues to the current day since the discovery that it forms a
glassy material when rapidly quenched from the melt.  The eutectic
composition (60.1\% Ag, 39.9\% Cu) was the first glassy metal reported
by Duwez in 1960.\cite{duwez:1136} The Duwez experiments give us an upper
bound on the cooling rate required ($\approx 10^{6} K/s$) in this
alloy, while other metallic glasses (based on Ti-Zr alloys) can be
formed with cooling rates of only $1 K/s$.\cite{Peker93} There is
clearly a great wealth of information waiting to be discovered in the
vitrification and crystallization in liquid alloys.

The Ag-Cu alloy is also of great theoretical interest because it
resembles a model glass-forming Lennard-Jones system which has
correlation functions that decay according to the the famous
Kohlrausch-Williams-Watts ({\sc kww}) law,
\begin{equation}
C(t) \approx  A \exp\left[-(t/\tau)^{\beta}\right],
\label{eq:kww2}
\end{equation}
where $\tau$ is the characteristic relaxation time for the system.

Kob and Andersen observed stretched exponential decay of the van Hove
correlation function in a system comprising an 80-20 mixture of
particles with different well depths $(\epsilon_{AA} \neq
\epsilon_{BB})$ and different range parameters $(\sigma_{AA} \neq
\sigma_{BB})$.\cite{Kob95a,Kob95b} They found that at low
temperatures, the best fitting value for the stretching parameter
$\beta$ was approximately $0.8$.  The same system was later
investigated by Sastry, Debenedetti and Stillinger.\cite{Stillinger98}
They observed nearly identical stretching parameters in the time
dependence of the self intermediate scattering function
$F_{s}(k,t)$.\cite{Hansen86}

In single component Lennard-Jones systems, the stretching parameters
appear to be somewhat lower.  Angelani {\em et al.}\cite{Angelani98} have reported
$\beta \approx 1/2$ for relatively low temperature Lennard-Jones
clusters, and Rabani {\em et al.} have reported similar values for the
decay of correlation functions in defective Lennard-Jones crystals.\cite{Rabani99}

There have been a few recent studies of amorphous metals using fairly
realistic potentials and molecular dynamics methodologies.  Gaukel and
Schober studied diffusion in
$\mbox{Zr}_{67}\mbox{Cu}_{33}$,\cite{Gaukel98} and Qi {\em et al.}
have looked at the static properties ($g(r)$, $s({\bf k})$, etc.) for
alloys of Cu with both Ag and Ni.\cite{Qi99} None of the metallic
liquid simulation studies have observed time correlation functions in
the super-cooled liquid near $T_{g}$.  It would therefore be
interesting to know whether these realistic models have the same
anomalous time dependence as the Lennard-Jones systems, and if so,
whether their stretching behavior is similar to that observed in the
two-component Lennard-Jones systems.

Additionally, bi-metallic alloys present an ideal opportunity
to apply the cage-correlation function methodology which was
developed to study the hopping rate in supercooled
liquids.\cite{Rabani97,Gezelter99,Rabani99,Rabani2000} In particular,
it will be used to test two models for diffusion, both of which use
hopping times which are easily calculated by observing the long-time
decay of the cage correlation function.  The two models are Zwanzig's
model,\cite{Zwanzig83} which is based on the periodic interruption of
harmonic motions around inherent structures, and the Continuous Time
Random Walk ({\sc ctrw})
model,\cite{Blumen83,Klafter94,Klafter96,Shlesinger99} which can be
used to derive transport properties from a random walk on a regular
lattice where the time between jumps is non-uniform.

\section{Theory}

In this section a brief introduction will be given to the models for
diffusion that we will be comparing, as well as a brief description of
how one can use the the cage-correlation function to obtain hopping
rates in liquids.

\subsection{Zwanzig's Model}
In his 1983 paper~\cite{Zwanzig83} on self-diffusion in liquids,
Zwanzig proposed a model for diffusion which consisted of ``cells'' or
basins in which the liquid's configuration oscillates until it
suddenly finds a saddle point on the potential energy surface and
jumps to another basin.  This model was based on and supported by
simulations done by Stillinger and Weber
\cite{Stillinger82,Stillinger83,Weber84,Stillinger85} in which the
liquid configurations generated by molecular dynamics were quenched
periodically by following the steepest descent path to the nearest
local minima on the potential energy surface. Stillinger and Weber
found that as their simulations progressed, the quenched
configurations were stable for short periods of time and then suddenly
jumped (with some re-crossing) to other
configurations.\cite{Stillinger83}

The starting point in Zwanzig's theory is the Green-Kubo
formula~\cite{Berne90,Hansen86} for the self-diffusion constant,
\begin{equation}
D=\frac{1}{3} \int_{0}^{\infty}dt \langle {\bf v}(t) \cdot {\bf v}(0)
\rangle,
\label{eq:kubo}
\end{equation} 
where $\langle {\bf v}(t) \cdot {\bf v}(0) \rangle$ is the
time-dependent velocity autocorrelation function.  Zwanzig's model
predicts the diffusion constant using $\tau$, the lifetime which
characterizes the distribution ($\exp(-t/\tau)$) of residence times in
the cells.

Following a jump, the coherence of the harmonic oscillations is
disrupted, so all correlations between velocities will be destroyed
after each jump.  Zwanzig writes the velocity autocorrelation function
in terms of the velocities of the normal modes in the nearest cell.
The normal mode frequencies are characterized by the normalized
distribution function $\rho(\omega)$, and the time integral can be
solved assuming the time dependence of a damped harmonic oscillator
for each of the normal modes.  In the continuum limit of normal mode
frequencies, one obtains
\begin{equation}
D = \frac{kT}{M}\int_{0}^{\infty} d\omega \rho(\omega)
\frac{\tau}{(1 + \omega^{2}\tau^{2})},
\label{eq:zwanz}
\end{equation} 
where $M$ is the mass of the particles. 

Zwanzig does not explicitly derive the inherent structure normal modes
from the potential energy surface (he used the Debye spectrum for
$\rho(\omega)$.) Moreover, the theory avoids the
problem of how to estimate the lifetime $\tau$ for cell jumps that
destroy the coherent oscillations in the sub-volume.  Nevertheless, the
model fits the experimental results quite well for the self-diffusion
of tetramethylsilane (TMS) and benzene over large ranges in
temperature.\cite{Parkhurst75a,Parkhurst75b}

\subsection{The {\sc ctrw} Model}
In the {\sc ctrw}
model,\cite{Blumen83,Klafter94,Klafter96,Shlesinger99} random walks
take place on a regular lattice but with a distribution of waiting
times,
\begin{equation}
\psi(t) \approx  \left\{ \begin{array}{ll} 
(t/\tau)^{-1-\gamma}, & 0 < \gamma < 1 \\
\frac{1}{\tau} e^{-t/\tau}, & \gamma = 1
\end{array}
\right.,
\label{eq:ctrw}
\end{equation}
between jumps.  $\gamma=1$ represents normal transport, while $\gamma
< 1$ represents transport in a ``fractal time'' regime.  In this
regime, the number of jumps grows in time $t$ as $t^\gamma$ and not
linearly with time.  This dependence implies that the jumps do not
possess a well-defined mean waiting time, and that there is no way to
define a hopping rate for systems that are operating in the fractal
time regime.

Klafter and Zumofen have derived probability distributions for
transport in these systems.\cite{Klafter94} In the $\gamma=1$ limit,
they obtain the standard diffusive behavior in which the diffusion
constant is inversely proportional to $\tau$ and proportional to the
square of the lattice spacing for the random walk ($\sigma_{0}$):
\begin{equation}
\label{eq:CTRW_diff}
D = \frac {\sigma_{0}^{2}} {6\tau}
\end{equation}
This behavior is also suggested by estimates of $\tau$ in
$\mbox{CS}_{2}$ and in Lennard-Jones systems using the cage
correlation function.\cite{Gezelter99,Rabani99}

When $\gamma < 1$, the situation is somewhat more complex.  The
mean-square displacement can be approximated at long times as
\begin{equation}
\label{eq:ctrw_msd}
\langle r^{2}(t) \rangle = 
\frac{4 \sigma_{0}^{2}}{3 \sqrt{\pi}}
\left(\frac{t}{\tau}\right)^{\gamma}
\left(\frac{2-\gamma}{\gamma}\right)^{\gamma}
\frac{1}{\left(2-\gamma\right)^{2}} \Gamma(\frac{7}{2} - \gamma).
\end{equation}
Note that there is no well-defined diffusion constant for transport
that behaves according to this expression.

The waiting time distribution in Eq. (\ref{eq:ctrw}) can be used to
derive a sticking probability,
\begin{equation}
\label{eq:ctrw_phi}
\Phi(t) \approx 1 - \int_{0}^{t} dt' \psi(t'),
\end{equation}
which is the probability of not having made a jump until time t.  The
long-time behavior of this function should be directly related to the
long-time behavior of the cage correlation function, which is
essentially a measurement of the fraction of atoms that are still
within their initial surroundings.  For the $\gamma < 1$ case, the
distribution of waiting times in Eq. (\ref{eq:ctrw}) results in
sticking probability that decays as $t^{-\gamma}$.  This is a very
different behavior than the {\sc kww} law (Eq. (\ref{eq:kww2}) that
has been observed in the defective Lennard-Jones
crystals.\cite{Rabani99}  

Ngai and Liu have asserted that the {\sc kww} decay law
(Eq. (\ref{eq:kww2}) cannot be explained through the use of the
waiting time distribution in Eq. (\ref{eq:ctrw}).\cite{Ngai81}
Instead, they propose a distribution,
\begin{equation}
\psi(t) = c \alpha t^{\alpha -1} e^{-c t^\alpha}
\label{eq:ngai}
\end{equation}
which can explain the experimentally-observed decay.  The form of the
waiting time distribution in Eq. (\ref{eq:ctrw}) is more amenable to
analytical approaches to calculating transport behavior.  Blumen, {\it
et al.} have investigated both proposed functional forms for $\psi(t)$
for recombination phenomena,\cite{Blumen83} and found that
Eq. (\ref{eq:ngai}) leads to time-independent rates at longer times.
They conclude that in order to observe anomalous transport, waiting
time distributions with pathological long-time tails are required.

\subsection{The Cage Correlation Function}

In a recent series of
papers,\cite{Gezelter97,Rabani97,Gezelter98a,Gezelter99,Rabani99,Rabani2000}
Gezelter {\it et al.} have investigated approaches to calculating the hopping
rate ($k_{h} = 1/\tau$) in liquids, supercooled liquids and defective
crystals.  To obtain this estimate, they introduced the {\it cage
correlation function} which measures the rate of change of atomic
surroundings, and associate the long-time decay of this function with
the basin hopping rate for diffusion.  The details on calculating the
cage correlation function can be found in Refs. \cite{Rabani97} and
\cite{Gezelter99}, but the concept will be briefly reviewed here.

An atom's immediate surroundings are best described by the list of
other atoms in the liquid that make up the first solvation shell.
When a diffusive barrier crossing involving the atom has occurred, the
atom has left its immediate surroundings.  Following the barrier
crossing, it will have a slightly different group of atoms surrounding
it.  If one were able to paint identifying numbers on each of the
atoms in a simulation, and kept track of the list of numbers that each
atom could see at any time, then the barrier crossing event would be
evident as a substantial change in this list of neighbors.

The cage correlation function uses a generalized neighbor list to keep
track of each atom's neighbors.  If the list of an atom's neighbors at
time $t$ is identical to the list of neighbors at time $0$, the cage
correlation function has a value of $1$ for that atom.  If any of the
original neighbors are {\em missing} at time $t$, it is assumed that
the atom participated in a hopping event, and the cage correlation
function is $0$.  The atom's surroundings can also change due to
vibrational motion, but at longer times, the cage will reconstitute
itself to include the original members.  Only those events which
result in irreversible changes to the surroundings will cause the cage
to de-correlate at long times.  The mathematical formulation of the
cage correlation function is given in Refs. \cite{Rabani97} and
\cite{Gezelter99}.

Averaging over all atoms in the simulation, and studying the decay of
the cage correlation function gives us a way to measure the hopping
rates directly from relatively short simulations.  Gezelter {\it et al.} have used the
cage correlation function to predict the hopping rates in
atomic~\cite{Rabani97} and molecular~\cite{Gezelter99} liquids, as
well as in defective crystals.\cite{Rabani99,Rabani2000} In the
defective crystals, they found that the cage correlation function, after
being corrected for the initial vibrational behavior at short times,
decayed according to the {\sc kww} law (Eq. (\ref{eq:kww2})) with a
stretching parameter $\beta \approx 1/2$.  Angelani {\em et al.} have
also reported $\beta \approx 1/2$ for relatively low temperature
Lennard-Jones clusters.  This is notably different behavior from the
correlation functions calculated for the Lennard-Jones mixtures that
have been studied by Kob and Andersen,\cite{Kob95a,Kob95b} and Sastry
{\em et al.}\cite{Stillinger98} In this system, the stretching
parameter was closer to 0.8.

This paper will concentrate on the use of the cage correlation
function to obtain hopping times for use by the Zwanzig and {\sc ctrw}
models for diffusion in the $\mbox{Ag}_{6}\mbox{Cu}_{4}$ melt.  We are
also looking for the presence of anomalous dynamical behavior in the
supercooled liquid at temperatures just above the glass transition to
confirm the behavior observed in the simpler Lennard-Jones system.
Section \ref{metglass:sec:details} will outline the computational methods used
to perform the simulations.  Section \ref{metglass:sec:results} contains our
results, and section \ref{metglass:sec:discuss} concludes.

\section{Computational Details}
\label{metglass:sec:details} 

The Sutton-Chen potential with the same
parameterization as Qi {\it et al.}\cite{Qi99} has been chosen for this set of simulations. 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}

In order to study the long-time portions of the correlation functions
in this system without interference from the simulation methodology,
a set of molecular dynamics simulations in the constant-NVE
ensemble was conducted.  The density of the system was taken to be $8.742
\mbox{g}/\mbox{cm}^3$. This density was chosen immediately to the
liquid side of the melting transition from the constant thermodynamic
tension simulations of Qi {\it et al.}.\cite{Qi99} Their simulations gave
excellent estimates of phase and structural behavior, and should be
seen as a starting point for investigations of these materials.

The equations of motion were integrated using the velocity Verlet
integrator with a timestep of 1 fs.  A cutoff radius was used in
which $r^{cut}_{ij}=2\alpha_{ij}$.  Nine independent configurations of
500 atoms were generated by starting from an FCC lattice and randomly
choosing the identities of the particles to match the
$\mbox{Ag}_{6}\mbox{Cu}_{4}$ composition which is near the eutectic
composition for this alloy.\cite{Banhart:1992sv}

Each configuration was started at a temperature of 1350 K (well in
excess of the melting temperature at this density) and then cooled in
50 K increments to 400 K.  At each temperature increment, the systems
re-sampled velocities from a Maxwell-Boltzmann distribution every ps
for 20 ps, then were allowed to equilibrate for 50 ps.  Following the
equilibration period, particle positions and velocities were collected
every ps for 250 ps.  The lower temperature runs (375 K, 350 K, 325 K,
and 300 K) were sampled for 500 ps, 1 ns, 1 ns and 7 ns respectively
to accumulate more accurate long-time statistics. Cooling in this
manner leads to a effective quenching rate of approximately 10$^{11}$
K/s.  This quenching rate is much larger then those that are achieved
through experimental methods which typically are on the order of
10$^{6}$ or 10$^{7}$ K/s depending on the quenching method.\cite{duwez:1136}

\section{Results}
\label{metglass:sec:results}

The simulation results were analyzed for several different structural
and dynamical properties.  Fig \ref{fig:gofr}. shows the pair
correlation function,
\begin{equation}
\label{eq:gofr}
g(r) = \frac {V} {N^2} \left \langle \sum _i  \sum _{j\ne i}\delta(\mathbf
{r-r}_{ij}) \right \rangle,
\end{equation}
at seven temperatures ranging from liquid to supercooled liquid.  The
appearance of a split second peak in the radial distribution function
has been proposed as a signature of the glass transition.\cite{Nagel96} This
behavior is evident in Fig. \ref{fig:gofr} for temperatures below 500
K.

\begin{figure}

        \centering
                \includegraphics[height=3in]{images/gofr_all.pdf}
        \caption{Radial distribution functions for Ag$_6$Cu$_4$. The appearance
                of the split second peak at 500 K indicates the onset of a structural
                change in the super-cooled liquid at temperatures above $T_{g}$.}
        \label{fig:gofr}
\end{figure}
        
Wendt and Abraham have proposed another structural feature to estimate
the location of the glass transition.  Their measure,
\begin{equation}
\label{eq:wendt}
R_{WA} = \frac {g_{min}} {g_{max}},
\end{equation}
is the ratio of the magnitude of the first minimum, $g_{min}$ to
the first maximum, $g_{max}$ in the radial distribution function.
According to their estimates, when the value of $R_{WA}$ reaches 0.14,
the system has passed through the glass transition.\cite{Wendt78} A $T_{g}^{WA}$ of 547 K was observed given a cooling rate of $1.56 \times
10^{11}$ K/s. Goddard, {\it et al.} observed a $T_g^{WA}\approx 500
K$ for a cooling rate of $2\times 10^{12}$ K/s in constant temperature,
constant thermodynamic tension (TtN) simulations\cite{Qi99}.

It should be noted that the split second peak in $g(r)$ appears at $T_g^{WA}$ (Figure \ref{fig:Wendt_abraham}), a temperature for which the diffusion constant is still an appreciable
fraction of it's value in the liquid phase. In addition, the
diffusive behavior continues at the lowest simulated temperature of
300 K indicating that glass transition for our system lies below 300
K. Taking the operational definition of the glass transition to be the
temperature at which the viscosity exceeds $10^{13}$
poise,\cite{Nagel96} it is obvious that neither of these structural
factors should be considered a definitive marker of $T_{g}$. It has
also been noted in previous papers that changes in these structural
factors can occur in certain supercooled metallic liquids
independently of glassy behavior.\cite{Lewis91,Liu92} However, both
structural features seem to mark the onset of anomalous thermodynamic
and dynamical behavior, and should at best be taken as evidence of the
location of the mode coupling theory ({\sc mct}) critical temperature,
$T_{c}$, and not $T_{g}$ itself.
\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/w_a.pdf}
        \caption{Wendt-Abraham parameter ($R_{WA}$) as a
                function of temperature. $T^{WA}_{g} \approx$ 540 K for a cooling rate
                of $1.56 \times 10^{11}$ K/s.}
        \label{fig:Wendt_abraham}
\end{figure}

In order to compute diffusion constants using Zwanzig's model
(Eq. (\ref{eq:zwanz})), one must obtain an estimate of the density of
vibrational states ($\rho(\omega)$) of the liquid.  $\rho(\omega)$ has been obtained in two different ways, first by quenching twenty
high-temperature (1350 K) configurations to the nearest local minimum
on the potential energy surface.  These structures correspond to
inherent structures on the liquid state potential energy surface.
Normal modes were then calculated for each of these inherent
structures by diagonalizing the mass-weighted Hessian to give the
squares of the normal mode frequencies. The second method we used for
determining the vibrational density of states was to compute the
normalized power spectrum of the velocity autocorrelation function
\begin{equation}
\label{eq:vcorr}
\rho_0(\omega) = \int_{-\infty}^\infty \left \langle \mathbf{v}(t) \cdot
\mathbf{v}(0) \right \rangle e^{-i \omega t} dt
\end{equation}
for trajectories which hover just above the inherent structures. These are displayed in Figure \ref{fig:rho}
To calculate the power spectrum, a small amount of kinetic energy (8 K)
was given to each of the twenty inherent structures and the system was
allowed to equilibrate for 30 ps.  After equilibration, the velocity
autocorrelation functions were calculated from relatively short (30
ps) trajectories and the density of states was obtained from a simple
discrete Fourier transform.  At this low temperature, the particles do
not diffuse, so the velocity correlation function doesn't decay due to
hopping.  This implies that the Fourier transform of $\langle
\mathbf{v}(t) \cdot \mathbf{v}(0) \rangle$ is $0$ at $\omega=0$

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/rho.pdf}
        \caption{Density of states calculated from quenched normal modes,
                $\rho_q(\omega)$, and from the Fourier transform of the
                velocity autocorrelation function $\rho_o(\omega)$.}
        \label{fig:rho}
\end{figure}


Both of these methods result in similar gross features for the density
of states, but there are some important differences in their estimates
at low frequencies.  Notably the normal mode density of states
($\rho_{q}(\omega)$) is missing some of the low frequency modes at
frequencies below $10 \mbox{cm}^{-1}$.  The power spectrum
($\rho_{0}(\omega)$) recovers these modes, but gives a much noisier
estimate for the density of states.

It was determined that with a radial cutoff of $2 \alpha_{ij}$ the
potential (and forces) exhibited a large number of discontinuities
which made the minimization into the inherent structures somewhat
challenging.  In this type of potential, the discontinuities at the
cutoff radius {\em cannot} be fixed by shifting the potential as is
done with the Lennard-Jones potential.  This limitation is due to the
non-additive nature of the density portion of the potential function.
In order to address this problem, the minimizations and frequencies
were obtained with a radial cutoff of $3 \alpha_{ij}$, which provided
a surface of adequate smoothness.  The larger cutoff radius requires a
larger system size, and so for the minimizations and the calculation
of the densities of states, we used a total of 1372 atoms.

It should be noted that one possible way to provide a surface without
discontinuities would be to use a shifted form of the density,
\begin{equation}
\rho_{i}=\sum_{j\neq
i}\left(\frac{\alpha_{ij}}{r_{ij}}\right)^{m_{ij}} -
\left(\frac{\alpha_{ij}}{r^{cut}_{ij}}\right)^{m_{ij}}.
\end{equation}
This modification substantially alters the attractive portion of the
potential, and has a deleterious effect on the forces.  A similar
shifting in the Lennard-Jones or other pairwise-additive potential
would not exhibit this problem.  Although this modification does
provide a much smoother potential, it would properly require a
re-parameterization of $c_{i}$ and $D_{ii}$, tasks which are outside
the purview of the current work.

\subsection{Diffusive Transport and Exponential Decay}


Translational diffusion constants were calculated via the Einstein
relation\cite{Berne90}
\begin{equation}
\label{eq:einstein}
D = \lim _{t \rightarrow \infty} \frac {1} {6t} \left \langle \vert
\mathbf {r} _{i} (t) - \mathbf {r} _{i} (0)
\vert ^2 \right \rangle
\end{equation}
using the slope of the long-time portion of the mean square
displacement.  We show an Arrhenius plot of the natural logarithm of
the diffusion constant vs. $1/T$ in Fig. \ref{fig:arrhenius}

The structural features that we noted previously (i.e. the
split-second peak in $g(r)$ and $T_{g}^{WA}$) appear at the same
temperature at which the log of the diffusion constant leaves the
traditional straight-line Arrhenius plot.  


\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/arrhenius.pdf}
        \caption{Arrhenius plot of the self-diffusion constants indicating
                significant deviation from Arrhenius behavior at temperatures below 450 K.}
        \label{fig:arrhenius}
\end{figure}


Truhlar and Kohen have suggested an interpretation of this type of
plot based on Tolman's interpretation of the activation
energy.\cite{Truhlar00,Tolman20,Tolman27} Specifically, positive
convexity in Arrhenius plots requires the average energy of all
diffusing particles rise more slowly then the average energy of all
particles in the ensemble with increasing temperature.  This only
occurs if the microcanonical-ensemble diffusion constant, D(E),
decreases with increasing energy. Truhlar and Kohen explain that one
possible explanation of this behavior is if configuration space can be
decomposed into two types of conformation, one which is ``reactive''
in which diffusive hopping occurs with a rate constant of $K_{R}(E)$,
and one in which there is no diffusion.  If the probability of being
found in the reactive region is given by $P_{R}(E)$, then the overall
diffusion constant,
\begin{equation}
\label{eq:truhlar}
D(E) = P_{R}(E) K_{R}(E),
\end{equation}
would allow $P_{R}(E)$ to decrease with $E$ faster than $K_{R}(E)$ can
increase.  In other words, higher energy systems explore a larger
segment of phase space and spend a smaller fraction of their time in
regions where there is a probability of diffusive hopping. As the
energy is decreased, more time is spent in regions of configuration
space where diffusive hopping is allowed so $D(E)$ should increase
with decreasing $E$.

In Fig. \ref{fig:diffplot} we show the diffusion constants calculated
via the Einstein relation Eq.(\ref{eq:einstein}) along with results
for the Zwanzig and {\sc ctrw} ($\gamma=1$) models using hopping times
obtained from simple linear fits to the long-time decay of the cage
correlation function.  For Zwanzig's model we show results using the
two different estimates for the vibrational density of states.  For
the {\sc ctrw} model, we have used a jump length ($\sigma_0$) of
1.016 {\AA} for {\em all} temperatures shown.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/diffplot.pdf}
        \caption{Self-Diffusion constant for the two models under
                consideration compared to the values computed via standard techniques.}
        \label{fig:diffplot}
\end{figure}


Zwanzig's model is extremely sensitive to the low-$\omega$ portion of
$\rho(\omega)$,\cite{Gezelter99} so it is no surprise that the choice
of the density of states gives a large variation in the predicted
results.  The agreement with the diffusion constants is better at lower
temperatures, but an obviously incorrect temperature
dependence is observed in the higher temperature liquid regime.

The {\sc ctrw} model with $\gamma=1$ and an assumption of fixed jump
distances gives much better agreement in the liquid regime, and the
trend with changing temperature seems to be in excellent agreement
with the Einstein relation.

Delving a bit more deeply into the {\sc ctrw} predictions, it can be
assumed that the distribution of hopping times is well behaved
(i.e. $\gamma=1$) but that the jump {\it distance} is temperature
dependent.  To obtain estimates of the jump distance as a function of
temperature, $\sigma_0(T)$, we can invert (Eq.(\ref{eq:CTRW_diff})) by
multiplying the cage-correlation hopping time by the self-diffusion
constant. The temperature-dependent hopping distances is shown in
Fig. \ref{fig:r0_ctrw}.  Note that these assumptions would lead one to
believe that the average jump distance is increasing sharply as one
approaches the glass transition, which could be an indicator of motion
dominated by Levy flights.\cite{Klafter96,Shlesinger99}

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/sig0_ctrw.pdf}
        \caption{{\sc ctrw} hopping distance, $\sigma_0(T)$, as a function of
                temperature assuming $\gamma=1$ and using the cage correlation hopping
                times.}
        \label{fig:r0_ctrw}
\end{figure}


\subsection{Non-Diffusive Transport     and Non-Exponential Decay}


A much more realistic scenario is that the distribution of hopping
{\it times} changes while the jump distance remains the same at all
temperatures.  In Figs. \ref{fig:exponent} and
\ref{fig:tauhop}, the effects of relaxing the linear fits of 
$\langle r^{2}(t) \rangle$ and of the long time portion of
$ln[C_{cage}(t)]$ are shown.  To fit the mean square displacements, a 
weighted non-linear least-squared fits using the {\sc ctrw} ($\gamma <
1$) expression for the mean square displacement
(Eq. (\ref{eq:ctrw_msd})) was performed.  The only free parameter in the fit is the
jump distance, which is fixed at a value of 1.016 \AA, the optimal
jump distance for the $\gamma=1$ case.  These fits allow estimates of $\gamma$ and the hopping times $\tau$ to be obtained directly from the
mean square displacements.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/exponent.pdf}
        \caption{Comparison of exponential stretching coefficients, $\beta$
                from the cage-correlation function and $\gamma$ from CTRW theory.}
        \label{fig:exponent}
\end{figure}


The cage correlation functions were fit to the {\sc kww} law
(Eq. (\ref{eq:kww2})), after correcting for the short-time (i.e. $< 1$
ps) vibrational decay.  The values of $\gamma$ from the $\langle
r^{2}(t) \rangle$ fits are shown with the {\sc kww} stretching
parameters in Fig. \ref{fig:exponent}.  Note that fits for $\gamma$
and the {\sc kww} stretching parameters both begin to show deviation
from their more normal high-temperature behavior at approximately the
same temperature as the appearance of the structural features
($T_{g}^{WA}$, and the split second peak in $g(r)$) noted above.

There is a small discrepancy between the stretching parameters
($\beta$) for the {\sc kww} fits and the values of $\gamma$ for the
{\sc ctrw} fits to the time dependence of the mean-square
displacement.  Although there is not enough data to make a firm
conclusion, it would appear that there is a small constant offset,
\begin{equation}
\beta \approx \gamma - a
\label{eq:exprel}
\end{equation}
with $a \approx 0.06$ for this system.  This could be due in part to
the sensitivity of the cage correlation function to intermediate
timescales (i.e. 10-30 ps).  Observations at these times will be
susceptible to the effects of short-lived inhomogeneities, while the
{\sc ctrw} fits (done over timescales from 10 ps - 1 ns) will show
only the effects of inhomogeneities which persist for longer times.
Since inhomogeneity can lead to the stretching behavior,\cite{Nagel96}
the values of the stretching parameter obtained from the longer time
fits are the most relevant for understanding the process of
vitrification.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/tauhop.pdf}
        \caption{The characteristic hopping time, $\tau_{hop}$, which
                characterizes the waiting time distribution.  The solid circles
                represent hopping times predicted from {\sc kww} fits to the
                cage-correlation function.  The open triangles represent the
                characteristic time calculated via the {\sc ctrw} model for $\langle
                r^{2}(t) \rangle$.}
        \label{fig:tauhop}
\end{figure}

In Fig. \ref{fig:tauhop}, the hopping times for both
models are shown.  As expected, the hopping times diverge as the temperature is
lowered closer to the glass transition.  There is relatively good
agreement between the hopping times calculated via {\sc ctrw} approach
with those calculated via the cage correlation function, although the
error is as much as a factor of 3 discrepancy at the lowest
temperatures.  

\section{Discussion}
\label{metglass:sec:discuss}

It is relatively clear from using the cage correlation function to
obtain hopping times that the {\sc ctrw} approach to dispersive
transport gives substantially better agreement than Zwanzig's model
which is based on interrupted harmonic motion in the inherent
structures.  The missing piece of the {\sc ctrw} theory is the average
hopping distance $\sigma_{0}$.  A method for calculating $\sigma_{0}$
from short trajectories would make it much more useful for the
calculation of diffusion constants in the $\gamma=1$ limit.  

In the low-temperature supercooled liquid, there are substantial
deviations from the $\gamma=1$ limit of the theory at temperatures
below 500 K.  Below this temperature, the mean square displacement
cannot be fit well with a linear function in time, and the cage
correlation function no longer has a simple exponential decay.  In
this paper, we have attempted to use the waiting time distribution due
to Klafter {\it et
al.}\cite{Blumen83,Klafter94,Klafter96,Shlesinger99} The Laplace
transform of Eq. (\ref{eq:ctrw}) lends itself to the derivation of an
analytical form for the propagator for dispersive
transport.\cite{Klafter94} However, it is somewhat troubling that the
cage correlation function (which has proven itself to be a good
indicator of the average decay of atoms from their initial sites)
cannot be fit with a sticking probability that is commensurate with
Eq. (\ref{eq:ctrw}).

Instead, the long-time decay of the cage
correlation function has been fit with the more familiar Kohlrausch-Williams-Watts
law (Eq. (\ref{eq:kww2})), which appears to be a more accurate model
for the sticking probability.  This point has been addressed by Ngai
and Liu.\cite{Ngai81} It is of some interest that there is a
relatively simple relationship between the stretching parameter in the
fit to the cage correlation function and the value of $\gamma$ from
the {\sc ctrw} model.  One hopes that there is some fundamental
relationship waiting to be discovered between {\sc kww}-like decay and
a {\sc ctrw}-analysis of the propagator for dispersive transport.
Revision:	3496
Committed:	Wed Apr 8 19:13:41 2009 UTC (16 years, 3 months ago) by chuckv
Content type:	application/x-tex
File size:	34811 byte(s)
Log Message:	Final Version