trunk/chainLength/GoldThiolsPaper.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{times}
\usepackage{mathptm}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
\usepackage{tabularx}
\usepackage{longtable}
\usepackage{graphicx}
\usepackage{multirow}
\usepackage{multicol}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures
%\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

\bibpunct{}{}{,}{s}{}{;}

\citestyle{nature}
\bibliographystyle{achemso}

\begin{document}

\newcolumntype{A}{p{1.5in}}
\newcolumntype{B}{p{0.75in}}

\title{Simulations of Heat Conduction at Thiolate-Capped Gold
  Surfaces: The Role of Chain Length and Solvent Penetration}

\author{Kelsey M. Stocker and J. Daniel
  Gezelter\footnote{Corresponding author. \ Electronic mail:
    gezelter@nd.edu} \\
  251 Nieuwland Science Hall, \\
        Department of Chemistry and Biochemistry,\\
        University of Notre Dame\\
        Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}

  We report on simulations of heat conduction through Au(111) / hexane
  interfaces in which the surface has been protected by a mixture of
  short and long chain alkanethiolate ligands.  Reverse
  non-equilibrium molecular dynamics (RNEMD) was used to create a
  thermal flux between the metal and solvent, and thermal conductance
  was computed using the resulting thermal profiles near the
  interface.  We find a non-monotonic dependence of the interfacial
  thermal conductance on the fraction of long chains present at the
  interface, and correlate this behavior to both solvent ordering and
  the rate of solvent escape from the thiolate layer immediately in
  contact with the metal surface.  Our results suggest that a mixed
  vibrational transfer / convection model is necessary to explain the
  features of heat transfer at this interface. The alignment of the
  solvent chains with the ordered ligand allows rapid transfer of
  energy to the trapped solvent and is the dominant feature for
  ordered ligand layers. Diffusion of the vibrationally excited
  solvent into the bulk also plays a significant role when the ligands
  are less tightly packed.

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **INTRODUCTION**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

The structural details of the interfaces of metal nanoparticles
determine how energy flows between these particles and their
surroundings. Understanding this energy flow is essential in designing
and functionalizing metallic nanoparticles for use in plasmonic photothermal
therapies,\cite{Jain:2007ux,Petrova:2007ad,Gnyawali:2008lp,Mazzaglia:2008to,Huff:2007ye,Larson:2007hw}
which rely on the ability of metallic nanoparticles to absorb light in
the near-IR, a portion of the spectrum in which living tissue is very
nearly transparent.  The relevant physical property controlling the
transfer of this energy as heat into the surrounding tissue is the
interfacial thermal conductance, $G$, which can be somewhat difficult
to determine experimentally.\cite{Wilson:2002uq,Plech:2005kx}

Metallic particles have also been proposed for use in  efficient
thermal-transfer fluids, although careful experiments by Eapen {\it et
  al.}  have shown that metal-particle-based nanofluids have thermal
conductivities that match Maxwell predictions.\cite{Eapen:2007th} The
likely cause of previously reported non-Maxwell
behavior\cite{Eastman:2001wb,Keblinski:2002bx,Lee:1999ct,Xue:2003ya,Xue:2004oa}
is percolation networks of nanoparticles exchanging energy via the
solvent,\cite{Eapen:2007mw} so it is important to get a detailed
molecular picture of particle-ligand and ligand-solvent interactions
in order to understand the thermal behavior of complex fluids. To
date, there have been few reported values (either from theory or
experiment) of $G$ for ligand-protected nanoparticles embedded in
liquids, and there is a significant gap in knowledge about how
chemically distinct ligands or protecting groups will affect heat
transport from the particles.

Experimentally, the thermal properties of various nanostructured
interfaces have been investigated by a number of groups. Cahill and
coworkers studied thermal transport from metal nanoparticle/fluid
interfaces, epitaxial TiN/single crystal oxide interfaces, and
hydrophilic and hydrophobic interfaces between water and solids with
different self-assembled
monolayers.\cite{cahill:793,Wilson:2002uq,PhysRevB.67.054302,doi:10.1021/jp048375k,PhysRevLett.96.186101} Schmidt {\it et al.} studied the role of
cetyltrimethylammonium bromide (CTAB) on the thermal transport between
gold nanorods and solvent.\cite{doi:10.1021/jp8051888}
Wang {\it et al.} studied heat transport through long-chain
hydrocarbon monolayers on unsolvated gold substrate at the individual molecular
level.\cite{Wang10082007} The introduction of solvent adds yet another interface and potential barrier for heat transfer. Juv\'e {\it
  et al.} studied the cooling dynamics of glass-embedded nanoparticles, which is controlled by thermal
