--- interfacial/interfacial.tex	2011/06/29 13:42:01	3728
+++ interfacial/interfacial.tex	2011/09/28 22:26:23	3764
@@ -23,13 +23,13 @@
 \setlength{\belowcaptionskip}{30 pt}
 
 %\renewcommand\citemid{\ } % no comma in optional reference note
-\bibpunct{[}{]}{,}{s}{}{;}
-\bibliographystyle{aip}
+\bibpunct{[}{]}{,}{n}{}{;}
+\bibliographystyle{achemso}
 
 \begin{document}
 
-\title{Simulating interfacial thermal conductance at metal-solvent
-  interfaces: the role of chemical capping agents}
+\title{Simulating Interfacial Thermal Conductance at Metal-Solvent
+  Interfaces: the Role of Chemical Capping Agents}
 
 \author{Shenyu Kuang and J. Daniel
 Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
@@ -44,22 +44,24 @@ Notre Dame, Indiana 46556}
 \begin{doublespace}
 
 \begin{abstract}
-
-We have developed a Non-Isotropic Velocity Scaling algorithm for
-setting up and maintaining stable thermal gradients in non-equilibrium
-molecular dynamics simulations. This approach effectively imposes
-unphysical thermal flux even between particles of different
-identities, conserves linear momentum and kinetic energy, and
-minimally perturbs the velocity profile of a system when compared with
-previous RNEMD methods. We have used this method to simulate thermal
-conductance at metal / organic solvent interfaces both with and
-without the presence of thiol-based capping agents.  We obtained
-values comparable with experimental values, and observed significant
-conductance enhancement with the presence of capping agents. Computed
-power spectra indicate the acoustic impedance mismatch between metal
-and liquid phase is greatly reduced by the capping agents and thus
-leads to higher interfacial thermal transfer efficiency.
-
+  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 the acoustic
+  impedance mismatch between the metal and liquid phases is
+  effectively reduced 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 {\bf vibrational heat diffusion} of solvent
+  molecules that have been in close contact with the capping agent.
 \end{abstract}
 
 \newpage
@@ -71,98 +73,181 @@ leads to higher interfacial thermal transfer efficienc
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Introduction}
-[BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM]
-Interfacial thermal conductance is extensively studied both
-experimentally and computationally, and systems with interfaces
-present are generally heterogeneous. Although interfaces are commonly
-barriers to heat transfer, it has been
-reported\cite{doi:10.1021/la904855s} that under specific circustances,
-e.g. with certain capping agents present on the surface, interfacial
-conductance can be significantly enhanced. However, heat conductance
-of molecular and nano-scale interfaces will be affected by the
-chemical details of the surface and is challenging to
-experimentalist. The lower thermal flux through interfaces is even
-more difficult to measure with EMD and forward NEMD simulation
-methods. Therefore, developing good simulation methods will be
-desirable in order to investigate thermal transport across interfaces.
+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.
 
-Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
+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}
+
+Theoretical and 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. {\bf And
+  Luo {\it et al.} studied the thermal conductance of Au-SAM-Au
+  junctions using the same approach, with comparison to a constant
+  temperature difference method\cite{Luo20101}. While this latter
+  approach establishes more thermal distributions compared to the
+  former RNEMD methods, it does not guarantee momentum or kinetic
+  energy conservations.}
+
+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 more effective thermal exchange
-between particles of different identities, and thus enables extensive
-study of interfacial conductance.
+the two slabs. Furthermore, it allows effective thermal exchange
+between particles of different identities, and thus makes the study of
+interfacial conductance much simpler.
 
