trunk/chuckDissertation/bulkmod.tex

%!TEX root = /Users/charles/Documents/chuckDissertation/dissertation.tex
\chapter{\label{chap:bulkmod}BREATHING MODE DYNAMICS AND ELASTIC PROPERTIES OF GOLD NANOPARTICLES}

In metallic nanoparticles, the relatively large surface area to volume ratio
induces a number of well-known size-dependent phenomena.  Notable
among these are the depression of the bulk melting
temperature,\cite{Buffat:1976yq,el-sayed00,el-sayed01} surface melting
transitions, increased room-temperature alloying
rates,\cite{ShibataT._ja026764r} changes in the breathing mode
frequencies,\cite{delfatti99,henglein99,hartland02a,hartland02c} and
rapid heat transfer to the surrounding
medium.\cite{HuM._jp020581+,hartland02d}

This paper reports on atomic-level simulations of the transient
response of metallic nanoparticles to the nearly instantaneous heating
undergone when photons are absorbed during ultrafast laser excitation
experiments.  The time scale for heating (determined by the {\it e-ph}
coupling constant) is faster than a single period of the breathing
mode for spherical nanoparticles.\cite{Simon2001,HartlandG.V._jp0276092}
Hot-electron pressure and direct lattice heating contribute to the
thermal excitation of the atomic degrees of freedom, and both
mechanisms are rapid enough to coherently excite the breathing mode of
the spherical particles.\cite{Hartland00}

There are questions posed by the experiments that may be easiest to
answer via computer simulation.  For example, the dephasing seen
following coherent excitation of the breathing mode may be due to
inhomogeneous size distributions in the sample, but it may also be due
to softening of the breathing mode vibrational frequency following a
melting transition.   Additionally, there are properties (such as the
bulk modulus) that may be nearly impossible to obtain experimentally,
but which are relatively easily obtained via simulation techniques.

We outline our simulation techniques in section \ref{bulkmod:sec:details}.
Results are presented in section \ref{bulkmod:sec:results}.  We discuss our
results in terms of Lamb's classical theory of elastic spheres in
section \ref{bulkmod:sec:discussion}.

\section{COMPUTATIONAL DETAILS}
\label{bulkmod:sec:details}

Spherical Au nanoparticles were created in a standard FCC lattice at
four different radii [20{\AA} (1926 atoms), 25{\AA} (3884 atoms),
30{\AA} (6602 atoms), and 35{\AA} (10606 atoms)].  To create spherical
nanoparticles, a large FCC lattice was built at the normal Au lattice
spacing (4.08 \AA) and any atoms outside the target radius were
excluded from the simulation.

\subsection{SIMULATION METHODOLOGY}

Potentials were calculated using the Voter-Chen
parameterization~\cite{Voter:87} of the Embedded Atom Method ({\sc
eam}), which has been widely used for MD simulations of metallic
particles.%\cite{Voter:87,Daw84,FBD86,johnson89,Lu97} 
Like other transition metal potentials,%\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02}
{\sc eam} describes bonding in metallic systems by including an
attractive interaction which models the stabilization of a positively
charged metal core ion in a sea of surrounding valence electrons.  A
repulsive pair potential describes the interactions of the core ions
with each other.  The {\sc eam} potential has the form:
\begin{eqnarray}
V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
\phi_{ij}({\bf r}_{ij})  \\
\rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij})
\end{eqnarray}
where $\phi_{ij}$ is the (primarily repulsive) pairwise interaction
between atoms $i$ and $j$.  $F_{i}[\rho_{i}]$ is an embedding function
that determines the energy to embed the positively-charged core, $i$,
in the electron density, $\rho_{i}$, due to the valence electrons of
the surrounding $j$ atoms, and $f_{j}(r)$ describes the radial
dependence of the field due to atom $j$. There is a cutoff distance,
$r_{cut}$, which limits the summations in the {\sc eam} equation to
the few dozen atoms surrounding atom $i$.  In these simulations, a
cutoff radius of 10~\AA\ was used.

