group/stokes/stokes.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{multirow}
%\usepackage{booktabs}
%\usepackage{bibentry}
%\usepackage{mathrsfs}
%\usepackage[ref]{overcite}
\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{\renewcommand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

%\renewcommand\citemid{\ } % no comma in optional reference note
\bibpunct{[}{]}{,}{n}{}{;}
\bibliographystyle{achemso}

\begin{document}

\title{ENTER TITLE HERE}

\author{Shenyu Kuang and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}
  REPLACE ABSTRACT HERE
  With the Non-Isotropic Velocity Scaling (NIVS) approach to Reverse
  Non-Equilibrium Molecular Dynamics (RNEMD) it is possible to impose
  an unphysical thermal flux between different regions of
  inhomogeneous systems such as solid / liquid interfaces.  We have
  applied NIVS to compute the interfacial thermal conductance at a
  metal / organic solvent interface that has been chemically capped by
  butanethiol molecules.  Our calculations suggest that coupling
  between the metal and liquid phases is enhanced by the capping
  agents, leading to a greatly enhanced conductivity at the interface.
  Specifically, the chemical bond between the metal and the capping
  agent introduces a vibrational overlap that is not present without
  the capping agent, and the overlap between the vibrational spectra
  (metal to cap, cap to solvent) provides a mechanism for rapid
  thermal transport across the interface. Our calculations also
  suggest that this is a non-monotonic function of the fractional
  coverage of the surface, as moderate coverages allow diffusive heat
  transport of solvent molecules that have been in close contact with
  the capping agent.

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
[DO THIS LATER]

[IMPORTANCE OF NANOSCALE TRANSPORT PROPERTIES STUDIES]

Due to the importance of heat flow (and heat removal) in
nanotechnology, interfacial thermal conductance has been studied
extensively both experimentally and computationally.\cite{cahill:793}
Nanoscale materials have a significant fraction of their atoms at
interfaces, and the chemical details of these interfaces govern the
thermal transport properties.  Furthermore, the interfaces are often
heterogeneous (e.g. solid - liquid), which provides a challenge to
computational methods which have been developed for homogeneous or
bulk systems.

Experimentally, the thermal properties of a number of interfaces have
been investigated.  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{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 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}

The acoustic mismatch model for interfacial conductance utilizes the
acoustic impedance ($Z_a = \rho_a v^s_a$) on both sides of the
interface.\cite{swartz1989} Here, $\rho_a$ and $v^s_a$ are the density
and speed of sound in material $a$.  The phonon transmission
probability at the $a-b$ interface is
\begin{equation}
t_{ab} = \frac{4 Z_a Z_b}{\left(Z_a + Z_b \right)^2},
\end{equation}
and the interfacial conductance can then be approximated as
\begin{equation}
G_{ab} \approx \frac{1}{4} C_D v_D t_{ab}
\end{equation}
where $C_D$ is the Debye heat capacity of the hot side, and $v_D$ is
the Debye phonon velocity ($1/v_D^3 = 1/3v_L^3 + 2/3 v_T^3$) where
$v_L$ and $v_T$ are the longitudinal and transverse speeds of sound,
respectively.  For the Au/hexane and Au/toluene interfaces, the
acoustic mismatch model predicts room-temperature $G \approx 87 \mbox{
  and } 129$ MW m$^{-2}$ K$^{-1}$, respectively.  However, it is not
clear how to apply the acoustic mismatch model to a
chemically-modified surface, particularly when the acoustic properties
of a monolayer film may not be well characterized.

[PREVIOUS METHODS INCLUDING NIVS AND THEIR LIMITATIONS]
[DIFFICULTY TO GENERATE JZKE AND JZP SIMUTANEOUSLY]

More precise computational models have also been used to study the
interfacial thermal transport in order to gain an understanding of
this phenomena at the molecular level. 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} 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. The comparatively low thermal flux through interfaces is
difficult to measure with Equilibrium
MD\cite{doi:10.1080/0026897031000068578} or forward NEMD simulation
methods. Therefore, the Reverse NEMD (RNEMD)
methods\cite{MullerPlathe:1997xw,kuang:164101} would be advantageous
in that they {\it apply} the difficult to measure quantity (flux),
while {\it measuring} the easily-computed quantity (the thermal
gradient).  This is particularly true for inhomogeneous interfaces
where it would not be clear how to apply a gradient {\it a priori}.
Garde and coworkers\cite{garde:nl2005,garde:PhysRevLett2009} applied
this approach to various liquid interfaces and studied how thermal
conductance (or resistance) is dependent on chemical details of a
number of hydrophobic and hydrophilic aqueous interfaces. And
recently, Luo {\it et al.} studied the thermal conductance of
Au-SAM-Au junctions using the same approach, comparing to a constant
temperature difference method.\cite{Luo20101} While this latter
approach establishes more ideal Maxwell-Boltzmann distributions than
previous RNEMD methods, it does not guarantee momentum or kinetic
energy conservation.

Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS)
algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
retains the desirable features of RNEMD (conservation of linear
momentum and total energy, compatibility with periodic boundary
conditions) while establishing true thermal distributions in each of
the two slabs. Furthermore, it allows effective thermal exchange
between particles of different identities, and thus makes the study of
interfacial conductance much simpler.

