--- interfacial/interfacial.tex	2011/07/14 21:55:24	3741
+++ interfacial/interfacial.tex	2011/09/30 19:37:13	3767
@@ -28,8 +28,8 @@
 
 \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,24 +44,27 @@ Notre Dame, Indiana 46556}
 \begin{doublespace}
 
 \begin{abstract}
+  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 vibrational
+  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.
 
-With the Non-Isotropic Velocity Scaling algorithm (NIVS) we have
-developed, an unphysical thermal flux can be effectively set up even
-for non-homogeneous systems like interfaces in non-equilibrium
-molecular dynamics simulations. In this work, this algorithm is
-applied for simulating thermal conductance at metal / organic solvent
-interfaces with various coverages of butanethiol capping
-agents. Different solvents and force field models were tested. Our
-results suggest that the United-Atom models are able to provide an
-estimate of the interfacial thermal conductivity comparable to
-experiments in our simulations with satisfactory computational
-efficiency. From our results, the acoustic impedance mismatch between
-metal and liquid phase is effectively reduced by the capping
-agents, and thus leads to interfacial thermal conductance
-enhancement. Furthermore, this effect is closely related to the
-capping agent coverage on the metal surfaces and the type of solvent
-molecules, and is affected by the models used in the simulations.
-
+Keywords: non-equilibrium, molecular dynamics, vibrational overlap,
+coverage dependent.
 \end{abstract}
 
 \newpage
@@ -73,50 +76,87 @@ Interfacial thermal conductance is extensively studied
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Introduction}
-Interfacial thermal conductance is extensively studied both
-experimentally and computationally\cite{cahill:793}, due to its
-importance in nanoscale science and technology. Reliability of
-nanoscale devices depends on their thermal transport
-properties. Unlike bulk homogeneous materials, nanoscale materials
-features significant presence of interfaces, and these interfaces
-could dominate the heat transfer behavior of these
-materials. Furthermore, these materials are generally heterogeneous,
-which challenges traditional research methods for homogeneous
-systems.
+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.
 
-Heat conductance of molecular and nano-scale interfaces will be
-affected by the chemical details of the surface. Experimentally,
-various interfaces have been investigated for their thermal
-conductance properties. 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 CTAB on thermal transport between gold nanorods
-and solvent\cite{doi:10.1021/jp8051888}; Juv\'e {\it et al.} studied
-the cooling dynamics, which is controlled by thermal interface
-resistence of glass-embedded metal
-nanoparticles\cite{PhysRevB.80.195406}. Although interfaces are
-commonly 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}.
+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
+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{schwartz} 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 one might apply the acoustic mismatch model to a
+chemically-modified surface, particularly when the acoustic properties
+of a monolayer film may not be well characterized.
+
+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
+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 has yet to be studied. 
-The comparatively low thermal flux through interfaces is
-difficult to measure with Equilibrium MD or forward NEMD simulation
-methods. Therefore, the Reverse NEMD (RNEMD) methods would have the
-advantage of having this difficult to measure flux known when studying
-the thermal transport across interfaces, given that the simulation
-methods being able to effectively apply an unphysical flux in
-non-homogeneous systems.
+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 the Non-Isotropic Velocity Scaling (NIVS)
+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
@@ -131,39 +171,39 @@ underlying mechanism for this phenomena was investigat
 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 this phenomena was investigated.
+underlying mechanism for the phenomena was investigated.
 
-[MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]
-
 \section{Methodology}
-\subsection{Imposd-Flux Methods in MD Simulations}
-For systems with low interfacial conductivity one must have a method
-capable of generating relatively small fluxes, compared to those
-required for bulk conductivity. This requirement makes the calculation
-even more difficult for those slowly-converging equilibrium
-methods\cite{Viscardy:2007lq}.
-Forward methods impose gradient, but in interfacail conditions it is
-not clear what behavior to impose at the boundary...
- Imposed-flux reverse non-equilibrium
-methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
-the thermal response becomes easier to
-measure than the flux. 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.
+\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 to
-non-equilibrium MD simulations is able to impose a wide range of
+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 arbitary identity, and
+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 energy flux satisfaction is solved. With the
 scaling operation applied to the system in a set frequency, bulk
@@ -171,151 +211,173 @@ momenta and energy and does not depend on an external 
 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, there are two