interfacial resistance and heat diffusion in the matrix.\cite{PhysRevB.80.195406} Hartland also notes the importance of heat transfer through diffusive solvent behavior.\cite{hartland2011} Although interfaces are
normally considered barriers for heat transport, Alper {\it et al.}
have suggested that specific ligands (capping agents) could completely
eliminate this barrier
($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}

Recently, Hase and coworkers employed Non-Equilibrium Molecular
Dynamics (NEMD) simulations to study thermal transport from hot
Au(111) substrate to a self-assembled monolayer of alkylthiol with
relatively long chain (8-20 carbon atoms).\cite{hase:2010,hase:2011}
These simulations explained many of the features of the experiments of
Wang {\it et al.}  However, ensemble averaged measurements for heat
conductance of interfaces between the capping monolayer on Au and a
solvent phase have yet to be studied with their approach.  In previous
simulations, our group applied a variant of reverse non-equilibrium
molecular dynamics (RNEMD) to calculate the interfacial thermal
conductance at a metal / organic solvent interface that had been
chemically protected by butanethiolate groups.\cite{kuang:AuThl} Our
calculations suggested an explanation for the very large thermal
conductivity at alkanethiol-capped metal surfaces when compared with
bare metal/solvent interfaces.  Specifically, the chemical bond
between the metal and the ligand introduces a vibrational overlap that
is not present without the protecting group, and the overlap between
the vibrational spectra (metal to ligand, ligand to solvent) provides
a mechanism for rapid thermal transport across the interface.

A notable result of the previous simulations was the non-monotonic
dependence of $G$ on the fractional coverage of the metal surface by
the chemical protecting group.  Gaps in surface coverage allow the
solvent molecules to come into direct contact with ligands that have
been heated by the metal surface, absorb thermal energy from the
ligands, and then diffuse away.  Quantifying the role of overall
surface coverage was difficult because the ligands have lateral
mobility on the surface and can aggregate to form domains on the
timescale of the simulation.

To isolate the effect of ligand/solvent coupling while avoiding
lateral mobility of the surface ligands, the current work utilizes
monolayers of mixed chain-lengths in which the length mismatch between
long and short chains is sufficient to accommodate solvent
molecules. These completely covered (but mixed-chain) surfaces are
significantly less prone to surface domain formation, and allow us to
further investigate the mechanism of heat transport to the solvent.  A
thermal flux is created using velocity shearing and scaling reverse
non-equilibrium molecular dynamics (VSS-RNEMD), and the resulting
temperature profiles are analyzed to yield information about the
interfacial thermal conductance.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **METHODOLOGY**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Methodology}

There are many ways to compute bulk transport properties from
classical molecular dynamics simulations.  Equilibrium Molecular
Dynamics (EMD) simulations can be used to compute the relevant time
correlation functions and transport coefficients can be calculated
assuming that linear response theory
holds.\cite{PhysRevB.37.5677,MASSOBRIO:1984bl,PhysRev.119.1,Viscardy:2007rp,che:6888,kinaci:014106}
For some transport properties (notably thermal conductivity), EMD
approaches are subject to noise and poor convergence of the relevant
correlation functions. Traditional Non-equilibrium Molecular Dynamics
(NEMD) methods impose a gradient (e.g. thermal or momentum) on a
simulation.\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2002dp,PhysRevA.34.1449,JiangHao_jp802942v}
However, the resulting flux is often difficult to
measure. Furthermore, problems arise for NEMD simulations of
heterogeneous systems, such as phase-phase boundaries or interfaces,
where the type of gradient to enforce at the boundary between
materials is unclear.

{\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt
a different approach in that an unphysical {\it flux} is imposed
between different regions or ``slabs'' of the simulation
box.\cite{MullerPlathe:1997xw,ISI:000080382700030,Kuang2010} The
system responds by developing a temperature or momentum {\it gradient}
between the two regions. Since the amount of the applied flux is known
exactly, and the measurement of a gradient is generally less
complicated, imposed-flux methods typically take shorter simulation
times to obtain converged results for transport properties.  The
corresponding temperature or velocity gradients which develop in
response to the applied flux are then related (via linear response
theory) to the transport coefficient of interest.  These methods are
quite efficient, in that they do not need many trajectories to provide
information about transport properties. To date, they have been
utilized in computing thermal and mechanical transfer of both
homogeneous
liquids~\cite{MullerPlathe:1997xw,ISI:000080382700030,Maginn:2010} as
well as heterogeneous
systems.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% VSS-RNEMD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{VSS-RNEMD}
The original ``swapping'' approaches by M\"{u}ller-Plathe {\it et
  al.}\cite{ISI:000080382700030,MullerPlathe:1997xw} can be understood