+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. {\bf To our knowledge, few previous MD inverstigations
+  have been found to address to these systems yet.} 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}
-\subsection{Algorithm}
-[BACKGROUND FOR MD METHODS]
-There have been many algorithms for computing thermal conductivity
-using molecular dynamics simulations. However, interfacial conductance
-is at least an order of magnitude smaller. This would make the
-calculation even more difficult for those slowly-converging
-equilibrium methods. Imposed-flux non-equilibrium
-methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
-the response of temperature or momentum gradients are easier to
-measure than the flux, if unknown, and thus, is a preferable way to
-the forward NEMD methods. Although the momentum swapping approach for
-flux-imposing 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.
+\subsection{Imposed-Flux Methods in MD Simulations}
+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 in
-non-equilibrium MD simulations is able to impose relatively large
-kinetic energy flux without obvious perturbation to the velocity
-distribution 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 arbitary identity, and
-the flux quantity is not restricted by particle mass difference.
+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 matricies. To
-determine these scaling factors in the matricies, a set of equations
+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 momentum/energy flux satisfaction is
-solved. With the scaling operation applied to the system in a set
-frequency, corresponding momentum/temperature gradients can be built,
-which can be used for computing transportation properties and other
-applications related to momentum/temperature gradients. The NIVS
-algorithm conserves momenta and energy and does not depend on an
-external thermostat. 
+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 interfaces with a relatively low interfacial conductance, the bulk
-regions on either side of an 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:
+\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 -
+  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 imposed non-physical kinetic energy
-transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle
-  T_\mathrm{cold}\rangle}$ are the average observed temperature of the
-two separated phases.
+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} = J_z *(t)*(2 L_x L_y)$.
 
-When the interfacial conductance is {\it not} small, two ways can be
-used to define $G$.
-
-One way is to assume the temperature is discretely different on two
-sides of the interface, $G$ can be calculated with the thermal flux
-applied $J$ and the maximum temperature difference measured along the
-thermal gradient max($\Delta T$), which occurs at the interface, as:
+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}
+  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}
+
+{\bf We attempt another approach by assuming that temperature is
+  continuous and differentiable throughout the space. Given that
+  $\lambda$ is also differentiable, $G$ can be defined as its
+  gradient. This quantity has the same unit as the commonly known $G$,
+  and the maximum of its magnitude denotes where thermal conductivity
+  has the largest change, i.e. the interface. And vector
+  $\nabla\lambda$ is normal to the interface. In a simplified
+  condition here, we have both $\vec{J}$ and the thermal gradient
+  paralell to the $z$ axis and yield the formula used in our
+  computations.}
+(original text)
 The other approach is to assume a continuous temperature profile along
 the thermal gradient axis (e.g. $z$) and define $G$ at the point where
-the magnitude of thermal conductivity $\lambda$ change reach its
+the magnitude of thermal conductivity ($\lambda$) change reaches its
 maximum, given that $\lambda$ is well-defined throughout the space:
 \begin{equation}
 G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
@@ -173,101 +258,108 @@ With the temperature profile obtained from simulations
 \label{derivativeG}
 \end{equation}
 
-With the temperature profile obtained from simulations, one is able to
+With temperature profiles obtained from simulation, one is able to
 approximate the first and second derivatives of $T$ with finite
-difference method and thus calculate $G^\prime$.
+difference methods and calculate $G^\prime$. In what follows, both
+definitions have been used, and are compared in the results.
 
-In what follows, both definitions are used for calculation and comparison.
+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}.
 
-[IMPOSE G DEFINITION INTO OUR SYSTEMS]
-To facilitate the use of the above definitions in calculating $G$ and
-$G^\prime$, we have a metal slab with its (111) surfaces perpendicular
-to the $z$-axis of our simulation cells. With or withour capping
-agents on the surfaces, the metal slab is solvated with organic
-solvents, as illustrated in Figure \ref{demoPic}.
+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]{demoPic}
-\caption{A sample showing how a metal slab has its (111) surface
-  covered by capping agent molecules and solvated by hexane.}
-\label{demoPic}
-\end{figure}
-
-With a simulation cell setup following the above manner, one is able
-to equilibrate the system and impose an unphysical thermal flux
-between the liquid and the metal phase with the NIVS algorithm. Under
-a stablized thermal gradient induced by periodically applying the
-unphysical flux, one is able to obtain a temperature profile and the
-physical thermal flux corresponding to it, which equals to the
-unphysical flux applied by NIVS. These data enables the evaluation of
-the interfacial thermal conductance of a surface. Figure \ref{gradT}
-is an example how those stablized thermal gradient can be used to
-obtain the 1st and 2nd derivatives of the temperature profile.
-
-\begin{figure}
 \includegraphics[width=\linewidth]{gradT}
