trunk/chainLength/GoldThiolsPaper.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{times}
\usepackage{mathptm}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
\usepackage{graphicx}
\usepackage{multirow}
\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}

\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}

\end{abstract}

\newpage

%\narrowtext

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

The structural and dynamical details of interfaces between metal
nanoparticles and solvents determines how energy flows between these
particles and their surroundings. Understanding this energy flow is
essential in designing and functionalizing metallic nanoparticles for
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
principle of this therapy is to pump the particles at high power at
the plasmon resonance and to allow heat dissipation to kill targeted
(e.g. cancerous) cells.  The relevant physical property controlling
this transfer of energy 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 highly 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 vital to get a detailed molecular
picture of particle-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) for $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. 

The thermal properties of various nanostructured interfaces have been
investigated experimentally by a number of groups: Cahill and
coworkers studied nanoscale thermal transport from metal
nanoparticle/fluid interfaces, to epitaxial TiN/single crystal oxides
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}
Wang {\it et al.} studied heat transport through long-chain
hydrocarbon monolayers on gold substrate at the individual molecular
level,\cite{Wang10082007} 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} and Juv\'e {\it
  et al.} studied the cooling dynamics, which is controlled by thermal
interface resistance of glass-embedded metal
nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
normally considered barriers for heat transport, Alper {\it et al.}
suggested that specific ligands (capping agents) could completely
eliminate this barrier
($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}

In previous simulations, we 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.  Our
calculations suggest 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.

One interesting result of our previous work was the observation of
{\it non-monotonic dependence} of the thermal conductance on the
coverage of a metal surface by a chemical protecting group.  Our
explanation for this behavior was that gaps in surface coverage
allowed solvent to penetrate close to the capping molecules that had
been heated by the metal surface, to absorb thermal energy from these
molecules, and then diffuse away.  The effect of surface coverage is
relatively difficult to study as the individual protecting groups have
lateral mobility on the surface and can aggregate to form domains on
the timescale of the simulation.

The work reported here involves the use of velocity shearing and
scaling reverse non-equilibrium molecular dynamics (VSS-RNEMD) to
study surfaces composed of mixed-length chains which collectively form
a complete monolayer on the surfaces.  These complete (but
mixed-chain) surfaces are significantly less prone to surface domain
formation, and would allow us to further investigate the mechanism of
heat transport to the solvent.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **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 by computing relevant time
correlation functions and 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,Kuang:2010uq} 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, and this has
undesirable side-effects when the applied flux becomes
large.\cite{Maginn:2010}        

Instead of having momentum exchange between {\it individual particles}
in each slab, the NIVS algorithm applies velocity scaling to all of
the selected particles in both slabs.\cite{Kuang:2010uq} A combination
of linear momentum, kinetic energy, and flux constraint equations
governs the amount of velocity scaling performed at each step.  NIVS
has been shown to be very effective at producing thermal gradients and
for computing thermal conductivities, particularly for heterogeneous
interfaces.  Although the NIVS algorithm can also be applied to impose
a directional momentum flux, thermal anisotropy was observed in
relatively high flux simulations, and the method is not suitable for
imposing a shear flux or for computing shear viscosities.

The most useful RNEMD
approach developed so far utilizes a series of simultaneous velocity
shearing and scaling exchanges between the two
slabs.\cite{2012MolPh.110..691K} 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, $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 total mass of particles within each slab 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 ({\it 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 of 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, and thermal anisotropy is
kept to a minimum.  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 {\it single 1 ns simulation}.\cite{2012MolPh.110..691K} 

\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.  We have been able to
demonstrate 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.}
  \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 gold/thiolate interface and the thiolate/solvent interface. We
could 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}
  \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. Force field parameters are needed for interactions both
between the same type of particles and between particles of different
species.

\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 metal lattice slab is 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 are used.  We have previously
utilized both united atom (UA) and all-atom (AA) force fields for
thermal conductivity in previous work,\cite{} 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 (200K) 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 butanethiol 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.  Table
\ref{table:pars} in the supporting information summarizes the
``metal/non-metal'' parameters utilized in our simulations.

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

We have implemented the VSS-RNEMD algorithm in OpenMD, our
group molecular dynamics code. A gold slab 11 atoms thick was
equilibrated at 1 atm and 200 K. The periodic box was expanded
in the z direction to expose two Au(111) faces.

A full monolayer of thiolates (1/3 the number of surface gold atoms) was placed on three-fold hollow sites on the gold interfaces. To efficiently test the effect of thiolate binding sites on the thermal conductance, all systems had one gold interface with thiolates placed only on fcc hollow sites and the other interface with thiolates only on hcp hollow sites. To test the effect of thiolate chain length on interfacial thermal conductance, full coverages of five chain lengths were tested: butanethiolate, hexanethiolate, octanethiolate, decanethiolate, and dodecanethiolate. To test the effect of mixed chain lengths, full coverages of butanethiolate/decanethiolate and butanethiolate/dodecanethiolate 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. Hexane, a straight chain flexible alkane, is very 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{Kuang2011} This overlap should provide a mechanism for thermal energy transfer from the thiolates to the solvent.