as a sequence of imaginary elastic collisions between particles in
regions separated by half of the simulation cell.  In each collision,
the entire momentum vectors of both particles may be exchanged to
generate a thermal flux. Alternatively, a single component of the
momentum vectors may be exchanged to generate a shear flux.  This
algorithm turns out to be quite useful in many simulations of bulk
liquids. However, the M\"{u}ller-Plathe swapping approach perturbs the
system away from ideal Maxwell-Boltzmann distributions, which has
undesirable side-effects when the applied flux becomes
large.\cite{Maginn:2010}        

The most useful alternative RNEMD approach developed so far utilizes a
series of simultaneous velocity shearing and scaling (VSS) exchanges between
the two slabs.\cite{Kuang2012} This method provides a set of
conservation constraints while simultaneously creating a desired flux
between the two slabs.  Satisfying the constraint equations ensures
that the new configurations are sampled from the same NVE ensemble.

The VSS moves are applied periodically to scale and shift the particle
velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and
$C$) which are separated by half of the simulation box,
\begin{displaymath}
\begin{array}{rclcl}

 & \underline{\mathrm{shearing}} & &
 \underline{~~~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~~~} \\  \\
\mathbf{v}_i \leftarrow & 
  \mathbf{a}_c\ & + & c\cdot\left(\mathbf{v}_i - \langle\mathbf{v}_c
  \rangle\right)  +  \langle\mathbf{v}_c\rangle \\
\mathbf{v}_j \leftarrow & 
  \mathbf{a}_h & + & h\cdot\left(\mathbf{v}_j - \langle\mathbf{v}_h
    \rangle\right) + \langle\mathbf{v}_h\rangle .

\end{array}
\end{displaymath}
Here $\langle\mathbf{v}_c\rangle$ and $\langle\mathbf{v}_h\rangle$ are
the center of mass velocities in the $C$ and $H$ slabs, respectively.
Within the two slabs, particles receive incremental changes or a
``shear'' to their velocities.  The amount of shear is governed by the
imposed momentum flux, $\mathbf{j}_z(\mathbf{p})$
\begin{eqnarray}
\mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\
\mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2}
\end{eqnarray}
where $M_{\{c,h\}}$ is the total mass of particles within each of the
slabs and $\Delta t$ is the interval between two separate operations.

To simultaneously impose a thermal flux ($J_z$) between the slabs we
use energy conservation constraints,
\begin{eqnarray}
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\mathbf{v}_c
\rangle^2) + \frac{1}{2}M_c (\langle \mathbf{v}_c \rangle + \mathbf{a}_c)^2 \label{vss3}\\
K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\mathbf{v}_h
\rangle^2) + \frac{1}{2}M_h (\langle \mathbf{v}_h \rangle +
\mathbf{a}_h)^2 \label{vss4}.
\label{constraint}
\end{eqnarray}
Simultaneous solution of these quadratic formulae for the scaling
coefficients, $c$ and $h$, will ensure that the simulation samples from
the original microcanonical (NVE) ensemble.  Here $K_{\{c,h\}}$ is the
instantaneous translational kinetic energy of each slab.  At each time
interval, it is a simple matter to solve for $c$, $h$, $\mathbf{a}_c$,
and $\mathbf{a}_h$, subject to the imposed momentum flux,
$j_z(\mathbf{p})$, and thermal flux, $J_z$, values.  Since the VSS
operations do not change the kinetic energy due to orientational
degrees of freedom or the potential energy of a system, configurations
after the VSS move have exactly the same energy (and linear
momentum) as before the move.

As the simulation progresses, the VSS moves are performed on a regular
basis, and the system develops a thermal or velocity gradient in
response to the applied flux.  Using the slope of the temperature or
velocity gradient, it is quite simple to obtain the thermal
conductivity ($\lambda$),
\begin{equation}
J_z = -\lambda \frac{\partial T}{\partial z},
\end{equation}
and shear viscosity ($\eta$),
\begin{equation}
j_z(p_x) = -\eta \frac{\partial \langle v_x \rangle}{\partial z}.
\end{equation}
Here, the quantities on the left hand side are the imposed flux
values, while the slopes are obtained from linear fits to the
gradients that develop in the simulated system.