-\caption{The 1st and 2nd derivatives of temperature profile can be
-  obtained with finite difference approximation.}
+\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{System Geometry}
-In our simulations, Au is used to construct a metal slab with bare
-(111) surface perpendicular to the $z$-axis. Different slab thickness
-(layer numbers of Au) are simulated. This metal slab is first
-equilibrated under normal pressure (1 atm) and a desired
-temperature. After equilibration, butanethiol is used as the capping
-agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
-atoms in the butanethiol molecules would occupy the three-fold sites
-of the surfaces, and the maximal butanethiol capacity on Au surface is
-$1/3$ of the total number of surface Au atoms[CITATION]. A series of
-different coverage surfaces is investigated in order to study the
-relation between coverage and conductance.
+\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.
 
-[COVERAGE DISCRIPTION] However, since the interactions between surface
-Au and butanethiol is non-bonded, the capping agent molecules are
-allowed to migrate to an empty neighbor three-fold site during a
-simulation. Therefore, the initial configuration would not severely
-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. [MAY NEED FIGURES]
+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 are equilibrated
-under canonical ensemble, Packmol\cite{packmol} is used to pack
-organic solvent molecules in the previously vacuum part of the
-simulation cells, which guarantees that short range repulsive
-interactions do not disrupt the simulations. Two solvents are
-investigated, one which has little vibrational overlap with the
-alkanethiol and plane-like shape (toluene), and one which has similar
-vibrational frequencies and chain-like shape ({\it n}-hexane). The
-spacing filled by solvent molecules, i.e. the gap between periodically
-repeated Au-butanethiol surfaces should be carefully chosen so that it
-would not be too short to affect the liquid phase structure, nor too
-long, leading to over cooling (freezing) or heating (boiling) when a
-thermal flux is applied. In our simulations, this spacing is usually
-$35 \sim 60$\AA.
+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 initial configurations generated by Packmol are further
-equilibrated with the $x$ and $y$ dimensions fixed, only allowing
-length scale change in $z$ dimension. This is to ensure that the
-equilibration of liquid phase does not affect the metal crystal
-structure in $x$ and $y$ dimensions. Further equilibration are run
-under NVT and then NVE ensembles.
+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.
 
-After the systems reach equilibrium, NIVS is implemented to impose a
-periodic unphysical thermal flux between the metal and the liquid
-phase. Most of our simulations are under an average temperature of
-$\sim$200K. Therefore, this flux usually comes from the metal to the
+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 excessively low temperature. This induced temperature
-gradient is stablized and the simulation cell is devided evenly into
-N slabs along the $z$-axis and the temperatures of each slab are
-recorded. 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|
+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
@@ -276,243 +368,685 @@ G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}
 \label{derivativeG2}
 \end{equation}
 
+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 various components. Therefore, force field
-parameter descriptions are needed for interactions both between the
-same type of particles and between particles of different species.
+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 \ref{MnM}.}
+\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}. 
+Sutton-Chen potentials.\cite{Chen90}
 
-For both solvent molecules, straight chain {\it n}-hexane and aromatic
-toluene, United-Atom (UA) and All-Atom (AA) models are used
-respectively. The TraPPE-UA
+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, pseudo-atoms are
-located at the carbon centers for alkyl groups. By eliminating
-explicit hydrogen atoms, these models are simple and computationally
-efficient, while maintains good accuracy. [LOW BOILING POINT IS A
-KNOWN PROBLEM FOR TRAPPE-UA ALKANES, NEED MORE DISCUSSION]
-for
-toluene,  force fields are
-used with rigid body constraints applied.[MORE DETAILS NEEDED]
+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.
 
-Besides the TraPPE-UA models, AA models are included in our studies as
-well. For hexane, the OPLS all-atom\cite{OPLSAA} force field is
-used. [MORE DETAILS] 
-For toluene,
+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.
 
-Buatnethiol molecules are used as capping agent for some of our
-simulations. United-Atom\cite{TraPPE-UA.thiols} and All-Atom models
-are respectively used corresponding to the force field type of
-solvent.
+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.
 
-To describe the interactions between metal Au and non-metal capping
-agent and solvent, we refer to Vlugt\cite{vlugt:cpc2007154} and derive
-other interactions which are not parametrized in their work. (can add
-hautman and klein's paper here and more discussion; need to put
-aromatic-metal interaction approximation here)
+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. 
 
-[TABULATED FORCE FIELD PARAMETERS NEEDED]
+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}
 