%Mixing rules as outlined by Johnson~\cite{johnson89} were used to
%compute the heterogenous pair potential,
%
%\begin{eqnarray}
%\label{eq:johnson}
%\phi_{ab}(r)=\frac{1}{2}\left(
%\frac{f_{b}(r)}{f_{a}(r)}\phi_{aa}(r)+
%\frac{f_{a}(r)}{f_{b}(r)}\phi_{bb}(r)
%\right).
%\end{eqnarray}

Before starting the molecular dynamics runs, a relatively short
steepest-descent minimization was performed to relax the lattice in
the initial configuration.  To facilitate the study of the larger
particles, simulations were run in parallel over 16 processors using
Plimpton's force-decomposition method.\cite{plimpton93}

To mimic the events following the absorption of light in the ultrafast
laser heating experiments, we have used a simple two-step process to
prepare the simulations.  Instantaneous heating of the lattice was
performed by sampling atomic velocities from a Maxwell-Boltzmann
distribution set to twice the target temperature for the simulation.
By equipartition, approximately half of the initial kinetic energy of
the system winds up in the potential energy of the system.  The system
was a allowed a very short (10 fs) evolution period with the new
velocities.  

Following this excitation step, the particles evolved under
Nos\'{e}-Hoover NVT dynamics~\cite{hoover85} for 40 ps.  Given the
mass of the constituent metal atoms, time steps of 5 fs give excellent
energy conservation in standard NVE integrators, so the same time step
was used in the NVT simulations.  Target temperatures for these
particles spanned the range from 300 K to 1350 K in 50 K intervals.
Five independent samples were run for each particle and temperature.
The results presented below are averaged properties for each of the
five independent samples.

\subsection{ANALYSIS}

Of primary interest when comparing our simulations to experiments is
the dynamics of the low-frequency breathing mode of the particles.  To
study this motion, we need access to accurate measures of both the
volume and surface area as a function of time.  Throughout the
simulations, we monitored both quantities using Barber {\it et al.}'s
very fast quickhull algorithm to obtain the convex hull for the
collection of 3-d coordinates of all of atoms at each point in
time.~\cite{barber96quickhull,qhull} The convex hull is the smallest convex
polyhedron which includes all of the atoms, so the volume of this
polyhedron is an excellent estimate of the volume of the nanoparticle.
The convex hull estimate of the volume will be problematic if the
nanoparticle breaks into pieces (i.e. if the bounding surface becomes
concave), but for the relatively short trajectories of this study, it
provides an excellent measure of particle volume as a function of
time.

The bulk modulus, which is the inverse of the compressibility,
\begin{equation}
K = \frac{1}{\kappa} = - V \left(\frac{\partial P}{\partial
V}\right)_{T} 
\label{eq:Kdef}
\end{equation}
can be related to quantities that are relatively easily accessible
from our molecular dynamics simulations.  We present here four
different approaches for estimating the bulk modulus directly from
basic or derived quantities:

\begin{itemize}
\item The traditional ``Energetic'' approach

Using basic thermodynamics and one of the Maxwell relations, we write
\begin{equation}
P = T\left(\frac{\partial S}{\partial V}\right)_{T} -
\left(\frac{\partial U}{\partial V}\right)_{T}.
\label{eq:Pdef}
\end{equation}
It follows that 
\begin{equation}
K = V \left(T \left(\frac{\partial^{2} S}{\partial V^{2}}\right)_{T} -
\left(\frac{\partial^{2} U}{\partial V^{2}}\right)_{T} \right).
\label{eq:Ksub}
\end{equation}
The standard practice in solid state physics is to assume the low
temperature limit ({\it i.e.} to neglect the entropic term), which
means the bulk modulus is usually expressed
\begin{equation}
K \approx V \left(\frac{\partial^{2} U }{\partial V^{2}}\right)_{T}.
\label{eq:Kuse}
\end{equation}
When the relationship between the total energy $U$ and volume $V$ of
the system are available at a fixed temperature (as it is in these
simulations), it is a simple matter to compute the bulk modulus from
the response of the system to the perturbation of the instantaneous
heating.  Although this information would, in theory, be available
from longer constant temperature simulations, the ranges of volumes
and energies explored by a nanoparticle under equilibrium conditions
are actually quite small.  Instantaneous heating, since it excites
coherent oscillations in the breathing mode, allows us to sample a
much wider range of volumes (and energies) for the particles.  The
problem with this method is that it neglects the entropic term near
the melting transition, which gives spurious results (negative bulk
moduli) at higher temperatures.