The VSS-RNEMD approach is versatile in that it may be used to
implement both thermal and shear transport either separately or
simultaneously.  Perturbations of velocities away from the ideal
Maxwell-Boltzmann distributions are minimal, as is thermal anisotropy.
This ability to generate simultaneous thermal and shear fluxes has
been previously utilized to map out the shear viscosity of SPC/E water
over a wide range of temperatures (90~K) with a single 1 ns
simulation.\cite{Kuang2012}

\begin{figure}
  \includegraphics[width=\linewidth]{figures/rnemd}
  \caption{The VSS-RNEMD approach imposes unphysical transfer of
    linear momentum or kinetic energy between a ``hot'' slab and a
    ``cold'' slab in the simulation box.  The system responds to this
    imposed flux by generating velocity or temperature gradients.  The
    slope of the gradients can then be used to compute transport
    properties (e.g. shear viscosity or thermal conductivity).}
  \label{fig:rnemd}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% INTERFACIAL CONDUCTANCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Reverse Non-Equilibrium Molecular Dynamics approaches
  to interfacial transport}

Interfaces between dissimilar materials have transport properties
which can be defined as derivatives of the standard transport
coefficients in a direction normal to the interface. For example, the
{\it interfacial} thermal conductance ($G$) can be thought of as the
change in the thermal conductivity ($\lambda$) across the boundary
between materials:
\begin{align}
  G &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\
  &= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/
  \left(\frac{\partial T}{\partial z}\right)_{z_0}^2.
  \label{derivativeG}
\end{align}
where $z_0$ is the location of the interface between two materials and
$\mathbf{\hat{n}}$ is a unit vector normal to the interface (assumed
to be the $z$ direction from here on).  RNEMD simulations, and
particularly the VSS-RNEMD approach, function by applying a momentum
or thermal flux and watching the gradient response of the
material. This means that the {\it interfacial} conductance is a
second derivative property which is subject to significant noise and
therefore requires extensive sampling.  Previous work has demonstrated
the use of the second derivative approach to compute interfacial
conductance at chemically-modified metal / solvent interfaces.
However, a definition of the interfacial conductance in terms of a
temperature difference ($\Delta T$) measured at $z_0$,
\begin{displaymath}
G = \frac{J_z}{\Delta T_{z_0}},
\end{displaymath}
is useful once the RNEMD approach has generated a stable temperature
gap across the interface.

\begin{figure}
  \includegraphics[width=\linewidth]{figures/resistor_series}
  \caption{The inverse of the interfacial thermal conductance, $G$, is
    the Kapitza resistance, $R_K$.  Because the gold / thiolate /
    solvent interface extends a significant distance from the metal
    surface, the interfacial resistance $R_K$ can be computed by
    summing a series of temperature drops between adjacent temperature
    bins along the $z$ axis. The depicted temperature profile is from a RNEMD simulation of 100\% butanethiolate (C$_4$) coverage.}
  \label{fig:resistor_series}
\end{figure}

In the particular case we are studying here, there are two interfaces
involved in the transfer of heat from the gold slab to the solvent:
the metal / thiolate interface and the thiolate / solvent interface. We
can treat the temperature on each side of an interface as discrete,
making the interfacial conductance the inverse of the Kaptiza
resistance, or $G = \frac{1}{R_k}$. To model the total conductance
across multiple interfaces, it is useful to think of the interfaces as
a set of resistors wired in series. The total resistance is then
additive, $R_{total} = \sum_i R_{i}$, and the interfacial conductance
is the inverse of the total resistance, or $G = \frac{1}{\sum_i
  R_i}$).  In the interfacial region, we treat each bin in the
VSS-RNEMD temperature profile as a resistor with resistance
$\frac{T_{2}-T_{1}}{J_z}$, $\frac{T_{3}-T_{2}}{J_z}$, etc.  The sum of
the set of resistors which spans the gold / thiolate interface, thiolate
chains, and thiolate / solvent interface simplifies to
\begin{equation}
  R_{K} = \frac{T_{n}-T_{1}}{J_z},
  \label{eq:finalG}
\end{equation} 
or the temperature difference between the gold side of the
gold / thiolate interface and the solvent side of the thiolate / solvent
interface over the applied flux.  
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **COMPUTATIONAL DETAILS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\section{Computational Details}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% FORCE-FIELD PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Force-Field Parameters}

Our simulations include a number of chemically distinct components.
Figure \ref{fig:structures} demonstrates the sites defined for both
the {\it n}-hexane and alkanethiolate ligands present in our
simulations.