The system was equilibrated to 220 K in the NVT ensemble, allowing the thiolates and solvent to warm gradually. Pressure correction to 1 atm was done in an NPT ensemble that allowed expansion or contraction only in the z direction, so as not to disrupt the crystalline structure of the gold. The diagonal elements of the pressure tensor were monitored during the pressure correction step. If the xx and/or yy elements had a mean above zero throughout the simulation -- indicating residual pressure 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 NVT simulation was performed. All systems were equilibrated in the microcanonical (NVE) ensemble before proceeding with the VSS-RNEMD step.

A kinetic energy flux was applied using VSS-RNEMD in the NVE ensemble. The total simulation time was 3 nanoseconds, with velocity scaling 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 average temperature was 220 K, with a temperature difference between the hot and cold slabs of approximately 30 K. The average temperature and kinetic energy flux were carefully selected with two considerations in mind: 1) if the cold bin gets too cold (below ~180 K) the solvent may freeze or undergo a glassy transition, and 2) due to the deep sulfur-gold potential well, sulfur atoms routinely embed into the gold slab, particularly at temperatures above 250 K. Simulation conditions were chosen to avoid both of these undesirable situations. A reversed flux direction resulted in frozen long chain thiolates and solvent too near its boiling point.

A temperature profile of the system was 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}


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

We examined full coverages of five chain lengths, n = 4, 6, 8, 10, 12. 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. The trend of interfacial conductance is mostly flat 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. There is, however, a peak in conductance for a chain length of 6 (hexanethiolate). This may be due to the equivalent chain lengths of the hexane solvent and the alkyl chain of the capping agent, leading to an especially high degree of vibrational overlap between the thiolate and solvent. Strong vibrational overlap would allow for efficient thermal energy transfer across the thiolate/solvent interface.
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% MIXED CHAINS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Effect of Mixed Chain Lengths}

Previous work demonstrated that for butanethiolate monolayers on a Au(111) surface, the interfacial conductance was a non-monotonic function of the percent coverage. This is believed to be due to enhanced solvent-thiolate coupling through greater penetration of solvent molecules into the thiolate layer. At lower coverages, hexane solvent can more easily line up lengthwise with the thiolate tails by fitting into gaps between the thiolates. However, a side effect of low coverages is surface aggregation of the thiolates. To simulate the effect of low coverages while preventing aggregation we maintain 100\% thiolate coverage while varying the proportions of short (butanethiolate, n = 4) and long (decanethiolate, n = 10 or dodecanethiolate, n = 12). In systems where there is a mix of short and long chain thiolates, interfacial conductance is a non-monotonic function of the percent of long chains. The depth of the gaps between the long chains is $n_{long} - n_{short}$, which has implications for the ability of the hexane solvent to fill in the gaps between the long chains.
        
        \subsubsection{Butanethiolate/Decanethiolate}
Mixtures of butanethiolate/decanethiolate (n = 4, 10) have a peak interfacial condutance for 25\%/75\% short/long chains.
        
        \subsubsection{Butanethiolate/Dodecanethiolate} 
Mixtures of butanethiolate/dodecanethiolate (n = 4, 12) have a peak interfacial conductance for 12.5\%/87.5\% short/long chains. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **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 orientational
ordering of the C-C chains as a mechanism for transferring vibrational
energy to these entrapped solvent molecules.

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

We use a simple survival correlation function, $C(t)$, to quantify the
residence time of a solvent molecule in the long thiolate chain
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. % The correlation function is easily fit
% using a biexponential,
% \begin{equation}
%  C(t) = A \, e^{-t/\tau_{short}} + (1-A)  e^{-t/\tau_{long}}
%   \label{eq:biexponential}
% \end{equation}
% to determine short and long residence timescales and the relative populations of solvent molecules that can escape rapidly.  
We define the mobility of solvent molecules at the interface as 
\begin{equation}
 M = \int_{0}^{T} 1 - C(t) dt,
  \label{eq:mobility}
\end{equation}
where T is the length of the simulation.
In figure \ref{figure:res} we show that interfacial solvent mobility decreases as the percentage of long thiolate chains increases.

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

\subsection{Vibrational coupling via orientational ordering}

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

To quantify this cooperative
ordering, we computed 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}$ was the $\alpha = x,y,z$ component of
the unit vector $\mathbf{e}_{i}$ along the long axis of molecule $i$.
For both kinds of molecules, the $\mathbf{e}$ vector is defined using
the terminal atoms of the chains.

The largest eigenvalue of $\overleftrightarrow{\mathsf{Q}}$ is
traditionally used to obtain 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}
          \left \langle \mathbf{d}_{thiolates} \left( t \right) \cdot
            \mathbf{d}_{solvent} \left( t \right) \right \rangle, 
          \label{eq:orientation}
        \end{equation}  
and reported in table \ref{table:ordering}.

Once the 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.

Although the alignment of the chains with the entrapped solvent is one
possible mechanism for the non-monotonic increase in the conductance
as a function 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          **ACKNOWLEDGMENTS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgments}
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}

\end{doublespace}
\end{document}
                
                        
Revision:	3822
Committed:	Tue Dec 18 16:22:13 2012 UTC (12 years, 7 months ago) by gezelter
Content type:	application/x-tex
File size:	34462 byte(s)
Log Message:	A few changes