-ways to define $G$.
-
-One 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 deviding surface,
-as:
+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}
 
-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
-maximum, given that $\lambda$ is well-defined throughout the space:
-\begin{equation}
-G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
-         = \Big|\frac{\partial}{\partial z}\left(-J_z\Big/
-           \left(\frac{\partial T}{\partial z}\right)\right)\Big|
-         = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
-         \Big/\left(\frac{\partial T}{\partial z}\right)^2
-\label{derivativeG}
-\end{equation}
-
-With the temperature profile obtained from simulations, one is able to
-approximate the first and second derivatives of $T$ with finite
-difference methods and thus calculate $G^\prime$.
-
-In what follows, both definitions have been used for calculation and
-are compared in the results.
-
-To compare the above definitions ($G$ and $G^\prime$), we have modeled
-a metal slab with its (111) surfaces perpendicular to the $z$-axis of
-our simulation cells. Both with and withour capping agents on the
-surfaces, the metal slab is solvated with simple organic solvents, as
-illustrated in Figure \ref{demoPic}.
-
 \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 response or a 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.}
+  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 are able to obtain a temperature profile and
-its spatial derivatives. These quantities enable the evaluation of the
-interfacial thermal conductance of a surface. Figure \ref{gradT} is an
-example how those applied thermal fluxes can be used to obtain the 1st
-and 2nd derivatives of the temperature profile.
+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{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{Simulation Protocol}
 The NIVS algorithm has been implemented in our MD simulation code,
-OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
-simulations. Different slab thickness (layer numbers of Au) were
-simulated. Metal slabs were first equilibrated under atmospheric
-pressure (1 atm) and a desired temperature (e.g. 200K). After
-equilibration, butanethiol capping agents were placed at three-fold
-sites on the Au(111) surfaces. The maximum butanethiol capacity on Au
-surface is $1/3$ of the total number of surface Au
-atoms\cite{vlugt:cpc2007154}. A series of different coverages was
-investigated in order to study the relation between coverage and
-interfacial conductance.
+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 would not noticeably affect the sampling of a
+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. [MAY NEED FIGURES]
+configurations explored in the simulations.
 
-After the modified Au-butanethiol surface systems were equilibrated
-under canonical ensemble, organic solvent molecules were packed in the
-previously empty part of the simulation cells\cite{packmol}. Two
+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 a planar shape (toluene), and one which has
-similar vibrational frequencies and chain-like shape ({\it n}-hexane).
+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 space filled by solvent molecules, i.e. the gap between
-periodically repeated Au-butanethiol surfaces should be carefully
-chosen. A very long length scale for the thermal gradient axis ($z$)
-may cause excessively hot or cold temperatures in the middle of the
+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. And the
-corresponding spacing is usually $35 \sim 60$\AA.
+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 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 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 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
+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
@@ -323,100 +385,115 @@ G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}
          \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 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.
 
-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}. 
-
-Figure \ref{demoMol} demonstrates how we name our pseudo-atoms of the
-organic solvent molecules in our simulations.
-
 \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
-  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} (UA) and
-  \protect\cite{OPLSAA} (AA).  Cross-interactions with the Au atoms are given
-  in Table \ref{MnM}.}
+  (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}
 
-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
+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, 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. However, the TraPPE-UA for
-alkanes is known to predict a lower boiling point than experimental
-values. Considering that after an unphysical thermal flux is applied
-to a system, the temperature of ``hot'' area in the liquid phase would be
-significantly higher than the average, to prevent over heating and
-boiling of the liquid phase, the average temperature in our
-simulations should be much lower than the liquid boiling point. [MORE DISCUSSION]
-For UA-toluene model, rigid body constraints are applied, so that the
-benzene ring and the methyl-CRar bond are kept rigid. This would save
-computational time.[MORE DETAILS]
+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. For hexane, the OPLS-AA\cite{OPLSAA}
-force field is used. [MORE DETAILS]
-For toluene, the United Force Field developed by Rapp\'{e} {\it et
-  al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]
+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 capping agent in our simulations, the butanethiol molecules can
-either use UA or AA model. The TraPPE-UA force fields includes
+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 adapt this
-change and derive suitable parameters for butanethiol adsorbed on
-Au(111) surfaces, we adopt the S parameters from Luedtke and
-Landman\cite{landman:1998} and modify parameters for its neighbor C
-atom for charge balance in the molecule. Note that the model choice
-(UA or AA) of 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:
+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}}
+  \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 capping
-agent and solvent particles, 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\cite{hautman:4994} for the
-Au(111) surface. As our simulations require the gold lattice slab to
-be non-rigid so that it could accommodate kinetic energy for thermal
-transport study purpose, the pair-wise form of potentials is
-preferred.
+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.
 