\begin{figure}
  \includegraphics[width=\linewidth]{figures/structures}
  \caption{Topologies of the thiolate capping agents and solvent
    utilized in the simulations. The chemically-distinct sites (S,
    \ce{CH2}, and \ce{CH3}) are treated as united atoms. Most
    parameters are taken from references \bibpunct{}{}{,}{n}{}{,}
    \protect\cite{TraPPE-UA.alkanes} and
    \protect\cite{TraPPE-UA.thiols}. Cross-interactions with the Au
    atoms were adapted from references
    \protect\cite{landman:1998},~\protect\cite{vlugt:cpc2007154},~and
    \protect\cite{hautman:4994}.}
  \label{fig:structures}
\end{figure}

The Au-Au interactions in the metal lattice slab were described by the
quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
potentials include zero-point quantum corrections and are
reparametrized for accurate surface energies compared to the
Sutton-Chen potentials.\cite{Chen90}

For the {\it n}-hexane solvent molecules, the TraPPE-UA
parameters\cite{TraPPE-UA.alkanes} were utilized.  In this model,
sites are located at the carbon centers for alkyl groups. Bonding
interactions, including bond stretches and bends and torsions, were
used for intra-molecular sites closer than 3 bonds. For non-bonded
interactions, Lennard-Jones potentials were used.  We have previously
utilized both united atom (UA) and all-atom (AA) force fields for
thermal conductivity,\cite{kuang:AuThl,Kuang2012} and since the united
atom force fields cannot populate the high-frequency modes that are
present in AA force fields, they appear to work better for modeling
thermal conductivity.  The TraPPE-UA model for alkanes is known to
predict a slightly lower boiling point than experimental values. This
is one of the reasons we used a lower average temperature (220 K) for
our simulations.

The TraPPE-UA force field includes parameters for thiol
molecules\cite{TraPPE-UA.thiols} which were used for the
alkanethiolate molecules in our simulations.  To derive suitable
parameters for butanethiolate adsorbed on Au(111) surfaces, we adopted
the S parameters from Luedtke and Landman\cite{landman:1998} and
modified the parameters for the CTS atom to maintain charge neutrality
in the molecule. 

To describe the interactions between metal (Au) and non-metal atoms,
we refer to an adsorption study of alkyl thiols on gold surfaces by
Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
Lennard-Jones form of potential parameters for the interaction between
Au and pseudo-atoms CH$_x$ and S based on a well-established and
widely-used effective potential of Hautman and Klein for the Au(111)
surface.\cite{hautman:4994} As our simulations require the gold slab
to be flexible to accommodate thermal excitation, the pair-wise form
of potentials they developed was used for our study.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% SIMULATION PROTOCOL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Simulation Protocol}

We have implemented the VSS-RNEMD algorithm in OpenMD, our group
molecular dynamics code.\cite{openmd} An 1188 atom gold slab was
prepared and equilibrated at 1 atm and 200 K.  The periodic box was
then expanded in the $z$ direction to expose two Au(111) faces on
either side of the 11-layer slab.

A full monolayer of thiolates (1/3 the number of surface gold atoms)
was placed on three-fold hollow sites on the gold interfaces. The
effect of thiolate binding sites on the thermal conductance was tested
by placing thiolates at both fcc and hcp hollow sites.  No appreciable
difference in the temperature profile due to the location of
thiolate binding was noted.  

To test the role of thiolate chain length on interfacial thermal
conductance, full coverages of each of five chain lengths were tested:
butanethiolate (C$_4$), hexanethiolate (C$_6$), octanethiolate
(C$_8$), decanethiolate (C$_{10}$), and dodecanethiolate
(C$_{12}$). To test the effect of mixed chain lengths, full coverages
of C$_4$ / C$_{10}$ and C$_4$ / C$_{12}$ mixtures were created in
short/long chain percentages of 25/75, 50/50, 62.5/37.5, 75/25, and
87.5/12.5. The short and long chains were placed on the surface hollow
sites in a random configuration.

The simulation box $z$-dimension was set to roughly double the length
of the gold / thiolate block. Hexane solvent molecules were placed in
the vacant portion of the box using the packmol
algorithm.\cite{packmol} Figure \ref{fig:timelapse} shows two of the
mixed chain length interfaces both before and after the RNEMD simulation.

\begin{figure}
  \includegraphics[width=\linewidth]{figures/timelapse}
  \caption{Images of 75\%~C$_4$~/~25\%~C$_{12}$ (top panel) and 25\%~C$_4$~/~75\%~C$_{12}$ (bottom panel) interfaces at the beginning and end of 3 ns simulations. Solvent molecules that were initially present in the thiolate layer are colored light blue.  Diffusion of the initially-trapped solvent into the bulk is apparent in the interface with fewer long chains.  Trapped solvent is orientationally locked to the ordered ligands (and is less able to diffuse into the bulk) when the fraction of long chains increases.}
  \label{fig:timelapse}