-\section{Results}
-\subsection{Toluene Solvent}
+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 results (Table \ref{AuThiolToluene}) show a
-significant conductance enhancement compared to the gold/water
-interface without capping agent and agree with available experimental
-data. This indicates that the metal-metal potential, though not
-predicting an accurate bulk metal thermal conductivity, does not
-greatly interfere with the simulation of the thermal conductance
-behavior across a non-metal interface. The solvent model is not
-particularly volatile, so the simulation cell does not expand
-significantly under higher temperature. We did not observe a
-significant conductance decrease when the temperature was increased to
-300K. The results show that the two definitions used for $G$ yield
-comparable values, though $G^\prime$ tends to be smaller.
+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 \ref{MnM} summarizes the ``metal/non-metal'' parameters used in
+our simulations.
 
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
-      \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the Au/butanethiol/toluene interface at
-        different temperatures using a range of energy fluxes.}
-      
-      \begin{tabular}{cccc}
+      \caption{Non-bonded interaction parameters (including cross
+        interactions with Au atoms) for both force fields used in this
+        work.}      
+      \begin{tabular}{lllllll}
         \hline\hline
-        $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
-        (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
+        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
+        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
         \hline
-        200 & 1.86 & 180 & 135 \\
-            & 2.15 & 204 & 113 \\
-            & 3.93 & 175 & 114 \\
-        300 & 1.91 & 143 & 125 \\
-            & 4.19 & 134 & 113 \\
+        United Atom (UA)
+        &CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
+        &CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
+        &CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
+        &CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
+        \hline
+        All Atom (AA) 
+        &CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
+        &CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
+        &CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
+        &HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
+        &CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
+        &HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
+        \hline
+        Both UA and AA 
+        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
         \hline\hline
       \end{tabular}
-      \label{AuThiolToluene}
+      \label{MnM}
     \end{center}
   \end{minipage}
 \end{table*}
 
-\subsection{Hexane Solvent}
 
-Using the united-atom model, different coverages of capping agent,
-temperatures of simulations and numbers of solvent molecules were all
-investigated and Table \ref{AuThiolHexaneUA} shows the results of
-these computations. The number of hexane molecules in our simulations
-does not affect the calculations significantly. However, a very long
-length scale for the thermal gradient axis ($z$) 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, while too
-few solvent molecules would change the normal behavior of the liquid
-phase. Our $N_{hexane}$ values were chosen to ensure that these
-extreme cases did not happen to our simulations.
+\section{Results}
+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.
 
-Table \ref{AuThiolHexaneUA} enables direct comparison between
-different coverages of capping agent, when other system parameters are
-held constant. With high coverage of butanethiol on the gold surface,
-the interfacial thermal conductance is enhanced
-significantly. Interestingly, a slightly lower butanethiol coverage
-leads to a moderately higher conductivity. This is probably due to
-more solvent/capping agent contact when butanethiol molecules are
-not densely packed, which enhances the interactions between the two
-phases and lowers the thermal transfer barrier of this interface.
-% [COMPARE TO AU/WATER IN PAPER]
+\subsection{Effects due to capping agent coverage}
 
-It is also noted that the overall simulation temperature is another
-factor that affects the interfacial thermal conductance. One
-possibility of this effect may be rooted in the decrease in density of
-the liquid phase. We observed that when the average temperature
-increases from 200K to 250K, the bulk hexane density becomes lower
-than experimental value, as the system is equilibrated under NPT
-ensemble. This leads to lower contact between solvent and capping
-agent, and thus lower conductivity.
+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. 
 
-Conductivity values are more difficult to obtain under higher
-temperatures. This is because the Au surface tends to undergo
-reconstructions in relatively high temperatures. Surface Au atoms can
-migrate outward to reach higher Au-S contact; and capping agent
-molecules can be embedded into the surface Au layer due to the same
-driving force. This phenomenon agrees with experimental
-results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. A surface
-fully covered in capping agent is more susceptible to reconstruction,
-possibly because fully coverage prevents other means of capping agent
-relaxation, such as migration to an empty neighbor three-fold site.
+\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}
 
-%MAY ADD MORE DATA TO TABLE
+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.
+
+{\bf MAY NOT NEED $J_z$ IN TABLE}
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
-      \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the Au/butanethiol/hexane interface
-        with united-atom model and different capping agent coverage
-        and solvent molecule numbers at different temperatures using a
-        range of energy fluxes.}
       
-      \begin{tabular}{cccccc}
+      \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; ``Avg.'' denotes results
+        that are averages of simulations under different applied
+        thermal flux $(J_z)$ values. Error estimates are indicated in
+        parentheses.}
+      
+      \begin{tabular}{llccc}
         \hline\hline
-        Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\
-        coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) &
+        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
+        (or bare surface) & model & (GW/m$^2$) &
         \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        0.0   & 200 & 200 & 0.96 & 43.3 & 42.7 \\
-              &     &     & 1.91 & 45.7 & 42.9 \\
-              &     & 166 & 0.96 & 43.1 & 53.4 \\
-        88.9  & 200 & 166 & 1.94 & 172  & 108  \\
-        100.0 & 250 & 200 & 0.96 & 81.8 & 67.0 \\
-              &     & 166 & 0.98 & 79.0 & 62.9 \\
-              &     &     & 1.44 & 76.2 & 64.8 \\
-              & 200 & 200 & 1.92 & 129  & 87.3 \\
-              &     &     & 1.93 & 131  & 77.5 \\
-              &     & 166 & 0.97 & 115  & 69.3 \\
-              &     &     & 1.94 & 125  & 87.1 \\
+        UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
+              & UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
+              & AA hexane    & Avg. & 131(6)    & 122(10)   \\
+              & UA toluene   & 1.96 & 187(16)   & 151(11)   \\
+              & AA toluene   & 1.89 & 200(36)   & 149(53)   \\
+        \hline
+        AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
+              & AA hexane    & Avg. & 442(14)   & 356(31)   \\
+              & AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
+              & UA toluene   & 1.98 & 125(25)   & 97(60)    \\
+              & AA toluene   & 3.79 & 487(56)   & 290(42)   \\
+        \hline
+        AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
+              & AA hexane    & 1.92 & 243(29)   & 191(11)   \\
+              & AA toluene   & 1.93 & 364(36)   & 322(67)   \\
+        \hline
+        bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
+              & UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
+              & AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
+              & UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
         \hline\hline
       \end{tabular}
-      \label{AuThiolHexaneUA}
+      \label{modelTest}
     \end{center}
   \end{minipage}
 \end{table*}
 
-For the all-atom model, the liquid hexane phase was not stable under NPT
-conditions. Therefore, the simulation length scale parameters are
-adopted from previous equilibration results of the united-atom model
-at 200K. Table \ref{AuThiolHexaneAA} shows the results of these
-simulations. The conductivity values calculated with full capping
-agent coverage are substantially larger than observed in the
-united-atom model, and is even higher than predicted by
-experiments. It is possible that our parameters for metal-non-metal
-particle interactions lead to an overestimate of the interfacial
-thermal conductivity, although the active C-H vibrations in the
-all-atom model (which should not be appreciably populated at normal
-temperatures) could also account for this high conductivity. The major
-thermal transfer barrier of Au/butanethiol/hexane interface is between
-the liquid phase and the capping agent, so extra degrees of freedom
-such as the C-H vibrations could enhance heat exchange between these
-two phases and result in a much higher conductivity.
+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, and
+is a typical case for acoustic impedance mismatch between these two
+phases.  
+
+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 (Table \ref{AuThiolHexaneUA}).
+
+\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
+\ref{AuThiolHexaneUA} and \ref{AuThiolToluene}. 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.  
+
+{\bf MAY MOVE TO SUPPORT INFO}
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
+      \caption{In the hexane-solvated interfaces, the system size has
+        little effect on the calculated values for interfacial
+        conductance ($G$ and $G^\prime$), but the direction of heat
+        flow (i.e. the sign of $J_z$) can alter the average
+        temperature of the liquid phase and this can alter the
+        computed conductivity.}
       
-      \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the Au/butanethiol/hexane interface
-        with all-atom model and different capping agent coverage at
-        200K using a range of energy fluxes.}
-      
-      \begin{tabular}{cccc}
+      \begin{tabular}{ccccccc}
         \hline\hline