-Besides, 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.[MORE DETAILS]
+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.
 
-However, the Lennard-Jones parameters between Au and other types of
-particles in our simulations are not yet well-established. For these
-interactions, we attempt to derive their parameters using the Mixing
-Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters
-for Au is first extracted from the Au-CH$_x$ parameters by applying
-the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
-parameters in our simulations.
-
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
@@ -443,7 +520,8 @@ parameters in our simulations.
         &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   \\
+        Both UA and AA 
+        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
         \hline\hline
       \end{tabular}
       \label{MnM}
@@ -452,317 +530,131 @@ parameters in our simulations.
 \end{table*}
 
 
-\section{Results and Discussions}
-[MAY HAVE A BRIEF SUMMARY]
-\subsection{How Simulation Parameters Affects $G$}
-[MAY NOT PUT AT FIRST]
-We have varied our protocol or other parameters of the simulations in
-order to investigate how these factors would affect the measurement of
-$G$'s. It turned out that while some of these parameters would not
-affect the results substantially, some other changes to the
-simulations would have a significant impact on the measurement
-results.
+\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.
 
-In some of our simulations, we allowed $L_x$ and $L_y$ to change
-during equilibrating the liquid phase. Due to the stiffness of the Au
-slab, $L_x$ and $L_y$ would not change noticeably after
-equilibration. Although $L_z$ could fluctuate $\sim$1\% after a system
-is fully equilibrated in the NPT ensemble, this fluctuation, as well
-as those comparably smaller to $L_x$ and $L_y$, would not be magnified
-on the calculated $G$'s, as shown in Table \ref{AuThiolHexaneUA}. This
-insensivity to $L_i$ fluctuations allows reliable measurement of $G$'s
-without the necessity of extremely cautious equilibration process. 
+\subsection{Effects due to capping agent coverage}
 