\item The Linear Response Approach

Linear response theory gives us another approach to calculating the
bulk modulus.  This method relates the low-wavelength fluctuations in
the density to the isothermal compressibility,\cite{BernePecora}
\begin{equation}
\underset{k \rightarrow 0}{\lim} \langle | \delta \rho(\vec{k}) |^2
\rangle  =  \frac{k_B T \rho^2 \kappa}{V}
\end{equation}
where the frequency dependent density fluctuations are the Fourier
transform of the spatial fluctuations,
\begin{equation}
\delta \rho(\vec{k})  =  \underset{V}{\int} e^{i \vec{k}\cdot\vec{r}} \left( \rho(\vec{r}, t) - \langle \rho \rangle \right) dV
\end{equation}
This approach is essentially equivalent to using the volume
fluctuations directly to estimate the bulk modulus,
\begin{equation}
K  = \frac{V k_B T}{\langle \delta V^2 \rangle} =  k_B T \frac{\langle V \rangle}{\langle V^2 \rangle - \langle V \rangle^2}
\end{equation}
It should be noted that in these experiments, the particles are {\it
far from equilibrium}, and so a linear response approach will not be
the most suitable way to obtain estimates of the bulk modulus.

\item The Extended System Approach

Since we are performing these simulations in the NVT ensemble using
Nos\'e-Hoover thermostatting, the quantity conserved by our integrator
($H_{NVT}$) can be expressed as:
\begin{equation}
H_{NVT}  =  U  + f k_B T_{ext} \left( \frac{\tau_T^2 \chi(t)^2}{2} +
\int_0^t \chi(s) ds \right).
\end{equation}
Here $f$ is the number of degrees of freedom in the (real) system,
$T_{ext}$ is the temperature of the thermostat, $\tau_T$ is the time
constant for the thermostat, and $\chi(t)$ is the instantaneous value
of the extended system thermostat variable.  The extended Hamiltonian
for the system, $H_{NVT}$ is, to within a constant, the Helmholtz free
energy.\cite{melchionna93} Since the pressure is a simple derivative
of the Helmholtz free energy,
\begin{equation}
P  = -\left( \frac{\partial A}{\partial V} \right)_T , 
\end{equation}
the bulk modulus can be obtained (theoretically) by a quadratic fit of
the fluctuations in $H_{NVT}$ against fluctuations in the volume,
\begin{equation}
K =  -V \left( \frac{\partial^2 H_{NVT}}{\partial V^2} \right)_T.
\end{equation}
However, $H_{NVT}$ is essentially conserved during these simulations,
so fitting fluctuations of this quantity to obtain meaningful physical
quantities is somewhat suspect.  We also note that this method would
fail in periodic systems because the volume itself is fixed in
periodic NVT simulations.

\item The Direct Pressure Approach

Our preferred method for estimating the bulk modulus is to compute it
{\it directly} from the internal pressure in the nanoparticle.  The
pressure is obtained via the trace of the pressure tensor,
\begin{equation}
P  =  \frac{1}{3}
\mathrm{Tr}\left[\overleftrightarrow{\mathsf{P}}\right],
\end{equation}
which has a kinetic contribution as well as a contribution from the
stress tensor ($\overleftrightarrow{\mathsf{W}}$):
\begin{equation}
\overleftrightarrow{\mathsf{P}}  =  \frac{1}{V} \left(
\sum_{i=1}^{N} m_i \vec{v}_i \otimes \vec{v}_i  \right) +
\overleftrightarrow{\mathsf{W}}.
\end{equation}
Here the $\otimes$ symbol represents the {\it outer} product of the
velocity vector for atom $i$ to yield a $3 \times 3$ matrix. The
virial is computed during the simulation using forces between pairs of
particles,
\begin{equation}
\overleftrightarrow{\mathsf{W}}  =  \sum_{i} \sum_{j>i}  \vec{r}_{ij} \otimes \vec{f}_{ij} 
\end{equation}
During the simulation, we record the internal pressure, $P$, as well as
the total energy, $U$, the extended system's hamiltonian, $H_{NVT}$,
and the particle coordinates.  Once we have calculated the time
dependent volume of the nanoparticle using the convex hull, we can use
any of these four methods to estimate the bulk modulus.
\end{itemize}