[WHAT IS COVERED IN THIS MANUSCRIPT]
[MAY PUT FIGURE 1 HERE]
The work presented here deals with the Au(111) surface covered to
varying degrees by butanethiol, a capping agent with short carbon
chain, and solvated with organic solvents of different molecular
properties. Different models were used for both the capping agent and
the solvent force field parameters. Using the NIVS algorithm, the
thermal transport across these interfaces was studied and the
underlying mechanism for the phenomena was investigated.

\section{Methodology}
Similar to the NIVS methodology,\cite{kuang:164101} we consider a
periodic system divided into a series of slabs along a certain axis
(e.g. $z$). The unphysical thermal and/or momentum flux is designated
from the center slab to one of the end slabs, and thus the center slab
would have a lower temperature than the end slab (unless the thermal
flux is negative). Therefore, the center slab is denoted as ``$c$''
while the end slab as ``$h$''.

To impose these fluxes, we periodically apply separate operations to
velocities of particles {$i$} within the center slab and of particles
{$j$} within the end slab:
\begin{eqnarray}
\vec{v}_i & \leftarrow & c\cdot\left(\vec{v}_i - \langle\vec{v}_c
  \rangle\right) + \left(\langle\vec{v}_c\rangle + \vec{a}_c\right) \\
\vec{v}_j & \leftarrow & h\cdot\left(\vec{v}_j - \langle\vec{v}_h
  \rangle\right) + \left(\langle\vec{v}_h\rangle + \vec{a}_h\right)
\end{eqnarray}
where $\langle\vec{v}_c\rangle$ and $\langle\vec{v}_h\rangle$ denotes
the instantaneous bulk velocity of slabs $c$ and $h$ respectively
before an operation occurs. When a momentum flux $\vec{j}_z(\vec{p})$
presents, these bulk velocities would have a corresponding change
($\vec{a}_c$ and $\vec{a}_h$ respectively) according to Newton's
second law:
\begin{eqnarray}
M_c \vec{a}_c & = & -\vec{j}_z(\vec{p}) \Delta t \\
M_h \vec{a}_h & = &  \vec{j}_z(\vec{p}) \Delta t
\end{eqnarray}
where
\begin{eqnarray}
M_c & = & \sum_{i = 1}^{N_c} m_i \\
M_h & = & \sum_{j = 1}^{N_h} m_j
\end{eqnarray}
and $\Delta t$ is the interval between two operations.

The above operations conserve the linear momentum of a periodic
system. To satisfy total energy conservation as well as to impose a
thermal flux $J_z$, one would have
%SUPPORT INFO MIGHT BE NECESSARY TO PUT EXTRA MATH IN
\begin{eqnarray}
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c
\rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\
K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\vec{v}_h
\rangle^2) + \frac{1}{2}M_h (\langle \vec{v}_h \rangle + \vec{a}_h)^2
\end{eqnarray}
where $K_c$ and $K_h$ denotes translational kinetic energy of slabs
$c$ and $h$ respectively before an operation occurs. These
translational kinetic energy conservation equations are sufficient to
ensure total energy conservation, as the operations applied do not
change the potential energy of a system, given that the potential
energy does not depend on particle velocity.

The above sets of equations are sufficient to determine the velocity
scaling coefficients ($c$ and $h$) as well as $\vec{a}_c$ and
$\vec{a}_h$. Note that two roots of $c$ and $h$ exist
respectively. However, to avoid dramatic perturbations to a system,
the positive roots (which are closer to 1) are chosen.

By implementing these operations at a certain frequency, a steady
thermal and/or momentum flux can be applied and the corresponding
temperature and/or momentum gradients can be established.
[REFER TO NIVS PAPER]
[ADVANTAGES]

Steady state MD simulations have an advantage in that not many
trajectories are needed to study the relationship between thermal flux
and thermal gradients. For systems with low interfacial conductance,
one must have a method capable of generating or measuring relatively
small fluxes, compared to those required for bulk conductivity. This
requirement makes the calculation even more difficult for
slowly-converging equilibrium methods.\cite{Viscardy:2007lq} Forward
NEMD methods impose a gradient (and measure a flux), but at interfaces
it is not clear what behavior should be imposed at the boundaries
between materials.  Imposed-flux reverse non-equilibrium
methods\cite{MullerPlathe:1997xw} set the flux {\it a priori} and
the thermal response becomes an easy-to-measure quantity.  Although
M\"{u}ller-Plathe's original momentum swapping approach can be used
for exchanging energy between particles of different identity, the
kinetic energy transfer efficiency is affected by the mass difference
between the particles, which limits its application on heterogeneous
interfacial systems.

The non-isotropic velocity scaling (NIVS) \cite{kuang:164101} approach
to non-equilibrium MD simulations is able to impose a wide range of
kinetic energy fluxes without obvious perturbation to the velocity
distributions of the simulated systems. Furthermore, this approach has
the advantage in heterogeneous interfaces in that kinetic energy flux
can be applied between regions of particles of arbitrary identity, and
the flux will not be restricted by difference in particle mass.

The NIVS algorithm scales the velocity vectors in two separate regions
of a simulation system with respective diagonal scaling matrices. To
determine these scaling factors in the matrices, a set of equations
including linear momentum conservation and kinetic energy conservation
constraints and target energy flux satisfaction is solved. With the
scaling operation applied to the system in a set frequency, bulk
temperature gradients can be easily established, and these can be used
for computing thermal conductivities. The NIVS algorithm conserves
momenta and energy and does not depend on an external thermostat.