-As stated in our computational details, the spacing filled with
-solvent molecules can be chosen within a range. This allows some
-change of solvent molecule numbers for the same Au-butanethiol
-surfaces. We did this study on our Au-butanethiol/hexane
-simulations. Nevertheless, the results obtained from systems of
-different $N_{hexane}$ did not indicate that the measurement of $G$ is
-susceptible to this parameter. For computational efficiency concern,
-smaller system size would be preferable, given that the liquid phase
-structure is not affected.
-
-Our NIVS algorithm allows change of unphysical thermal flux both in
-direction and in quantity. This feature extends our investigation of
-interfacial thermal conductance. However, the magnitude of this
-thermal flux is not arbitary if one aims to obtain a stable and
-reliable thermal gradient. A temperature profile would be
-substantially affected by noise when $|J_z|$ has a much too low
-magnitude; while an excessively large $|J_z|$ that overwhelms the
-conductance capacity of the interface would prevent a thermal gradient
-to reach a stablized steady state. NIVS has the advantage of allowing
-$J$ to vary in a wide range such that the optimal flux range for $G$
-measurement can generally be simulated by the algorithm. Within the
-optimal range, we were able to study how $G$ would change according to
-the thermal flux across the interface. For our simulations, we denote
-$J_z$ to be positive when the physical thermal flux is from the liquid
-to metal, and negative vice versa. The $G$'s measured under different
-$J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These
-results do not suggest that $G$ is dependent 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
-couple $J_z$'s and do not need to test a large series of fluxes.
-
-%ADD MORE TO TABLE
-\begin{table*}
-  \begin{minipage}{\linewidth}
-    \begin{center}
-      \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the 100\% covered Au-butanethiol/hexane
-        interfaces with UA model and different hexane molecule numbers
-        at different temperatures using a range of energy fluxes.}
-      
-      \begin{tabular}{ccccccc}
-        \hline\hline
-        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
-        $J_z$ & $G$ & $G^\prime$ \\
-        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
-        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
-        \hline
-        200 & 266 & No  & 0.672 & -0.96 & 102()  & 80.0() \\
-            & 200 & Yes & 0.694 &  1.92 & 129()  & 87.3() \\
-            &     & Yes & 0.672 &  1.93 & 131()  & 77.5() \\
-            &     & No  & 0.688 &  0.96 & 125()  & 90.2() \\
-            &     &     &       &  1.91 & 139()  & 101()  \\
-            &     &     &       &  2.83 & 141()  & 89.9() \\
-            & 166 & Yes & 0.679 &  0.97 & 115()  & 69.3() \\
-            &     &     &       &  1.94 & 125()  & 87.1() \\
-            &     & No  & 0.681 &  0.97 & 141()  & 77.7() \\
-            &     &     &       &  1.92 & 138()  & 98.9() \\
-        \hline
-        250 & 200 & No  & 0.560 &  0.96 & 74.8() & 61.8() \\
-            &     &     &       & -0.95 & 49.4() & 45.7() \\
-            & 166 & Yes & 0.570 &  0.98 & 79.0() & 62.9() \\
-            &     & No  & 0.569 &  0.97 & 80.3() & 67.1() \\
-            &     &     &       &  1.44 & 76.2() & 64.8() \\
-            &     &     &       & -0.95 & 56.4() & 54.4() \\
-            &     &     &       & -1.85 & 47.8() & 53.5() \\
-        \hline\hline
-      \end{tabular}
-      \label{AuThiolHexaneUA}
-    \end{center}
-  \end{minipage}
-\end{table*}
-
-Furthermore, we also attempted to increase system average temperatures
-to above 200K. These simulations are first equilibrated in the NPT
-ensemble under normal pressure. As stated above, the TraPPE-UA model
-for hexane tends to predict a lower boiling point. In our simulations,
-hexane had diffculty to remain in liquid phase when NPT equilibration
-temperature is higher than 250K. Additionally, the equilibrated liquid
-hexane density under 250K becomes lower than experimental value. This
-expanded liquid phase leads to lower contact between hexane and
-butanethiol as well.[MAY NEED FIGURE] And this reduced contact would
-probably be accountable for a lower interfacial thermal conductance,
-as shown in Table \ref{AuThiolHexaneUA}. 
-
-A similar study for TraPPE-UA toluene agrees with the above result as
-well. Having a higher boiling point, toluene tends to remain liquid in
-our simulations even equilibrated under 300K in NPT
-ensembles. Furthermore, the expansion of the toluene liquid phase is
-not as significant as that of the hexane. This prevents severe
-decrease of liquid-capping agent contact and the results (Table
-\ref{AuThiolToluene}) show only a slightly decreased interface
-conductance. Therefore, solvent-capping agent contact should play an
-important role in the thermal transport process across the interface
-in that higher degree of contact could yield increased conductance.
-
-[ADD ERROR ESTIMATE TO TABLE]
-\begin{table*}
-  \begin{minipage}{\linewidth}
-    \begin{center}
-      \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
-        interface at different temperatures using a range of energy
-        fluxes.}
-      
-      \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 & -1.86 & 180() & 135() \\
-            &       &  2.15 & 204() & 113() \\
-            &       & -3.93 & 175() & 114() \\
-        \hline
-        300 & 0.855 & -1.91 & 143() & 125() \\
-            &       & -4.19 & 134() & 113() \\
-        \hline\hline
-      \end{tabular}
-      \label{AuThiolToluene}
-    \end{center}
-  \end{minipage}
-\end{table*}
+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. 
 
-Besides lower interfacial thermal conductance, surfaces in relatively
-high temperatures are susceptible to reconstructions, when
-butanethiols have a full coverage on the Au(111) surface. These
-reconstructions include surface Au atoms migrated outward to the S
-atom layer, and butanethiol molecules embedded into the original
-surface Au layer. The driving force for this behavior is the strong
-Au-S interactions in our simulations. And these reconstructions lead
-to higher ratio of Au-S attraction and thus is energetically
-favorable. Furthermore, this phenomenon agrees with experimental
-results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. Vlugt
-{\it et al.} had kept their Au(111) slab rigid so that their
-simulations can reach 300K without surface reconstructions. Without
-this practice, simulating 100\% thiol covered interfaces under higher
-temperatures could hardly avoid surface reconstructions. However, our
-measurement is based on assuming homogeneity on $x$ and $y$ dimensions
-so that measurement of $T$ at particular $z$ would be an effective
-average of the particles of the same type. Since surface
-reconstructions could eliminate the original $x$ and $y$ dimensional
-homogeneity, measurement of $G$ is more difficult to conduct under
-higher temperatures. Therefore, most of our measurements are
-undertaken at $\langle T\rangle\sim$200K.
+\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}
 