We find, however, that only the fourth method (the direct pressure
approach) gives meaningful results.  Bulk moduli for the 35 \AA\
particle were computed with the traditional (Energy vs. Volume)
approach as well as the direct pressure approach.  A comparison of the
Bulk Modulus obtained via both methods and are shown in
Fig. \ref{fig:Methods} Note that the second derivative fits in the
traditional approach can give (in the liquid droplet region) negative
curvature, and this results in negative values for the bulk modulus.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/Methods.pdf}
        \caption{Comparison of two of the methods for estimating the bulk
                modulus as a function of temperature for the 35\AA\ particle.}
        \label{fig:Methods}
\end{figure}

The Bulk moduli reported in the rest of this paper were computed using
the direct pressure method.

To study the frequency of the breathing mode, we have calculated the
power spectrum for volume ($V$) fluctuations,
\begin{equation}
\rho_{\Delta V}(\omega) = \int_{-\infty}^{\infty} \langle \Delta V(t)
\Delta V(0) \rangle e^{-i \omega t} dt
\label{eq:volspect}
\end{equation}
where $\Delta V(t) = V(t) - \langle V \rangle$.  Because the
instantaneous heating excites all of the vibrational modes of the
particle, the power spectrum will contain contributions from all modes
that perturb the total volume of the particle.  The lowest frequency
peak in the power spectrum should give the frequency (and period) for
the breathing mode, and these quantities are most readily compared
with the Hartland group experiments.\cite{HartlandG.V._jp0276092} Further
analysis of the breathing dynamics follows in section
\ref{bulkmod:sec:discussion}.

We have also computed the heat capacity for our simulations to verify
the location of the melting transition.  Calculations of the heat
capacity were performed on the non-equilibrium, instantaneous heating
simulations, as well as on simulations of nanoparticles that were at
equilibrium at the target temperature.
   
\section{RESULTS}
\label{bulkmod:sec:results}

\subsection{THE BULK MODULUS AND HEAT CAPACITY}

The upper panel in Fig. \ref{fig:BmCp} shows the temperature
dependence of the Bulk Modulus ($K$).  In all samples, there is a
dramatic (size-dependent) drop in $K$ at temperatures well below the
melting temperature of bulk polycrystalline gold.  This drop in $K$
coincides with the actual melting transition of the nanoparticles.
Surface melting has been confirmed at even lower temperatures using
the radial-dependent density, $\rho(r) / \rho$, which shows a merging
of the crystalline peaks in the outer layer of the nanoparticle.
However, the bulk modulus only has an appreciable drop when the
particle melts fully.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/Stacked_Bulk_modulus_and_Cp.pdf}
        \caption{The temperature dependence of the bulk modulus (upper panel)
                and heat capacity (lower panel) for nanoparticles of four different
                radii. Note that the peak in the heat capacity coincides with the {\em
                start} of the peak in the bulk modulus.}
        \label{fig:BmCp}
\end{figure}


Another feature of these transient (non-equilibrium) calculations is
the width of the peak in the heat capacity.  Calculation of $C_{p}$
from longer equilibrium trajectories should indicate {\it sharper}
features in $C_{p}$ for the larger particles.  Since we are initiating
and observing the melting process itself in these calculation, the
smaller particles melt more rapidly, and thus exhibit sharper features
in $C_{p}$.  Indeed, longer trajectories do show that $T_{m}$ occurs
at lower temperatures and with sharper transitions in larger particles
than can be observed from transient response calculations.
Fig. \ref{fig:Cp2} shows the results of 300 ps simulations which give
much sharper and lower temperature melting transitions than those
observed in the 40 ps simulations.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/Cp_vs_T.pdf}
        \caption{The dependence of the spike in the heat capacity on the
                length of the simulation.  Longer heating-response calculations result
                in melting transitions that are sharper and lower in temperature than
                the short-time transient response simulations.  Shorter runs don't
                allow the particles to melt completely.}
        \label{fig:Cp2}