\subsection{Defining Interfacial Thermal Conductivity ($G$)}

For an interface with relatively low interfacial conductance, and a
thermal flux between two distinct bulk regions, the regions on either
side of the interface rapidly come to a state in which the two phases
have relatively homogeneous (but distinct) temperatures. The
interfacial thermal conductivity $G$ can therefore be approximated as:
\begin{equation}
  G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
    \langle T_\mathrm{cold}\rangle \right)}
\label{lowG}
\end{equation}
where ${E_{total}}$ is the total imposed non-physical kinetic energy
transfer during the simulation and ${\langle T_\mathrm{hot}\rangle}$
and ${\langle T_\mathrm{cold}\rangle}$ are the average observed
temperature of the two separated phases.  For an applied flux $J_z$
operating over a simulation time $t$ on a periodically-replicated slab
of dimensions $L_x \times L_y$, $E_{total} = 2 J_z t L_x L_y$.

When the interfacial conductance is {\it not} small, there are two
ways to define $G$. One common way is to assume the temperature is
discrete on the two sides of the interface. $G$ can be calculated
using the applied thermal flux $J$ and the maximum temperature
difference measured along the thermal gradient max($\Delta T$), which
occurs at the Gibbs dividing surface (Figure \ref{demoPic}). This is
known as the Kapitza conductance, which is the inverse of the Kapitza
resistance.
\begin{equation}
  G=\frac{J}{\Delta T}
\label{discreteG}
\end{equation}

\begin{figure}
\includegraphics[width=\linewidth]{method}
\caption{Interfacial conductance can be calculated by applying an
  (unphysical) kinetic energy flux between two slabs, one located
  within the metal and another on the edge of the periodic box.  The
  system responds by forming a thermal gradient.  In bulk liquids,
  this gradient typically has a single slope, but in interfacial
  systems, there are distinct thermal conductivity domains.  The
  interfacial conductance, $G$ is found by measuring the temperature
  gap at the Gibbs dividing surface, or by using second derivatives of
  the thermal profile.}
\label{demoPic}
\end{figure}

Another approach is to assume that the temperature is continuous and
differentiable throughout the space. Given that $\lambda$ is also
differentiable, $G$ can be defined as its gradient ($\nabla\lambda$)
projected along a vector normal to the interface ($\mathbf{\hat{n}}$)
and evaluated at the interface location ($z_0$). This quantity,
\begin{align}
G^\prime &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\
         &= \frac{\partial}{\partial z}\left(-\frac{J_z}{
           \left(\frac{\partial T}{\partial z}\right)}\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}
has the same units as the common definition for $G$, and the maximum
of its magnitude denotes where thermal conductivity has the largest
change, i.e. the interface.  In the geometry used in this study, the
vector normal to the interface points along the $z$ axis, as do
$\vec{J}$ and the thermal gradient.  This yields the simplified
expressions in Eq. \ref{derivativeG}. 

With temperature profiles obtained from simulation, one is able to
approximate the first and second derivatives of $T$ with finite
difference methods and calculate $G^\prime$. In what follows, both
definitions have been used, and are compared in the results.

To investigate the interfacial conductivity at metal / solvent
interfaces, we have modeled a metal slab with its (111) surfaces
perpendicular to the $z$-axis of our simulation cells. The metal slab
has been prepared both with and without capping agents on the exposed
surface, and has been solvated with simple organic solvents, as
illustrated in Figure \ref{gradT}.

With the simulation cell described above, we are able to equilibrate
the system and impose an unphysical thermal flux between the liquid
and the metal phase using the NIVS algorithm. By periodically applying
the unphysical flux, we obtained a temperature profile and its spatial
derivatives. Figure \ref{gradT} shows how an applied thermal flux can
be used to obtain the 1st and 2nd derivatives of the temperature
profile.

\begin{figure}
\includegraphics[width=\linewidth]{gradT}
\caption{A sample of Au (111) / butanethiol / hexane interfacial
  system with the temperature profile after a kinetic energy flux has
  been imposed.  Note that the largest temperature jump in the thermal
  profile (corresponding to the lowest interfacial conductance) is at
  the interface between the butanethiol molecules (blue) and the
  solvent (grey).  First and second derivatives of the temperature
  profile are obtained using a finite difference approximation (lower
  panel).}
\label{gradT}
\end{figure}

\section{Computational Details}
\subsection{Simulation Protocol}
The NIVS algorithm has been implemented in our MD simulation code,
OpenMD\cite{Meineke:2005gd,openmd}, and was used for our simulations.
Metal slabs of 6 or 11 layers of Au atoms were first equilibrated
under atmospheric pressure (1 atm) and 200K. After equilibration,
butanethiol capping agents were placed at three-fold hollow sites on
the Au(111) surfaces. These sites are either {\it fcc} or {\it
  hcp} sites, although Hase {\it et al.} found that they are
equivalent in a heat transfer process,\cite{hase:2010} so we did not
distinguish between these sites in our study. The maximum butanethiol
capacity on Au surface is $1/3$ of the total number of surface Au
atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
series of lower coverages was also prepared by eliminating
butanethiols from the higher coverage surface in a regular manner. The
lower coverages were prepared in order to study the relation between
coverage and interfacial conductance.