-However, when the surface is not completely covered by butanethiols,
-the simulated system is more resistent to the reconstruction
-above. Our Au-butanethiol/toluene system did not see this phenomena
-even at $\langle T\rangle\sim$300K. The Au(111) surfaces have a 90\%
-coverage of butanethiols and have empty three-fold sites. These empty
-sites could help prevent surface reconstruction in that they provide
-other means of capping agent relaxation. It is observed that
-butanethiols can migrate to their neighbor empty sites during a
-simulation. Therefore, we were able to obtain $G$'s for these
-interfaces even at a relatively high temperature without being
-affected by surface reconstructions.
+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.
 
-\subsection{Influence of Capping Agent Coverage on $G$}
-To investigate the influence of butanethiol coverage on interfacial
-thermal conductance, a series of different coverage Au-butanethiol
-surfaces is prepared and solvated with various organic
-molecules. These systems are then equilibrated and their interfacial
-thermal conductivity are measured with our NIVS algorithm. Table
-\ref{tlnUhxnUhxnD} lists these results for direct comparison between
-different coverages of butanethiol. To study the isotope effect in
-interfacial thermal conductance, deuterated UA-hexane is included as
-well.
+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.
 
-It turned out that with partial covered butanethiol on the Au(111)
-surface, the derivative definition for $G$ (Eq. \ref{derivativeG}) has
-difficulty to apply, due to the difficulty in locating the maximum of
-change of $\lambda$. Instead, the discrete definition
-(Eq. \ref{discreteG}) is easier to apply, as max($\Delta T$) can still
-be well-defined. Therefore, $G$'s (not $G^\prime$) are used for this
-section. 
+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.
 
-From Table \ref{tlnUhxnUhxnD}, one can see the significance of the
-presence of capping agents. Even when a fraction of the Au(111)
-surface sites are covered with butanethiols, the conductivity would
-see an enhancement by at least a factor of 3. This indicates the
-important role cappping agent is playing for thermal transport
-phenomena on metal/organic solvent surfaces.
+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.
 
-Interestingly, as one could observe from our results, the maximum
-conductance enhancement (largest $G$) happens while the surfaces are
-about 75\% covered with butanethiols. This again indicates that
-solvent-capping agent contact has an important role of the thermal
-transport process. Slightly lower butanethiol coverage allows small
-gaps between butanethiols to form. And these gaps could be filled with
-solvent molecules, which acts like ``heat conductors'' on the
-surface. The higher degree of interaction between these solvent
-molecules and capping agents increases the enhancement effect and thus
-produces a higher $G$ than densely packed butanethiol arrays. However,
-once this maximum conductance enhancement is reached, $G$ decreases
-when butanethiol coverage continues to decrease. Each capping agent
-molecule reaches its maximum capacity for thermal
-conductance. Therefore, even higher solvent-capping agent contact
-would not offset this effect. Eventually, when butanethiol coverage
-continues to decrease, solvent-capping agent contact actually
-decreases with the disappearing of butanethiol molecules. In this
-case, $G$ decrease could not be offset but instead accelerated.
+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.
 
-A comparison of the results obtained from differenet organic solvents
-can also provide useful information of the interfacial thermal
-transport process. The deuterated hexane (UA) results do not appear to
-be much different from those of normal hexane (UA), given that
-butanethiol (UA) is non-deuterated for both solvents. These UA model
-studies, even though eliminating C-H vibration samplings, still have
-C-C vibrational frequencies different from each other. However, these
-differences in the infrared range do not seem to produce an observable
-difference for the results of $G$. [MAY NEED FIGURE]
+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.
 
-Furthermore, results for rigid body toluene solvent, as well as other
-UA-hexane solvents, are reasonable within the general experimental
-ranges[CITATIONS]. This suggests that explicit hydrogen might not be a
-required factor for modeling thermal transport phenomena of systems
-such as Au-thiol/organic solvent.
+\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.
 
-However, results for Au-butanethiol/toluene do not show an identical
-trend with those for Au-butanethiol/hexane in that $G$'s remain at
-approximately the same magnitue when butanethiol coverage differs from
-25\% to 75\%. This might be rooted in the molecule shape difference
-for plane-like toluene and chain-like {\it n}-hexane. Due to this
-difference, toluene molecules have more difficulty in occupying
-relatively small gaps among capping agents when their coverage is not
-too low. Therefore, the solvent-capping agent contact may keep
-increasing until the capping agent coverage reaches a relatively low
-level. This becomes an offset for decreasing butanethiol molecules on
-its effect to the process of interfacial thermal transport. Thus, one
-can see a plateau of $G$ vs. butanethiol coverage in our results.
-
-[NEED ERROR ESTIMATE]
-\begin{figure}
-\includegraphics[width=\linewidth]{coverage}
-\caption{Comparison of interfacial thermal conductivity ($G$) values
-  for the Au-butanethiol/solvent interface with various UA models and
-  different capping agent coverages at $\langle T\rangle\sim$200K
-  using certain energy flux respectively.}
-\label{coverage}
-\end{figure}
-
-\subsection{Influence of Chosen Molecule Model on $G$}
-[MAY COMBINE W MECHANISM STUDY]
-
-In addition to UA solvent/capping agent models, AA models are included
-in our simulations as well. Besides simulations of the same (UA or AA)
-model for solvent and capping agent, different models can be applied
-to different components. Furthermore, regardless of models chosen,
-either the solvent or the capping agent can be deuterated, similar to
-the previous section. Table \ref{modelTest} summarizes the results of
-these studies.
-
-[MORE DATA; ERROR ESTIMATE]
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
       
-      \caption{Computed interfacial thermal conductivity ($G$ and
+      \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. (D stands for deuterated solvent
-        or capping agent molecules; ``Avg.'' denotes results that are
-        averages of several simulations.)}
+        $\langle T\rangle\sim$200K.  Here ``D'' stands for deuterated
+        solvent or capping agent molecules. Error estimates are
+        indicated in parentheses.}
       
-      \begin{tabular}{ccccc}
+      \begin{tabular}{llccc}
         \hline\hline
-        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
-        (or bare surface) & model & (GW/m$^2$) &
+        Butanethiol model & Solvent & $G$ & $G^\prime$ \\
+        (or bare surface) & model & 
         \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        UA    & UA hexane    & Avg. & 131()  & 86.5() \\
-              & UA hexane(D) & 1.95 & 153()  & 136()  \\
-              & AA hexane    & 1.94 & 135()  & 129()  \\
-              &              & 2.86 & 126()  & 115()  \\
-              & UA toluene   & 1.96 & 187()  & 151()  \\
-              & AA toluene   & 1.89 & 200()  & 149()  \\
+        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    & 1.94 & 116()  & 129()  \\
-              & AA hexane    & Avg. & 442()  & 356()  \\
-              & AA hexane(D) & 1.93 & 222()  & 234()  \\
-              & UA toluene   & 1.98 & 125()  & 96.5() \\
-              & AA toluene   & 3.79 & 487()  & 290()  \\
+        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    & 1.94 & 158()  & 172()  \\
-              & AA hexane    & 1.92 & 243()  & 191()  \\
-              & AA toluene   & 1.93 & 364()  & 322()  \\
+        AA(D) & UA hexane    & 158(25)   & 172(4)    \\
+              & AA hexane    & 243(29)   & 191(11)   \\
+              & AA toluene   & 364(36)   & 322(67)   \\
         \hline
-        bare  & UA hexane    & Avg. & 46.5() & 49.4() \\
-              & UA hexane(D) & 0.98 & 43.9() & 43.0() \\
-              & AA hexane    & 0.96 & 31.0() & 29.4() \\
-              & UA toluene   & 1.99 & 70.1() & 65.8() \\
+        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}
@@ -770,129 +662,337 @@ To facilitate direct comparison, the same system with 
   \end{minipage}
 \end{table*}
 
-To facilitate direct comparison, the same system with differnt models
-for different components uses the same length scale for their
-simulation cells. Without the presence of capping agent, using
-different models for hexane yields similar results for both $G$ and
-$G^\prime$, and these two definitions agree with eath other very
-well. This indicates very weak interaction between the metal and the
-solvent, and is a typical case for acoustic impedance mismatch between
-these two phases.
+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.
 
-As for Au(111) surfaces completely covered by butanethiols, the choice
-of models for capping agent and solvent could impact the measurement
-of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
-interfaces, using AA model for both butanethiol and hexane yields
-substantially higher conductivity values than using UA model for at
-least one component of the solvent and capping agent, which exceeds
-the upper bond of experimental value range. This is probably due to
-the classically treated C-H vibrations in the AA model, which should
-not be appreciably populated at normal temperatures. In comparison,
-once either the hexanes or the butanethiols are deuterated, one can
-see a significantly lower $G$ and $G^\prime$. In either of these
-cases, the C-H(D) vibrational overlap between the solvent and the
-capping agent is removed. [MAY NEED FIGURE] Conclusively, the
-improperly treated C-H vibration in the AA model produced
-over-predicted results accordingly. Compared to the AA model, the UA
-model yields more reasonable results with higher computational
-efficiency.
+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.
 
-However, for Au-butanethiol/toluene interfaces, having the AA
-butanethiol deuterated did not yield a significant change in the
-measurement results. Compared to the C-H vibrational overlap between
-hexane and butanethiol, both of which have alkyl chains, that overlap
-between toluene and butanethiol is not so significant and thus does
-not have as much contribution to the ``Intramolecular Vibration
-Redistribution''[CITE HASE]. Conversely, extra degrees of freedom such
-as the C-H vibrations could yield higher heat exchange rate between
-these two phases and result in a much higher conductivity.
+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.
 
-Although the QSC model for Au is known to predict an overly low value
-for bulk metal gold conductivity\cite{kuang:164101}, our computational
-results for $G$ and $G^\prime$ do not seem to be affected by this
-drawback of the model for metal. Instead, our results suggest that the
-modeling of interfacial thermal transport behavior relies mainly on
-the accuracy of the interaction descriptions between components
-occupying the interfaces.
-
-\subsection{Mechanism of Interfacial Thermal Conductance Enhancement
-  by Capping Agent}
-%OR\subsection{Vibrational spectrum study on conductance mechanism}
-
-[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL, EQN'S]
-
-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.
-
-[MAY RELATE TO HASE'S]
- The gold surfaces covered by
-butanethiol molecules, compared to bare gold surfaces, 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.
-
-[REDO. MAY NEED TO CONVERT TO JPEG]
+\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 1 and 2 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}
+
+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{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}
 
-[COMPARISON OF TWO G'S; AU SLAB WIDTHS; ETC]
-% The results show that the two definitions used for $G$ yield
-% comparable values, though $G^\prime$ tends to be smaller.
+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 we developed has been applied to simulations of
-Au-butanethiol surfaces with organic solvents. This algorithm allows
-effective unphysical thermal flux transferred between the metal and
-the liquid phase. With the flux applied, we were able to measure the
-corresponding thermal gradient and to obtain interfacial thermal
-conductivities. Our simulations have seen significant conductance
-enhancement with the presence of capping agent, compared to the bare
-gold/liquid interfaces. The acoustic impedance mismatch between the
-metal and the liquid phase is effectively eliminated by proper capping
-agent. Furthermore, the coverage precentage of the capping agent plays
-an important role in the interfacial thermal transport process.
+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 measurement results, particularly of the UA models, agree with
-available experimental data. This indicates that our force field
-parameters have a nice description of the interactions between the
-particles at the interfaces. AA models tend to overestimate the
+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
-vibration would be overly sampled. Compared to the AA models, the UA
-models have higher computational efficiency with satisfactory
-accuracy, and thus are preferable in interfacial thermal transport
-modelings.
+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.
 
-Vlugt {\it et al.} has investigated the surface thiol structures for
-nanocrystal gold and pointed out that they differs from those of the
-Au(111) surface\cite{vlugt:cpc2007154}. This difference might lead to
-change of interfacial thermal transport behavior as well. To
-investigate this problem, an effective means to introduce thermal flux
-and measure the corresponding thermal gradient is desirable for
-simulating structures with spherical symmetry.
+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}