\end{figure}


\subsection{BREATHING MODE DYNAMICS}

Fig.\ref{fig:VolTime} shows representative samples of the volume
vs. time traces for the 20 \AA\ and 35 \AA\ particles at a number of
different temperatures.  We can clearly see that the period of the
breathing mode is dependent on temperature, and that the coherent
oscillations of the particles' volume are destroyed after only a few
ps in the smaller particles, while they live on for 10-20 ps in the
larger particles.  The de-coherence is also strongly temperature
dependent, with the high temperature samples decohering much more
rapidly than lower temperatures.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/Vol_vs_time.pdf}
        \caption{Sample Volume traces for the 20 \AA\ and 35 \AA\ particles at a
                range of temperatures.  Note the relatively rapid ($<$ 10 ps)
                decoherence due to melting in the 20 \AA\ particle as well as the
                difference between the 1100 K and 1200 K traces in the 35 \AA\
                particle.}
        \label{fig:VolTime}
\end{figure}


Although $V$ vs. $t$ traces can say a great deal, it is more
instructive to compute the autocorrelation function for volume
fluctuations to give more accurate short-time information.  Fig
\ref{fig:volcorr} shows representative autocorrelation functions for
volume fluctuations.  Although many traces exhibit a single frequency
with decaying amplitude, a number of the samples show distinct beat
patterns indicating the presence of multiple frequency components in
the breathing motion of the nanoparticles.  In particular, the 20 \AA\
particle shows a distinct beat in the volume fluctuations in the 800 K
trace.


\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/volcorr.pdf}
        \caption{Volume fluctuation autocorrelation functions for the 20 \AA\
                and 35 \AA\ particles at a range of temperatures.  Successive
                temperatures have been translated upwards by one unit.  Note the beat
                pattern in the 20 \AA\ particle at 800K.}
        \label{fig:volcorr}
\end{figure}


When the power spectrum of the volume autocorrelation functions are
analyzed (Eq. (\ref{eq:volspect})), the samples which exhibit beat
patterns do indeed show multiple peaks in the power spectrum.  We plot
the period corresponding to the two lowest frequency peaks in
Fig. \ref{fig:Period}.  The smaller particles have the most evident
splitting, particularly as the temperature rises above the melting
points for these particles.

\begin{figure}[htbp]
        \centering
                \includegraphics[height=3in]{images/Period_vs_T.pdf}
        \caption{The temperature dependence of the period of the breathing
                mode for the four different nanoparticles studied in this
                work.}
        \label{fig:Period}
\end{figure}


\section{DISCUSSION}
\label{bulkmod:sec:discussion}

Lamb's classical theory of elastic spheres~\cite{Lamb1882} provides
two possible explanations for the split peak in the vibrational
spectrum.  The periods of the longitudinal and transverse vibrations
in an elastic sphere of radius $R$ are given by:
\begin{equation}
\tau_{t} = \frac{2 \pi R}{\theta c_{t}}
\end{equation}
and
\begin{equation}
\tau_{l} = \frac{2 \pi R}{\eta c_{l}}
\end{equation}
where $\theta$ and $n$ are obtained from the solutions to the
transcendental equations
\begin{equation}
\tan \theta = \frac{3 \theta}{3 - \theta^{2}}
\end{equation}
\begin{equation}
\tan \eta = \frac{4 \eta}{4 - \eta^{2}\frac{c_{l}^{2}}{c_{t}^{2}}}
\end{equation}.

$c_{l}$ and $c_{t}$ are the longitudinal and transverse speeds
of sound in the material.  In an isotropic material, these speeds are
simply related to the elastic constants and the density ($\rho$),
\begin{equation}
c_{l} = \sqrt{c_{11}/\rho}
\end{equation}
\begin{equation}
c_{t} = \sqrt{c_{44}/\rho}.
\end{equation}