\end{figure}

The system was equilibrated to 220 K in the canonical (NVT) ensemble,
allowing the thiolates and solvent to warm gradually. Pressure
correction to 1 atm was done using an isobaric-isothermal (NPT)
integrator that allowed expansion or contraction only in the $z$
direction, maintaining the crystalline structure of the gold as close
to the bulk result as possible.  The diagonal elements of the pressure
tensor were monitored during the pressure equilibration stage.  If the
$xx$ and/or $yy$ elements had a mean above zero throughout the
simulation -- indicating residual surface tension in the plane of the
gold slab -- an additional short NPT equilibration step was performed
allowing all box dimensions to change.  Once the pressure was stable
at 1 atm, a final equilibration stage was performed at constant
temperature. All systems were equilibrated in the microcanonical (NVE)
ensemble before proceeding with the VSS-RNEMD and data collection
stages.

A kinetic energy flux was applied using VSS-RNEMD during a data
collection period of 3 nanoseconds, with velocity scaling moves
occurring every 10 femtoseconds.  The ``hot'' slab was centered in the
gold and the ``cold'' slab was placed in the center of the solvent
region.  The entire system had a (time-averaged) temperature of 220 K,
with a temperature difference between the hot and cold slabs of
approximately 30 K.  The average temperature and kinetic energy flux
were selected to avoid solvent freezing (or glass formation) and to
prevent the thiolates from burying in the gold slab.  The Au-S
interaction has a deep potential energy well, which allows sulfur
atoms to embed into the gold slab at temperatures above 250 K.
Simulation conditions were chosen to avoid both of these 
situations.

Temperature profiles of the system were created by dividing the box
into $\sim$ 3 \AA \, bins along the $z$ axis and recording the average
temperature of each bin.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **RESULTS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                         
\section{Results}

The solvent, hexane, is a straight chain flexible alkane that is structurally
similar to the thiolate alkane tails. Previous work has shown that UA models
of hexane and butanethiolate have a high degree of vibrational
overlap.\cite{kuang:AuThl} This overlap provides a mechanism for thermal
energy conduction from the thiolates to the solvent. Indeed, we observe that
the interfacial conductance is twice as large with the thiolate monolayers (of
all chain lengths) than with the bare metal surface.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% CHAIN LENGTH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Effect of Chain Length}

We examined full coverages of five alkyl chain lengths, C$_{4}$, C$_{6}$, C$_{8}$, C$_{10}$, and C$_{12}$. As shown in table \ref{table:chainlengthG}, the interfacial conductance is roughly constant as a function of chain length, indicating that the length of the thiolate alkyl chains does not play a significant role in the transport of heat across the gold/thiolate and thiolate/solvent interfaces. In all cases, the hexane solvent was unable to penetrate into the thiolate layer, leading to a persistent 2-4 \AA \, gap between the solvent region and the thiolates. However, while the identity of the alkyl thiolate capping agent has little effect on the interfacial thermal conductance, the presence of a full monolayer of capping agent provides a two-fold increase in the G value relative to a bare gold surface.
\begin{longtable}{p{4cm} p{3cm}}
        \caption{Computed interfacial thermal conductance ($G$) values for bare gold and 100\% coverages of various thiolate alkyl chain lengths.}
        \\
        \centering {\bf Chain Length} & \centering\arraybackslash {\bf G (MW/m$^2$/K)} \\ \hline
\endhead
\hline
\endfoot
\centering bare metal & \centering\arraybackslash 30.2 \\
$~~~~~~~~~~~~~~~~~~$ C$_{4}$ & \centering\arraybackslash 59.4 \\
$~~~~~~~~~~~~~~~~~~$ C$_{6}$ & \centering\arraybackslash 60.2 \\
$~~~~~~~~~~~~~~~~~~$ C$_{8}$ & \centering\arraybackslash 61.0 \\
$~~~~~~~~~~~~~~~~~~$ C$_{10}$ & \centering\arraybackslash 58.2 \\
$~~~~~~~~~~~~~~~~~~$ C$_{12}$ & \centering\arraybackslash 58.8
\label{table:chainlengthG}
\end{longtable}
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% MIXED CHAINS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Effect of Mixed Chain Lengths}