The capping agent molecules were allowed to migrate during the
simulations. They distributed themselves uniformly and sampled a
number of three-fold sites throughout out study. Therefore, the
initial configuration does not noticeably affect the sampling of a
variety of configurations of the same coverage, and the final
conductance measurement would be an average effect of these
configurations explored in the simulations.

After the modified Au-butanethiol surface systems were equilibrated in
the canonical (NVT) ensemble, organic solvent molecules were packed in
the previously empty part of the simulation cells.\cite{packmol} Two
solvents were investigated, one which has little vibrational overlap
with the alkanethiol and which has a planar shape (toluene), and one
which has similar vibrational frequencies to the capping agent and
chain-like shape ({\it n}-hexane).

The simulation cells were not particularly extensive along the
$z$-axis, as a very long length scale for the thermal gradient may
cause excessively hot or cold temperatures in the middle of the
solvent region and lead to undesired phenomena such as solvent boiling
or freezing when a thermal flux is applied. Conversely, too few
solvent molecules would change the normal behavior of the liquid
phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
these extreme cases did not happen to our simulations. The spacing
between periodic images of the gold interfaces is $45 \sim 75$\AA in
our simulations.

The initial configurations generated are further equilibrated with the
$x$ and $y$ dimensions fixed, only allowing the $z$-length scale to
change. This is to ensure that the equilibration of liquid phase does
not affect the metal's crystalline structure. Comparisons were made
with simulations that allowed changes of $L_x$ and $L_y$ during NPT
equilibration. No substantial changes in the box geometry were noticed
in these simulations. After ensuring the liquid phase reaches
equilibrium at atmospheric pressure (1 atm), further equilibration was
carried out under canonical (NVT) and microcanonical (NVE) ensembles.

After the systems reach equilibrium, NIVS was used to impose an
unphysical thermal flux between the metal and the liquid phases. Most
of our simulations were done under an average temperature of
$\sim$200K. Therefore, thermal flux usually came from the metal to the
liquid so that the liquid has a higher temperature and would not
freeze due to lowered temperatures. After this induced temperature
gradient had stabilized, the temperature profile of the simulation cell
was recorded. To do this, the simulation cell is divided evenly into
$N$ slabs along the $z$-axis. The average temperatures of each slab
are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
the same, the derivatives of $T$ with respect to slab number $n$ can
be directly used for $G^\prime$ calculations: \begin{equation}
  G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
         \Big/\left(\frac{\partial T}{\partial z}\right)^2
         = |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
         \Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
         = |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big|
         \Big/\left(\frac{\partial T}{\partial n}\right)^2
\label{derivativeG2}
\end{equation}
The absolute values in Eq. \ref{derivativeG2} appear because the
direction of the flux $\vec{J}$ is in an opposing direction on either
side of the metal slab. 

All of the above simulation procedures use a time step of 1 fs. Each
equilibration stage took a minimum of 100 ps, although in some cases,
longer equilibration stages were utilized.

\subsection{Force Field Parameters}
Our simulations include a number of chemically distinct components.
Figure \ref{demoMol} demonstrates the sites defined for both
United-Atom and All-Atom models of the organic solvent and capping
agents 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]{structures}
\caption{Structures of the capping agent and solvents utilized in
  these simulations. The chemically-distinct sites (a-e) are expanded
  in terms of constituent atoms for both United Atom (UA) and All Atom
  (AA) force fields.  Most parameters are from References
  \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes,TraPPE-UA.thiols}
  (UA) and \protect\cite{OPLSAA} (AA). Cross-interactions with the Au
  atoms are given in Table 1 in the supporting information.}
\label{demoMol}
\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 two solvent molecules, {\it n}-hexane and toluene, two
different atomistic models were utilized. Both solvents were modeled
using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
for our UA solvent molecules. In these models, 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.

By eliminating explicit hydrogen atoms, the TraPPE-UA models are
simple and computationally efficient, while maintaining good accuracy.
However, 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. If heat is transferred to the liquid phase during the
NIVS simulation, the liquid in the hot slab can actually be
substantially warmer than the mean temperature in the simulation. The
lower mean temperatures therefore prevent solvent boiling.

For UA-toluene, the non-bonded potentials between intermolecular sites
have a similar Lennard-Jones formulation. The toluene molecules were
treated as a single rigid body, so there was no need for
intramolecular interactions (including bonds, bends, or torsions) in
this solvent model.

Besides the TraPPE-UA models, AA models for both organic solvents are
included in our studies as well. The OPLS-AA\cite{OPLSAA} force fields
were used. For hexane, additional explicit hydrogen sites were
included. Besides bonding and non-bonded site-site interactions,
partial charges and the electrostatic interactions were added to each
CT and HC site. For toluene, a flexible model for the toluene molecule
was utilized which included bond, bend, torsion, and inversion
potentials to enforce ring planarity. 

The butanethiol capping agent in our simulations, were also modeled
with both UA and AA model. The TraPPE-UA force field includes
parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
UA butanethiol model in our simulations. The OPLS-AA also provides
parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
surfaces do not have the hydrogen atom bonded to sulfur. To derive
suitable parameters for butanethiol adsorbed on Au(111) surfaces, we
adopt the S parameters from Luedtke and Landman\cite{landman:1998} and
modify the parameters for the CTS atom to maintain charge neutrality
in the molecule.  Note that the model choice (UA or AA) for the capping
agent can be different from the solvent. Regardless of model choice,
the force field parameters for interactions between capping agent and
solvent can be derived using Lorentz-Berthelot Mixing Rule:
\begin{eqnarray}
  \sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
  \epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}} 
\end{eqnarray}

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.

The potentials developed from {\it ab initio} calculations by Leng
{\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
interactions between Au and aromatic C/H atoms in toluene. However,
the Lennard-Jones parameters between Au and other types of particles,
(e.g. AA alkanes) have not yet been established. For these
interactions, the Lorentz-Berthelot mixing rule can be used to derive
effective single-atom LJ parameters for the metal using the fit values
for toluene. These are then used to construct reasonable mixing
parameters for the interactions between the gold and other atoms.
Table 1 in the supporting information summarizes the
``metal/non-metal'' parameters utilized in our simulations.

\section{Results}
[L-J COMPARED TO RENMD NIVS; WATER COMPARED TO RNEMD NIVS;
SLIP BOUNDARY VS STICK BOUNDARY; ICE-WATER INTERFACES]

There are many factors contributing to the measured interfacial
conductance; some of these factors are physically motivated
(e.g. coverage of the surface by the capping agent coverage and
solvent identity), while some are governed by parameters of the
methodology (e.g. applied flux and the formulas used to obtain the
conductance). In this section we discuss the major physical and
calculational effects on the computed conductivity.

\subsection{Effects due to capping agent coverage}

A series of different initial conditions with a range of surface
coverages was prepared and solvated with various with both of the
solvent molecules. These systems were then equilibrated and their
interfacial thermal conductivity was measured with the NIVS
algorithm. Figure \ref{coverage} demonstrates the trend of conductance
with respect to surface coverage. 

\begin{figure}
\includegraphics[width=\linewidth]{coverage}
\caption{The interfacial thermal conductivity ($G$) has a
  non-monotonic dependence on the degree of surface capping.  This
  data is for the Au(111) / butanethiol / solvent interface with
  various UA force fields at $\langle T\rangle \sim $200K.}
\label{coverage}
\end{figure}

In partially covered surfaces, the derivative definition for
$G^\prime$ (Eq. \ref{derivativeG}) becomes difficult to apply, as the
location of maximum change of $\lambda$ becomes washed out.  The
discrete definition (Eq. \ref{discreteG}) is easier to apply, as the
Gibbs dividing surface is still well-defined. Therefore, $G$ (not
$G^\prime$) was used in this section.

From Figure \ref{coverage}, one can see the significance of the
presence of capping agents. When even a small fraction of the Au(111)
surface sites are covered with butanethiols, the conductivity exhibits
an enhancement by at least a factor of 3.  Capping agents are clearly
playing a major role in thermal transport at metal / organic solvent
surfaces.

We note a non-monotonic behavior in the interfacial conductance as a
function of surface coverage. The maximum conductance (largest $G$)
happens when the surfaces are about 75\% covered with butanethiol
caps.  The reason for this behavior is not entirely clear.  One
explanation is that incomplete butanethiol coverage allows small gaps
between butanethiols to form. These gaps can be filled by transient
solvent molecules.  These solvent molecules couple very strongly with
the hot capping agent molecules near the surface, and can then carry
away (diffusively) the excess thermal energy from the surface.

There appears to be a competition between the conduction of the
thermal energy away from the surface by the capping agents (enhanced
by greater coverage) and the coupling of the capping agents with the
solvent (enhanced by interdigitation at lower coverages).  This
competition would lead to the non-monotonic coverage behavior observed
here.

Results for rigid body toluene solvent, as well as the UA hexane, are
within the ranges expected from prior experimental
work.\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406} This suggests
that explicit hydrogen atoms might not be required for modeling
thermal transport in these systems.  C-H vibrational modes do not see
significant excited state population at low temperatures, and are not
likely to carry lower frequency excitations from the solid layer into
the bulk liquid.

The toluene solvent does not exhibit the same behavior as hexane in
that $G$ remains at approximately the same magnitude when the capping
coverage increases from 25\% to 75\%.  Toluene, as a rigid planar
molecule, cannot occupy the relatively small gaps between the capping
agents as easily as the chain-like {\it n}-hexane.  The effect of
solvent coupling to the capping agent is therefore weaker in toluene
except at the very lowest coverage levels.  This effect counters the
coverage-dependent conduction of heat away from the metal surface,
leading to a much flatter $G$ vs. coverage trend than is observed in
{\it n}-hexane.