In crystalline materials, the speeds depend on the direction of
propagation of the wave relative to the crystal plane.\cite{Kittel:1996fk}
For the remainder of our analysis, we assume the nanoparticles are
isotropic (which should be valid only above the melting transition).
A more detailed analysis of the lower temperature particles would take
the crystal lattice into account.

If we use the experimental values for the elastic constants for 30
\AA\ Au particles at 300K, the low-frequency longitudinal (breathing)
mode should have a period of 2.19 ps while the low-frequency
transverse (toroidal) mode should have a period of 2.11 ps.  Although
the actual calculated frequencies in our simulations are off of these
values, the difference in the periods (0.08 ps) is approximately half
of the splitting observed room-temperature simulations.  This,
therefore, may be an explanation for the low-temperature splitting in
Fig. \ref{fig:Period}.

We note that Cerullo {\it et al.} used a similar treatment to obtain
the low frequency longitudinal frequencies for crystalline
semiconductor nanoparticles,\cite{Cerullo1999} and Simon and Geller
have investigated the effects of ensembles of particle size on the
Lamb mode using the classical Lamb theory results for isotropic
elastic spheres.\cite{Simon2001}

\subsection{MELTED AND PARTIALLY-MELTED PARTICLES}

Hartland {\it et al.} have extended the Lamb analysis to include
surface stress ($\gamma$).\cite{HartlandG.V._jp0276092} In this case, the
transcendental equation that must be solved to obtain the
low-frequency longitudinal mode is
\begin{equation}
\eta \cot \eta = 1 - \frac{\eta^{2} c_{l}^{2}}{4 c_{l}^{2} -
2 \gamma / (\rho R)}.
\end{equation}
In ideal liquids, inclusion of the surface stress is vital since the
transverse speed of sound ($c_{t}$) vanishes.  Interested readers
should consult Hartland {\it et al.}'s paper for more details on the
extension to liquid-like particles, but the primary result is that the
vibrational period of the breathing mode for liquid droplets may be
written
\begin{equation}
\tau = \frac{2 R}{c_{l}(l)}
\end{equation}
where $c_{l}(l)$ is the longitudinal speed of sound in the liquid.
Iida and Guthrie list the speed of sound in liquid Au metal as 
\begin{equation}
c_{l}(l) = 2560 - 0.55 (T - T_{m})  (\mbox{m s}^{-1})
\end{equation}
where $T_{m}$ is the melting temperature.\cite{Iida1988} A molten 35
\AA\ particle just above $T_{m}$ would therefore have a vibrational
period of 2.73 ps, and this would be markedly different from the
vibrational period just below $T_{m}$ if the melting transition were
sharp.  

We know from our calculations of $C_{p}$ that the complete melting of
the particles is {\it not} sharp, and should take longer than the 40
ps observation time.  There are therefore two explanations which are
commensurate with our observations.
\begin{enumerate}
\item The melting may occur at some time partway through observation
of the response to instantaneous heating.  The early part 
of the simulation would then show a higher-frequency breathing mode
than would be evident during the latter parts of the simulation. 
\item The melting may take place by softening the outer layers of the
particle first, followed by a melting of the core at higher
temperatures.  The liquid-like outer layer would then contribute a
lower frequency component than the interior of the particle.
\end{enumerate}

The second of these explanations is consistent with the core-shell
melting hypothesis advanced by Hartland {\it et al.} to explain their
laser heating experiments.\cite{HartlandG.V._jp0276092} At this stage, our
simulations cannot distinguish between the two hypotheses.  One
possible avenue for future work would be the computation of a
radial-dependent order parameter to help evaluate whether the
solid-core/liquid-shell structure exists in our simulation.
Revision:	3483
Committed:	Tue Jan 13 14:39:50 2009 UTC (15 years, 8 months ago) by chuckv
Content type:	application/x-tex
File size:	24630 byte(s)
Log Message:	Chuck's dissertation for PhD Jan 2009