Previous simulations have demonstrated non-monotonic behavior for $G$ as a function of the surface coverage.  One difficulty with the previous study was the ability of butanethiolate ligands to migrate on the Au(111) surface and to form segregated domains.  To simulate the effect of low coverages while preventing thiolate domain formation, we maintain 100\% thiolate coverage while varying the proportions of short (butanethiolate, C$_4$) and long (decanethiolate, C$_{10}$, or dodecanethiolate, C$_{12}$) alkyl chains.  Data on the conductance trend as the fraction of long chains was varied is shown in figure \ref{fig:Gstacks}.  Note that as in the previous study, $G$ is dependent upon solvent accessibility to thermally excited ligands.  Our simulations indicate a similar (but less dramatic) non-monotonic dependence on the fraction of long chains.
\begin{figure}
  \includegraphics[width=\linewidth]{figures/Gstacks}
  \caption{Interfacial thermal conductivity of mixed-chains has a non-monotonic dependence on the fraction of long chains (lower panels).  At low fractions of long chains, the solvent escape rate ($k_{escape}$) dominates the heat transfer process, while the solvent-thiolate orientational ordering ($<d>$) dominates in systems with higher fractions of long chains (upper panels).}
  \label{fig:Gstacks}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **DISCUSSION**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion} 

In the mixed chain-length simulations, solvent molecules
can become temporarily trapped or entangled with the thiolate chains. Their
residence in close proximity to the higher temperature environment close to
the surface allows them to carry heat away from the surface quite efficiently.
There are two aspects of this behavior that are relevant to thermal
conductance of the interface: the residence time of solvent molecules in the
thiolate layer, and the alignment of solvent molecules with the ligand alkyl
chains as a mechanism for transferring vibrational energy to these entrapped
solvent molecules. To quantify these competing effects, we have computed
solvent escape rates from the thiolate layer as well as a joint orientational
order parameter between the trapped solvent and the thiolate ligands.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% RESIDENCE TIME
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Mobility of solvent in the interfacial layer}

We use a simple survival correlation function, $C(t)$, to measure the
residence time of a solvent molecule in the  thiolate 
layer. This function correlates the identity of all hexane molecules
within the $z$-coordinate range of the thiolate layer at two separate
times.  If the solvent molecule is present at both times, the
configuration contributes a $1$, while the absence of the molecule at
the later time indicates that the solvent molecule has migrated into
the bulk, and this configuration contributes a $0$.  A steep decay in
$C(t)$ indicates a high turnover rate of solvent molecules from the
chain region to the bulk. We define the escape rate for trapped solvent molecules at the interface as 
\begin{equation}
 k_{escape} = \left( \int_0^T C(t) dt \right)^{-1}
  \label{eq:mobility}
\end{equation}
where T is the length of the simulation.  This is a direct measure of
the rate at which solvent molecules initially entangled in the thiolate layer
can escape into the bulk.  As $k_{escape} \rightarrow 0$, the
solvent becomes permanently trapped in the thiolate layer.  In
figure \ref{fig:Gstacks} we show that interfacial solvent mobility
decreases as the fraction of long thiolate chains increases.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% ORDER PARAMETER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

\subsection{Vibrational coupling via orientational ordering}

As the fraction of long-chain thiolates increases, the entrapped
solvent molecules must find specific orientations relative to the mean
orientation of the thiolate chains.  This alignment allows for
efficient thermal energy exchange between the thiolate alkyl chain and
the solvent molecules.

Once the interfacial solvent molecules have picked up thermal energy from the
thiolates, they carry heat away from the gold as they diffuse back
into the bulk solvent. When the percentage of long chains decreases,
the tails of the long chains are much more disordered and do not
provide structured channels for the solvent to fill.

To measure this cooperative ordering, we compute the orientational order
parameters and director axes for both the thiolate chains and for the
entrapped solvent. The director axis can be easily obtained by diagonalization
of the order parameter tensor,
\begin{equation} 
\mathsf{Q}_{\alpha \beta} = \frac{1}{2 N} \sum_{i=1}^{N} \left( 3 \mathbf{e}_{i
    \alpha} \mathbf{e}_{i \beta} - \delta_{\alpha \beta} \right)  
\end{equation}
where $\mathbf{e}_{i \alpha}$ is the $\alpha = x,y,z$ component of
the unit vector $\mathbf{e}_{i}$ along the long axis of molecule $i$.
For both the solvent and the ligand, the $\mathbf{e}$ vector is defined using
the terminal atoms of the molecule.