\subsection{Effects due to Solvent \& Solvent Models}
In addition to UA solvent and capping agent models, AA models have
also been included in our simulations.  In most of this work, the same
(UA or AA) model for solvent and capping agent was used, but it is
also possible to utilize different models for different components.
We have also included isotopic substitutions (Hydrogen to Deuterium)
to decrease the explicit vibrational overlap between solvent and
capping agent. Table \ref{modelTest} summarizes the results of these
studies.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed interfacial thermal conductance ($G$ and
        $G^\prime$) values for interfaces using various models for
        solvent and capping agent (or without capping agent) at
        $\langle T\rangle\sim$200K.  Here ``D'' stands for deuterated
        solvent or capping agent molecules. Error estimates are
        indicated in parentheses.}
      
      \begin{tabular}{llccc}
        \hline\hline
        Butanethiol model & Solvent & $G$ & $G^\prime$ \\
        (or bare surface) & model & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        UA    & UA hexane    & 131(9)    & 87(10)    \\
              & UA hexane(D) & 153(5)    & 136(13)   \\
              & AA hexane    & 131(6)    & 122(10)   \\
              & UA toluene   & 187(16)   & 151(11)   \\
              & AA toluene   & 200(36)   & 149(53)   \\
        \hline
        AA    & UA hexane    & 116(9)    & 129(8)    \\
              & AA hexane    & 442(14)   & 356(31)   \\
              & AA hexane(D) & 222(12)   & 234(54)   \\
              & UA toluene   & 125(25)   & 97(60)    \\
              & AA toluene   & 487(56)   & 290(42)   \\
        \hline
        AA(D) & UA hexane    & 158(25)   & 172(4)    \\
              & AA hexane    & 243(29)   & 191(11)   \\
              & AA toluene   & 364(36)   & 322(67)   \\
        \hline
        bare  & UA hexane    & 46.5(3.2) & 49.4(4.5) \\
              & UA hexane(D) & 43.9(4.6) & 43.0(2.0) \\
              & AA hexane    & 31.0(1.4) & 29.4(1.3) \\
              & UA toluene   & 70.1(1.3) & 65.8(0.5) \\
        \hline\hline
      \end{tabular}
      \label{modelTest}
    \end{center}
  \end{minipage}
\end{table*}

To facilitate direct comparison between force fields, systems with the
same capping agent and solvent were prepared with the same length
scales for the simulation cells.

On bare metal / solvent surfaces, different force field models for
hexane yield similar results for both $G$ and $G^\prime$, and these
two definitions agree with each other very well. This is primarily an
indicator of weak interactions between the metal and the solvent.

For the fully-covered surfaces, the choice of force field for the
capping agent and solvent has a large impact on the calculated values
of $G$ and $G^\prime$.  The AA thiol to AA solvent conductivities are
much larger than their UA to UA counterparts, and these values exceed
the experimental estimates by a large measure.  The AA force field
allows significant energy to go into C-H (or C-D) stretching modes,
and since these modes are high frequency, this non-quantum behavior is
likely responsible for the overestimate of the conductivity.  Compared
to the AA model, the UA model yields more reasonable conductivity
values with much higher computational efficiency.

\subsubsection{Are electronic excitations in the metal important?}
Because they lack electronic excitations, the QSC and related embedded
atom method (EAM) models for gold are known to predict unreasonably
low values for bulk conductivity
($\lambda$).\cite{kuang:164101,ISI:000207079300006,Clancy:1992} If the
conductance between the phases ($G$) is governed primarily by phonon
excitation (and not electronic degrees of freedom), one would expect a
classical model to capture most of the interfacial thermal
conductance.  Our results for $G$ and $G^\prime$ indicate that this is
indeed the case, and suggest that the modeling of interfacial thermal
transport depends primarily on the description of the interactions
between the various components at the interface.  When the metal is
chemically capped, the primary barrier to thermal conductivity appears
to be the interface between the capping agent and the surrounding
solvent, so the excitations in the metal have little impact on the
value of $G$.

\subsection{Effects due to methodology and simulation parameters}

We have varied the parameters of the simulations in order to
investigate how these factors would affect the computation of $G$.  Of
particular interest are: 1) the length scale for the applied thermal
gradient (modified by increasing the amount of solvent in the system),
2) the sign and magnitude of the applied thermal flux, 3) the average
temperature of the simulation (which alters the solvent density during
equilibration), and 4) the definition of the interfacial conductance
(Eqs. (\ref{discreteG}) or (\ref{derivativeG})) used in the
calculation.

Systems of different lengths were prepared by altering the number of
solvent molecules and extending the length of the box along the $z$
axis to accomodate the extra solvent.  Equilibration at the same
temperature and pressure conditions led to nearly identical surface
areas ($L_x$ and $L_y$) available to the metal and capping agent,
while the extra solvent served mainly to lengthen the axis that was
used to apply the thermal flux.  For a given value of the applied
flux, the different $z$ length scale has only a weak effect on the
computed conductivities.

\subsubsection{Effects of applied flux}
The NIVS algorithm allows changes in both the sign and magnitude of
the applied flux.  It is possible to reverse the direction of heat
flow simply by changing the sign of the flux, and thermal gradients
which would be difficult to obtain experimentally ($5$ K/\AA) can be
easily simulated.  However, the magnitude of the applied flux is not
arbitrary if one aims to obtain a stable and reliable thermal gradient.
A temperature gradient can be lost in the noise if $|J_z|$ is too
small, and excessive $|J_z|$ values can cause phase transitions if the
extremes of the simulation cell become widely separated in
temperature.  Also, if $|J_z|$ is too large for the bulk conductivity
of the materials, the thermal gradient will never reach a stable
state.  

Within a reasonable range of $J_z$ values, we were able to study how
$G$ changes as a function of this flux.  In what follows, we use
positive $J_z$ values to denote the case where energy is being
transferred by the method from the metal phase and into the liquid.
The resulting gradient therefore has a higher temperature in the
liquid phase.  Negative flux values reverse this transfer, and result
in higher temperature metal phases.  The conductance measured under
different applied $J_z$ values is listed in Tables 2 and 3 in the
supporting information. These results do not indicate that $G$ depends
strongly on $J_z$ within this flux range. The linear response of flux
to thermal gradient simplifies our investigations in that we can rely
on $G$ measurement with only a small number $J_z$ values.

The sign of $J_z$ is a different matter, however, as this can alter
the temperature on the two sides of the interface. The average
temperature values reported are for the entire system, and not for the
liquid phase, so at a given $\langle T \rangle$, the system with
positive $J_z$ has a warmer liquid phase.  This means that if the
liquid carries thermal energy via diffusive transport, {\it positive}
$J_z$ values will result in increased molecular motion on the liquid
side of the interface, and this will increase the measured
conductivity.

\subsubsection{Effects due to average temperature}

We also studied the effect of average system temperature on the
interfacial conductance.  The simulations are first equilibrated in
the NPT ensemble at 1 atm.  The TraPPE-UA model for hexane tends to
predict a lower boiling point (and liquid state density) than
experiments.  This lower-density liquid phase leads to reduced contact
between the hexane and butanethiol, and this accounts for our
observation of lower conductance at higher temperatures.  In raising
the average temperature from 200K to 250K, the density drop of
$\sim$20\% in the solvent phase leads to a $\sim$40\% drop in the
conductance.

Similar behavior is observed in the TraPPE-UA model for toluene,
although this model has better agreement with the experimental
densities of toluene.  The expansion of the toluene liquid phase is
not as significant as that of the hexane (8.3\% over 100K), and this
limits the effect to $\sim$20\% drop in thermal conductivity.

Although we have not mapped out the behavior at a large number of
temperatures, is clear that there will be a strong temperature
dependence in the interfacial conductance when the physical properties
of one side of the interface (notably the density) change rapidly as a
function of temperature.

Besides the lower interfacial thermal conductance, surfaces at
relatively high temperatures are susceptible to reconstructions,
particularly when butanethiols fully cover the Au(111) surface. These
reconstructions include surface Au atoms which migrate outward to the
S atom layer, and butanethiol molecules which embed into the surface
Au layer. The driving force for this behavior is the strong Au-S
interactions which are modeled here with a deep Lennard-Jones
potential. This phenomenon agrees with reconstructions that have been
experimentally
observed.\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}.  Vlugt
{\it et al.} kept their Au(111) slab rigid so that their simulations
could reach 300K without surface
reconstructions.\cite{vlugt:cpc2007154} Since surface reconstructions
blur the interface, the measurement of $G$ becomes more difficult to
conduct at higher temperatures.  For this reason, most of our
measurements are undertaken at $\langle T\rangle\sim$200K where
reconstruction is minimized.

However, when the surface is not completely covered by butanethiols,
the simulated system appears to be more resistent to the
reconstruction. Our Au / butanethiol / toluene system had the Au(111)
surfaces 90\% covered by butanethiols, but did not see this above
phenomena even at $\langle T\rangle\sim$300K.  That said, we did
observe butanethiols migrating to neighboring three-fold sites during
a simulation.  Since the interface persisted in these simulations, we
were able to obtain $G$'s for these interfaces even at a relatively
high temperature without being affected by surface reconstructions.

\section{Discussion}
[COMBINE W. RESULTS]
The primary result of this work is that the capping agent acts as an
efficient thermal coupler between solid and solvent phases.  One of
the ways the capping agent can carry out this role is to down-shift
between the phonon vibrations in the solid (which carry the heat from
the gold) and the molecular vibrations in the liquid (which carry some
of the heat in the solvent).

To investigate the mechanism of interfacial thermal conductance, the
vibrational power spectrum was computed. Power spectra were taken for
individual components in different simulations. To obtain these
spectra, simulations were run after equilibration in the
microcanonical (NVE) ensemble and without a thermal
gradient. Snapshots of configurations were collected at a frequency
that is higher than that of the fastest vibrations occurring in the
simulations. With these configurations, the velocity auto-correlation
functions can be computed:
\begin{equation}
C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
\label{vCorr} 
\end{equation}
The power spectrum is constructed via a Fourier transform of the
symmetrized velocity autocorrelation function,
\begin{equation}
  \hat{f}(\omega) =
  \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
\label{fourier} 
\end{equation}

\subsection{The role of specific vibrations}
The vibrational spectra for gold slabs in different environments are
shown as in Figure \ref{specAu}. Regardless of the presence of
solvent, the gold surfaces which are covered by butanethiol molecules
exhibit an additional peak observed at a frequency of
$\sim$165cm$^{-1}$.  We attribute this peak to the S-Au bonding
vibration. This vibration enables efficient thermal coupling of the
surface Au layer to the capping agents. Therefore, in our simulations,
the Au / S interfaces do not appear to be the primary barrier to
thermal transport when compared with the butanethiol / solvent
interfaces.  This supports the results of Luo {\it et
  al.}\cite{Luo20101}, who reported $G$ for Au-SAM junctions roughly
twice as large as what we have computed for the thiol-liquid
interfaces.

\begin{figure}
\includegraphics[width=\linewidth]{vibration}
\caption{The vibrational power spectrum for thiol-capped gold has an
  additional vibrational peak at $\sim $165cm$^{-1}$.  Bare gold
  surfaces (both with and without a solvent over-layer) are missing
  this peak.   A similar peak at  $\sim $165cm$^{-1}$ also appears in
  the vibrational power spectrum for the butanethiol capping agents.}
\label{specAu}
\end{figure}

Also in this figure, we show the vibrational power spectrum for the
bound butanethiol molecules, which also exhibits the same
$\sim$165cm$^{-1}$ peak.

\subsection{Overlap of power spectra}
A comparison of the results obtained from the two different organic
solvents can also provide useful information of the interfacial
thermal transport process.  In particular, the vibrational overlap
between the butanethiol and the organic solvents suggests a highly
efficient thermal exchange between these components.  Very high
thermal conductivity was observed when AA models were used and C-H
vibrations were treated classically. The presence of extra degrees of
freedom in the AA force field yields higher heat exchange rates
between the two phases and results in a much higher conductivity than
in the UA force field. The all-atom classical models include high
frequency modes which should be unpopulated at our relatively low
temperatures.  This artifact is likely the cause of the high thermal
conductance in all-atom MD simulations.

The similarity in the vibrational modes available to solvent and
capping agent can be reduced by deuterating one of the two components
(Fig. \ref{aahxntln}).  Once either the hexanes or the butanethiols
are deuterated, one can observe a significantly lower $G$ and
$G^\prime$ values (Table \ref{modelTest}).

\begin{figure}
\includegraphics[width=\linewidth]{aahxntln}
\caption{Spectra obtained for all-atom (AA) Au / butanethiol / solvent
  systems. When butanethiol is deuterated (lower left), its
  vibrational overlap with hexane decreases significantly.  Since
  aromatic molecules and the butanethiol are vibrationally dissimilar,
  the change is not as dramatic when toluene is the solvent (right).}
\label{aahxntln}
\end{figure}

For the Au / butanethiol / toluene interfaces, having the AA
butanethiol deuterated did not yield a significant change in the
measured conductance. Compared to the C-H vibrational overlap between
hexane and butanethiol, both of which have alkyl chains, the overlap
between toluene and butanethiol is not as significant and thus does
not contribute as much to the heat exchange process. 

Previous observations of Zhang {\it et al.}\cite{hase:2010} indicate
that the {\it intra}molecular heat transport due to alkylthiols is
highly efficient.  Combining our observations with those of Zhang {\it
  et al.}, it appears that butanethiol acts as a channel to expedite
heat flow from the gold surface and into the alkyl chain.  The
vibrational coupling between the metal and the liquid phase can
therefore be enhanced with the presence of suitable capping agents.

Deuterated models in the UA force field did not decouple the thermal
transport as well as in the AA force field.  The UA models, even
though they have eliminated the high frequency C-H vibrational
overlap, still have significant overlap in the lower-frequency
portions of the infrared spectra (Figure \ref{uahxnua}).  Deuterating
the UA models did not decouple the low frequency region enough to
produce an observable difference for the results of $G$ (Table
\ref{modelTest}). 

\begin{figure}
\includegraphics[width=\linewidth]{uahxnua}
\caption{Vibrational power spectra for UA models for the butanethiol
  and hexane solvent (upper panel) show the high degree of overlap
  between these two molecules, particularly at lower frequencies.
  Deuterating a UA model for the solvent (lower panel) does not
  decouple the two spectra to the same degree as in the AA force
  field (see Fig \ref{aahxntln}).}
\label{uahxnua}
\end{figure}

\section{Conclusions}
The NIVS algorithm has been applied to simulations of
butanethiol-capped Au(111) surfaces in the presence of organic
solvents. This algorithm allows the application of unphysical thermal
flux to transfer heat between the metal and the liquid phase. With the
flux applied, we were able to measure the corresponding thermal
gradients and to obtain interfacial thermal conductivities. Under
steady states, 2-3 ns trajectory simulations are sufficient for
computation of this quantity.

Our simulations have seen significant conductance enhancement in the
presence of capping agent, compared with the bare gold / liquid
interfaces. The vibrational coupling between the metal and the liquid
phase is enhanced by a chemically-bonded capping agent. Furthermore,
the coverage percentage of the capping agent plays an important role
in the interfacial thermal transport process. Moderately low coverages
allow higher contact between capping agent and solvent, and thus could
further enhance the heat transfer process, giving a non-monotonic
behavior of conductance with increasing coverage.

Our results, particularly using the UA models, agree well with
available experimental data.  The AA models tend to overestimate the
interfacial thermal conductance in that the classically treated C-H
vibrations become too easily populated. Compared to the AA models, the
UA models have higher computational efficiency with satisfactory
accuracy, and thus are preferable in modeling interfacial thermal
transport.

Of the two definitions for $G$, the discrete form
(Eq. \ref{discreteG}) was easier to use and gives out relatively
consistent results, while the derivative form (Eq. \ref{derivativeG})
is not as versatile. Although $G^\prime$ gives out comparable results
and follows similar trend with $G$ when measuring close to fully
covered or bare surfaces, the spatial resolution of $T$ profile
required for the use of a derivative form is limited by the number of
bins and the sampling required to obtain thermal gradient information.

Vlugt {\it et al.} have investigated the surface thiol structures for
nanocrystalline gold and pointed out that they differ from those of
the Au(111) surface.\cite{landman:1998,vlugt:cpc2007154} This
difference could also cause differences in the interfacial thermal
transport behavior. To investigate this problem, one would need an
effective method for applying thermal gradients in non-planar
(i.e. spherical) geometries. 

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

\end{doublespace}
\end{document}

Revision:	3770
Committed:	Fri Dec 2 20:14:03 2011 UTC (13 years, 8 months ago) by skuang
Content type:	application/x-tex
File size:	52519 byte(s)
Log Message:	start from methodology right onto interfacial paper. add a bib file.