-        Thiol & $J_z$ & $G$ & $G^\prime$ \\
-        coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        $\langle T\rangle$ & $N_{hexane}$  & $\rho_{hexane}$ &
+        $J_z$ & $G$ & $G^\prime$ \\
+        (K) &  & (g/cm$^3$) & (GW/m$^2$) &
+        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        0.0   & 0.95 & 28.5 & 27.2 \\
-              & 1.88 & 30.3 & 28.9 \\
-        100.0 & 2.87 & 551  & 294  \\
-              & 3.81 & 494  & 193  \\
+        200 & 266 &  0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
+            & 200 &  0.688 &  0.96 & 125(16)   & 90.2(15)  \\
+            &     &        &  1.91 & 139(10)   & 101(10)   \\
+            &     &        &  2.83 & 141(6)    & 89.9(9.8) \\
+            & 166 &  0.681 &  0.97 & 141(30)   & 78(22)    \\
+            &     &        &  1.92 & 138(4)    & 98.9(9.5) \\
+        \hline
+        250 & 200 &  0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
+            &     &        & -0.95 & 49.4(0.3) & 45.7(2.1) \\
+            & 166 &  0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
+            &     &        &  1.44 & 76.2(5.0) & 64.8(3.8) \\
+            &     &        & -0.95 & 56.4(2.5) & 54.4(1.1) \\
+            &     &        & -1.85 & 47.8(1.1) & 53.5(1.5) \\
         \hline\hline
       \end{tabular}
-      \label{AuThiolHexaneAA}
+      \label{AuThiolHexaneUA}
     \end{center}
   \end{minipage}
 \end{table*}
 
-%subsubsection{Vibrational spectrum study on conductance mechanism}
-To investigate the mechanism of this interfacial thermal conductance,
-the vibrational spectra of various gold systems were obtained and are
-shown as in the upper panel of Fig. \ref{vibration}. To obtain these
-spectra, one first runs a simulation in the NVE ensemble and collects
-snapshots of configurations; these configurations are used to compute
-the velocity auto-correlation functions, which is used to construct a
-power spectrum via a Fourier transform. The gold surfaces covered by
-butanethiol molecules exhibit an additional peak observed at a
-frequency of $\sim$170cm$^{-1}$, which is attributed to the vibration
-of the S-Au bond. This vibration enables efficient thermal transport
-from surface Au atoms to the capping agents. Simultaneously, as shown
-in the lower panel of Fig. \ref{vibration}, the large overlap of the
-vibration spectra of butanethiol and hexane in the all-atom model,
-including the C-H vibration, also suggests high thermal exchange
-efficiency. The combination of these two effects produces the drastic
-interfacial thermal conductance enhancement in the all-atom model.
+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 as shown in
+Table \ref{AuThiolHexaneUA}.  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  (Table
+\ref{AuThiolToluene}).
+
+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.
+
+{\bf MAY MOVE TO SUPPORT INFO}
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      \caption{When toluene is the solvent, the interfacial thermal
+        conductivity is less sensitive to temperature, but again, the
+        direction of the heat flow can alter the solvent temperature
+        and can change the computed conductance values.}
+      
+      \begin{tabular}{ccccc}
+        \hline\hline
+        $\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
+        (K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        \hline
+        200 & 0.933 &  2.15 & 204(12) & 113(12) \\
+            &       & -1.86 & 180(3)  & 135(21) \\
+            &       & -3.93 & 176(5)  & 113(12) \\
+        \hline
+        300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
+            &       & -4.19 & 135(9)  & 113(12) \\
+        \hline\hline
+      \end{tabular}
+      \label{AuThiolToluene}
+    \end{center}
+  \end{minipage}
+\end{table*}
+
+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}
+
+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. {\bf This confirms the results from Luo {\it et
+    al.}\cite{Luo20101}, which reported $G$ for Au-SAM junctions
+  generally twice larger than what we have computed for the
+  thiol-liquid interfaces.}
+
 \begin{figure}
 \includegraphics[width=\linewidth]{vibration}
-\caption{Vibrational spectra obtained for gold in different
-  environments (upper panel) and for Au/thiol/hexane simulation in
-  all-atom model (lower panel).}
-\label{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}
-% 600dpi, letter size. too large?
 
+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. {\bf Due to the classical models used, this
+  even includes those high frequency modes which should be unpopulated
+  at our relatively low temperatures. This artifact causes high
+  frequency vibrations accountable for thermal transport in classical
+  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
+acoustic impedance mismatch between the metal and the liquid phase can
+therefore be effectively reduced 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 acoustic impedance mismatch between the metal and the
+liquid phase is effectively eliminated 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
+Dame. 
 
+\section{Supporting Information}
+This information is available free of charge via the Internet at
+http://pubs.acs.org.
+
+\newpage
+
 \bibliography{interfacial}
 
 \end{doublespace}