The largest eigenvalue of $\overleftrightarrow{\mathsf{Q}}$ is
traditionally used to obtain the orientational order parameter, while the
eigenvector corresponding to the order parameter yields the director
axis ($\mathbf{d}(t)$), which defines the average direction of
molecular alignment at any time $t$.  The overlap between the director
axes of the thiolates and the entrapped solvent is time-averaged,
        \begin{equation}
          \langle d \rangle = \langle \mathbf{d}_{thiolates} \left( t \right) \cdot
            \mathbf{d}_{solvent} \left( t \right) \rangle_t
          \label{eq:orientation}
        \end{equation}  
and reported in figure \ref{fig:Gstacks}. Values of $\langle d \rangle$ range from $0$ (solvent molecules in the ligand layer are perpendicular to the thiolate chains) to $1$ (solvent and ligand chains are aligned parallel to each other).

C$_4$ / C$_{10}$ mixed monolayers have a peak interfacial conductance with 75\% long chains. At this fraction of long chains, the cooperative orientational ordering of the solvent molecules and chains becomes the dominant effect while the solvent escape rate is quite slow. C$_4$ / C$_{12}$ mixtures have a peak interfacial conductance for 87.5\% long chains. The solvent-thiolate orientational ordering reaches its maximum value at this long chain fraction.  Long chain fractions of over $0.5$ for the C$_4$ / C$_{12}$ system are well ordered, but this effect is tempered by the exceptionally slow solvent escape rate ($\sim$ 1 molecule / 2 ns).
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                          **CONCLUSIONS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\section{Conclusions}
Our results suggest that a mixed vibrational transfer / convection
model may be necessary to explain the features of heat transfer at
this interface.  The alignment of the solvent chains with the ordered
ligand allows rapid transfer of energy to the trapped solvent and
becomes the dominant feature for ordered ligand layers.  Diffusion of
the vibrationally excited solvent into the bulk also plays a
significant role when the ligands are less tightly packed.

In the language of earlier continuum approaches to interfacial
conductance,\cite{RevModPhys.61.605} the alignment of the chains is an
important factor in the transfer of phonons from the thiolate layer to the
trapped solvent. The aligned solvent and thiolate chains have nearly identical
acoustic impedances and the phonons can scatter directly into a solvent
molecule that has been forced into alignment. When the entrapped solvent has
more configurations available, the likelihood of an impedance mismatch is
higher, and the phonon scatters into the solvent with  lower
efficiency. The fractional coverage of the long chains is therefore a simple
way of tuning the acoustic mismatch between the thiolate layer and the hexane
solvent.

Efficient heat transfer also can be accomplished via convective or diffusive motion of vibrationally excited solvent back into the bulk.  Once the entrapped solvent becomes too tightly aligned with the ligands, however, the convective avenue of heat transfer is cut off.  

Our simulations suggest a number of routes to make interfaces with high thermal conductance.  If it is possible to create an interface which forces the solvent into alignment with a ligand that shares many of the solvent's vibrational modes, while simultaneously preserving the ability of the solvent to diffuse back into the bulk, we would expect a significant jump in the interfacial conductance.  One possible way to accomplish this is to use polyene ligands with alternating unsaturated bonds. These are significantly more rigid than long alkanes, and could force solvent alignment (even at low relative coverages) while preserving mobility of solvent molecules within the ligand layer.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **ACKNOWLEDGMENTS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgments} 

We gratefully acknowledge conversations with Dr. Shenyu Kuang. Support for
this project was provided by the National Science Foundation under grant
CHE-0848243. Computational time was provided by the Center for Research
Computing (CRC) at the University of Notre Dame.

\newpage

\bibliography{thiolsRNEMD}

\newpage

\begin{figure}
 \centering{\includegraphics{figures/toc2}}
 \caption{Images of 75\%~C$_4$~/~25\%~C$_{12}$ (top panel) and 25\%~C$_4$~/~75\%~C$_{12}$ (bottom panel) interfaces at the beginning and end of 3 ns simulations. Solvent molecules that were initially present in the thiolate layer are colored light blue.  Diffusion of the initially-trapped solvent into the bulk is apparent in the interface with fewer long chains.  Trapped solvent is orientationally locked to the ordered ligands (and is less able to diffuse into the bulk) when the fraction of long chains increases.}
  \label{fig:toc}
\end{figure}

\end{doublespace}
\end{document}
                
                        
Revision:	3857
Committed:	Mon Feb 11 23:12:08 2013 UTC (12 years, 5 months ago) by kstocke1
Content type:	application/x-tex
File size:	41203 byte(s)
Log Message: