--- interfacial/interfacial.tex	2011/07/11 18:19:26	3734
+++ interfacial/interfacial.tex	2011/07/25 19:11:33	3749
@@ -23,8 +23,8 @@
 \setlength{\belowcaptionskip}{30 pt}
 
 %\renewcommand\citemid{\ } % no comma in optional reference note
-\bibpunct{[}{]}{,}{s}{}{;}
-\bibliographystyle{aip}
+\bibpunct{[}{]}{,}{n}{}{;}
+\bibliographystyle{achemso}
 
 \begin{document}
 
@@ -74,13 +74,15 @@ experimentally and computationally, due to its importa
 
 \section{Introduction}
 Interfacial thermal conductance is extensively studied both
-experimentally and computationally, 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
+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.
+which challenges traditional research methods for homogeneous
+systems.
 
 Heat conductance of molecular and nano-scale interfaces will be
 affected by the chemical details of the surface. Experimentally,
@@ -97,20 +99,20 @@ Theoretical and computational studies were also engage
 that specific ligands (capping agents) could completely eliminate this
 barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}.
 
-Theoretical and computational studies were also engaged in the
-interfacial thermal transport research in order to gain an
-understanding of this phenomena at the molecular level. 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 alkylthiolate with relatively long chain
-(8-20 carbon atoms)[CITE TWO PAPERS]. However, emsemble average measurements for heat
-conductance of interfaces between the capping monolayer on Au and a
-solvent phase has yet to be studied. The relatively 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
+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 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.
 
@@ -119,63 +121,65 @@ the two slabs. Furthermore, it allows more effective t
 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 under steady states.
+the two slabs. Furthermore, it allows effective thermal exchange
+between particles of different identities, and thus makes the study of
+interfacial conductance much simpler.
 
-Our work presented here investigated the Au(111) surface with various
-coverage of butanethiol, a capping agent with shorter carbon chain,
-solvated with organic solvents of different molecular shapes. And
-different models were used for both the capping agent and the solvent
-force field parameters. With the NIVS algorithm applied, the thermal
-transport through these interfacial systems was studied and the
-underlying mechanism for this phenomena was investigated.
+The work presented here deals with the Au(111) surface covered to
+varying degrees by butanethiol, a capping agent with short carbon
+chain, and solvated with organic solvents of different molecular
+properties. Different models were used for both the capping agent and
+the solvent force field parameters. Using the NIVS algorithm, the
+thermal transport across these interfaces was studied and the
+underlying mechanism for the phenomena was investigated.
 
-[WHY STUDY AU-THIOL SURFACE; MAY CITE SHAOYI JIANG]
+[MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]
 
 \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
+\subsection{Imposd-Flux Methods in MD Simulations}
+Steady state MD simulations has the advantage that not many
+trajectories are needed to study the relationship between thermal flux
+and thermal gradients. For systems including low conductance
+interfaces 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 those slowly-converging equilibrium
+methods\cite{Viscardy:2007lq}. Forward methods may impose gradient,
+but in interfacial conditions it is not clear what behavior to impose
+at the interfacial boundaries. Imposed-flux reverse 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.
+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.
 
-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 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
-the flux quantity is not restricted by particle mass difference.
+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
 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:
+Given a system with thermal gradients and the corresponding thermal
+flux, 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:
 \begin{equation}
 G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
     \langle T_\mathrm{cold}\rangle \right)}
@@ -186,18 +190,33 @@ When the interfacial conductance is {\it not} small, t
   T_\mathrm{cold}\rangle}$ are the average observed temperature of the
 two separated phases.
 
-When the interfacial conductance is {\it not} small, two ways can be
-used to define $G$.
+When the interfacial conductance is {\it not} small, there are two
+ways 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:
+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
+(Figure \ref{demoPic}):
 \begin{equation}
 G=\frac{J}{\Delta T}
 \label{discreteG}
 \end{equation}
 
+\begin{figure}
+\includegraphics[width=\linewidth]{method}
+\caption{Interfacial conductance can be calculated by applying an
+  (unphysical) kinetic energy flux between two slabs, one located
+  within the metal and another on the edge of the periodic box.  The
+  system responds by forming a thermal 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.}
+\label{demoPic}
+\end{figure}
+
 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
@@ -213,77 +232,70 @@ difference method and thus calculate $G^\prime$.
 
 With the temperature profile obtained from simulations, one is able to
 approximate the first and second derivatives of $T$ with finite
-difference method and thus calculate $G^\prime$.
+difference methods and thus calculate $G^\prime$.
 
-In what follows, both definitions are used for calculation and comparison.
+In what follows, both definitions have been used for calculation and
+are compared in the results.
 
-[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}.
+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 without capping agents on the
+surfaces, the metal slab is solvated with simple organic solvents, as
+illustrated in Figure \ref{gradT}.
 
-\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 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 of how an applied thermal flux can be used to obtain the 1st
+and 2nd derivatives of the temperature profile.
 
-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-butanethiol/hexane interfacial system and the
+  temperature profile after a kinetic energy flux is imposed to
+  it. The 1st and 2nd derivatives of the temperature profile can be
+  obtained with finite difference approximation (lower panel).}
 \label{gradT}
 \end{figure}
 
-[MAY INCLUDE POWER SPECTRUM PROTOCOL]
-
 \section{Computational Details}
 \subsection{Simulation Protocol}
-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.
-
-[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 NIVS algorithm has been implemented in our MD simulation code,
+OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
+simulations. Different metal slab thickness (layer numbers of Au) was
+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
+hollow sites on the Au(111) surfaces. These sites could be either a
+{\it fcc} or {\it hcp} site. However, Hase {\it et al.} found that
+they are equivalent in a heat transfer process\cite{hase:2010}, so
+they are not distinguished 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 different coverages was derived by evenly eliminating
+butanethiols on the surfaces, and was investigated in order to study
+the relation between coverage and interfacial conductance.
 
-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). [MAY
-EXPLAIN WHY WE CHOOSE THEM]
+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
+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 SNAPSHOTS]
 
-The spacing filled by solvent molecules, i.e. the gap between
+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
+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).
+
+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
@@ -292,26 +304,32 @@ corresponding spacing is usually $35 \sim 60$\AA.
 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.
+corresponding spacing is usually $35 \sim 75$\AA.
 
-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 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.
+To investigate this effect, comparisons were made with simulations
+that allow changes of $L_x$ and $L_y$ during NPT equilibration, and
+the results are shown in later sections. After ensuring the liquid
+phase reaches equilibrium at atmospheric pressure (1 atm), further
+equilibration are followed under NVT and then 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
 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:
+freeze due to excessively low temperature. After this induced
+temperature gradient is stablized, the temperature profile of the
+simulation cell is recorded. To do this, the simulation cell is
+devided evenly into $N$ slabs along the $z$-axis and $N$ is maximized
+for highest possible spatial resolution but not too many to have some
+slabs empty most of the time. 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
@@ -322,44 +340,74 @@ 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. And
+each equilibration / stabilization step usually takes 100 ps, or
+longer, if necessary.
+
 \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 various components. Figure \ref{demoMol}
+demonstrates the sites defined for both United-Atom and All-Atom
+models of the organic solvent and capping agent molecules in our
+simulations. Force field parameter descriptions 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
+  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}.}
+\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
+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] demonstrates how we name our pseudo-atoms of the
-molecules in our simulations.
-[FIGURE FOR MOLECULE NOMENCLATURE] 
-
 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
 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 not separated by more than 3 bonds. Otherwise, for non-bonded
+interactions, Lennard-Jones potentials are used. [CHECK CITATION]
 
+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 of the system, 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.
+
+For UA-toluene model, the non-bonded potentials between
+inter-molecular sites have a similar Lennard-Jones formulation. For
+intra-molecular interactions, considering the stiffness of the benzene
+ring, rigid body constraints are applied for further computational
+efficiency. All bonds in the benzene ring and between the ring and the
+methyl group remain rigid during the progress of simulations.
+
 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]
+force field is used. 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, the United Force Field developed by
+Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} is
+adopted. Without the rigid body constraints, bonding interactions were
+included. For the aromatic ring, improper torsions (inversions) were
+added as an extra potential for maintaining the planar shape.
+[CHECK CITATION]
 
 The capping agent in our simulations, the butanethiol molecules can
 either use UA or AA model. The TraPPE-UA force fields includes
@@ -368,60 +416,75 @@ Au(111) surfaces, we adopt the S parameters from [CITA
 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 [CITATION CF VLUGT]
-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: 
+Au(111) surfaces, we adopt the S parameters from Luedtke and
+Landman\cite{landman:1998}[CHECK CITATION]
+ 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:
+\begin{eqnarray}
+\sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
+\epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
+\end{eqnarray}
 
-
 To describe the interactions between metal Au and non-metal 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[CITATION] 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
+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.
 
 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]
+interactions between Au and aromatic C/H atoms in toluene. A set of
+pseudo Lennard-Jones parameters were provided for Au in their force
+fields. By using the Mixing Rule, this can be used to derive pair-wise
+potentials for non-bonded interactions between Au and non-metal sites.
 
 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''
+particles, such as All-Atom normal alkanes in our simulations are not
+yet well-established. For these interactions, we attempt to derive
+their parameters using the Mixing Rule. To do this, Au pseudo
+Lennard-Jones parameters for ``Metal-non-Metal'' (MnM) interactions
+were 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}
-      \caption{Lennard-Jones parameters for Au-non-Metal
-        interactions in our simulations.}
-      
-      \begin{tabular}{ccc}
+      \caption{Non-bonded interaction parameters (including cross
+        interactions with Au atoms) for both force fields used in this
+        work.}      
+      \begin{tabular}{lllllll}
         \hline\hline
-        Non-metal atom   & $\sigma$ & $\epsilon$ \\
-        (or pseudo-atom) & \AA      & kcal/mol  \\
+        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
+        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
+        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
         \hline
-        S    & 2.40   & 8.465   \\
-        CH3  & 3.54   & 0.2146  \\
-        CH2  & 3.54   & 0.1749  \\
-        CT3  & 3.365  & 0.1373  \\
-        CT2  & 3.365  & 0.1373  \\
-        CTT  & 3.365  & 0.1373  \\
-        HC   & 2.865  & 0.09256 \\
-        CHar & 3.4625 & 0.1680  \\
-        CRar & 3.555  & 0.1604  \\
-        CA   & 3.173  & 0.0640  \\
-        HA   & 2.746  & 0.0414  \\
+        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{MnM}
@@ -429,11 +492,35 @@ parameters in our simulations.
   \end{minipage}
 \end{table*}
 
+\subsection{Vibrational Spectrum}
+To investigate the mechanism of interfacial thermal conductance, the
+vibrational spectrum is utilized as a complementary tool. Vibrational
+spectra were taken for individual components in different
+simulations. To obtain these spectra, simulations were run after
+equilibration, in the NVE ensemble. Snapshots of configurations were
+collected at a frequency that is higher than that of the fastest
+vibrations occuring 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}
+Followed by Fourier transforms, the power spectrum can be constructed:
+\begin{equation}
+\hat{f}(\omega) = \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
+\label{fourier}
+\end{equation}
 
 \section{Results and Discussions}
-[MAY HAVE A BRIEF SUMMARY]
+In what follows, how the parameters and protocol of simulations would
+affect the measurement of $G$'s is first discussed. With a reliable
+protocol and set of parameters, the influence of capping agent
+coverage on thermal conductance is investigated. Besides, different
+force field models for both solvents and selected deuterated models
+were tested and compared. Finally, a summary of the role of capping
+agent in the interfacial thermal transport process is given.
+
 \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
@@ -442,14 +529,15 @@ during equilibrating the liquid phase. Due to the stif
 results.
 
 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. 
+during equilibrating the liquid phase. Due to the stiffness of the
+crystalline Au structure, $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 of $L_x$ and $L_y$ (which is significantly smaller),
+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.
 
 As stated in our computational details, the spacing filled with
 solvent molecules can be chosen within a range. This allows some
@@ -476,44 +564,47 @@ $J_z$ is listed in Table \ref{AuThiolHexaneUA} and [RE
 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.
+$J_z$ is listed in Table \ref{AuThiolHexaneUA} and
+\ref{AuThiolToluene}. 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.}
+        at different temperatures using a range of energy
+        fluxes. Error estimates indicated in parenthesis.}
       
-      \begin{tabular}{cccccccc}
+      \begin{tabular}{ccccccc}
         \hline\hline
-        $\langle T\rangle$ & & $L_x$ & $L_y$ & $L_z$ & $J_z$ &
-        $G$ & $G^\prime$ \\
-        (K) & $N_{hexane}$ & \multicolumn{3}{c}{(\AA)} & (GW/m$^2$) &
+        $\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 & 29.86 & 25.80 & 113.1 & -0.96 &
-        102()  & 80.0() \\
-            & 200 & 29.84 & 25.81 &  93.9 &  1.92 &
-        129()  & 87.3() \\
-            &     & 29.84 & 25.81 &  95.3 &  1.93 &
-        131()  & 77.5() \\
-            & 166 & 29.84 & 25.81 &  85.7 &  0.97 &
-        115()  & 69.3() \\
-            &     &       &       &       &  1.94 &
-        125()  & 87.1() \\
-        250 & 200 & 29.84 & 25.87 & 106.8 &  0.96 &
-        81.8() & 67.0() \\
-            & 166 & 29.87 & 25.84 &  94.8 &  0.98 &
-        79.0() & 62.9() \\
-            &     & 29.84 & 25.85 &  95.0 &  1.44 &
-        76.2() & 64.8() \\
+        200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
+            & 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
+            &     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
+            &     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
+            &     &     &       &  1.91 & 139(10)   & 101(10)   \\
+            &     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
+            & 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
+            &     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
+            &     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
+            &     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
+        \hline
+        250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
+            &     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
+            & 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
+            &     & No  & 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{AuThiolHexaneUA}
@@ -529,7 +620,8 @@ butanethiol as well.[MAY NEED FIGURE] And this reduced
 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
+butanethiol as well.[MAY NEED SLAB DENSITY FIGURE] 
+And this reduced contact would
 probably be accountable for a lower interfacial thermal conductance,
 as shown in Table \ref{AuThiolHexaneUA}. 
 
@@ -544,25 +636,25 @@ in that higher degree of contact could yield increased
 important role in the thermal transport process across the interface
 in that higher degree of contact could yield increased conductance.
 
-[ADD Lxyz AND 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.}
+        fluxes. Error estimates indicated in parenthesis.}
       
-      \begin{tabular}{cccc}
+      \begin{tabular}{ccccc}
         \hline\hline
-        $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
-        (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        $\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 & -1.86 & 180() & 135() \\
-            &  2.15 & 204() & 113() \\
-            & -3.93 & 175() & 114() \\
-        300 & -1.91 & 143() & 125() \\
-            & -4.19 & 134() & 113() \\
+        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}
@@ -594,40 +686,50 @@ above. Our Au-butanethiol/toluene system did not see t
 
 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.
+above. 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. The empty three-fold sites not occupied by
+capping agents 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.
 
 \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.
+thermal conductivity are measured with our NIVS algorithm. Figure
+\ref{coverage} demonstrates the trend of conductance change with
+respect to different coverages of butanethiol. To study the isotope
+effect in interfacial thermal conductance, deuterated UA-hexane is
+included as well.
 
+\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}
+
 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. 
+surface, the derivative definition for $G^\prime$
+(Eq. \ref{derivativeG}) was difficult 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 the Gibbs
+deviding surface can still be well-defined. Therefore, $G$ (not
+$G^\prime$) was used for this section.
 
-From Table \ref{tlnUhxnUhxnD}, one can see the significance of the
+From Figure \ref{coverage}, 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.
+phenomena on metal / organic solvent surfaces.
 
 Interestingly, as one could observe from our results, the maximum
 conductance enhancement (largest $G$) happens while the surfaces are
@@ -646,7 +748,8 @@ case, $G$ decrease could not be offset but instead acc
 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.
+case, $G$ decrease could not be offset but instead accelerated. [NEED
+SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
 
 A comparison of the results obtained from differenet organic solvents
 can also provide useful information of the interfacial thermal
@@ -656,19 +759,29 @@ difference for the results of $G$. [MAY NEED FIGURE]
 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]
+difference for the results of $G$ (Figure \ref{uahxnua}).
 
+\begin{figure}
+\includegraphics[width=\linewidth]{uahxnua}
+\caption{Vibrational spectra obtained for normal (upper) and
+  deuterated (lower) hexane in Au-butanethiol/hexane
+  systems. Butanethiol spectra are shown as reference. Both hexane and
+  butanethiol were using United-Atom models.}
+\label{uahxnua}
+\end{figure}
+
 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.
+ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
+suggests that explicit hydrogen might not be a required factor for
+modeling thermal transport phenomena of systems such as
+Au-thiol/organic solvent.
 
 However, results for Au-butanethiol/toluene do not show an identical
-trend with those for Au-butanethiol/hexane in that $G$'s remain at
+trend with those for Au-butanethiol/hexane in that $G$ remains 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
+for planar 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
@@ -676,37 +789,8 @@ can see a plateau of $G$ vs. butanethiol coverage in o
 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, MAY ALSO PUT J HERE]
-\begin{table*}
-  \begin{minipage}{\linewidth}
-    \begin{center}
-      \caption{Computed 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.}
-      
-      \begin{tabular}{cccc}
-        \hline\hline
-        Thiol & \multicolumn{3}{c}{$G$(MW/m$^2$/K)} \\
-        coverage (\%) & hexane & hexane(D) & toluene \\
-        \hline
-        0.0   & 46.5() & 43.9() & 70.1() \\
-        25.0  & 151()  & 153()  & 249()  \\
-        50.0  & 172()  & 182()  & 214()  \\
-        75.0  & 242()  & 229()  & 244()  \\
-        88.9  & 178()  & -      & -      \\
-        100.0 & 137()  & 153()  & 187()  \\
-        \hline\hline
-      \end{tabular}
-      \label{tlnUhxnUhxnD}
-    \end{center}
-  \end{minipage}
-\end{table*}
 
 \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
@@ -715,7 +799,6 @@ these studies.
 the previous section. Table \ref{modelTest} summarizes the results of
 these studies.
 
-[MORE DATA; ERROR ESTIMATE]
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
@@ -723,23 +806,37 @@ these studies.
       \caption{Computed interfacial thermal conductivity ($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.}
+        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
+        or capping agent molecules; ``Avg.'' denotes results that are
+        averages of simulations under different $J_z$'s. Error
+        estimates indicated in parenthesis.)}
       
-      \begin{tabular}{ccccc}
+      \begin{tabular}{llccc}
         \hline\hline
         Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
         (or bare surface) & model & (GW/m$^2$) &
         \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        UA    & AA hexane  & 1.94 & 135()  & 129()  \\
-              &            & 2.86 & 126()  & 115()  \\
-              & AA toluene & 1.89 & 200()  & 149()  \\
-        AA    & UA hexane  & 1.94 & 116()  & 129()  \\
-              & AA hexane  & 3.76 & 451()  & 378()  \\
-              &            & 4.71 & 432()  & 334()  \\
-              & AA toluene & 3.79 & 487()  & 290()  \\
-        AA(D) & UA hexane  & 1.94 & 158()  & 172()  \\
-        bare  & AA hexane  & 0.96 & 31.0() & 29.4() \\
+        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{modelTest}
@@ -762,71 +859,81 @@ the upper bond of experimental value range. This is pr
 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.
+the general range of experimental measurement results. 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 (Figure
+\ref{aahxntln}). 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.
 
+\begin{figure}
+\includegraphics[width=\linewidth]{aahxntln}
+\caption{Spectra obtained for All-Atom model Au-butanethil/solvent
+  systems. When butanethiol is deuterated (lower left), its
+  vibrational overlap with hexane would decrease significantly,
+  compared with normal butanethiol (upper left). However, this
+  dramatic change does not apply to toluene as much (right).}
+\label{aahxntln}
+\end{figure}
+
 However, for Au-butanethiol/toluene interfaces, having the AA
 butanethiol deuterated did not yield a significant change in the
-measurement results.
-. , 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.
+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 heat exchange
+process. 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.
 
-
 Although the QSC model for Au is known to predict an overly low value
-for bulk metal gold conductivity[CITE NIVSRNEMD], our computational
+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, the modeling of interfacial
-thermal transport behavior relies mainly on an accurate description of
-the interactions between components occupying the interfaces.
+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}
+\subsection{Role of Capping Agent in Interfacial Thermal Conductance}
+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 covered by butanethiol molecules, compared
+to bare gold surfaces, exhibit an additional peak observed at the
+frequency of $\sim$170cm$^{-1}$, which is attributed to the S-Au
+bonding vibration. This vibration enables efficient thermal transport
+from surface Au layer to the capping agents. Therefore, in our
+simulations, the Au/S interfaces do not appear major heat barriers
+compared to the butanethiol / solvent interfaces.
 
-[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL, EQN'S]
+Simultaneously, the vibrational overlap between butanethiol and
+organic solvents suggests higher thermal exchange efficiency between
+these two components. Even exessively high heat transport was observed
+when All-Atom models were used and C-H vibrations were treated
+classically. Compared to metal and organic liquid phase, the heat
+transfer efficiency between butanethiol and organic solvents is closer
+to that within bulk liquid phase.
 
-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.
+Furthermore, our observation validated previous
+results\cite{hase:2010} that the intramolecular heat transport of
+alkylthiols is highly effecient. As a combinational effects of these
+phenomena, butanethiol acts as a channel to expedite thermal transport
+process. The acoustic impedance mismatch between the metal and the
+liquid phase can be effectively reduced with the presence of suitable
+capping agents.
 
- 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.
-
-[MAY NEED TO CONVERT TO JPEG]
 \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}
+  environments.}
+\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.
+[MAY ADD COMPARISON OF AU SLAB WIDTHS]
 
 \section{Conclusions}
 The NIVS algorithm we developed has been applied to simulations of
@@ -834,12 +941,20 @@ conductivities. Our simulations have seen significant 
 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
+conductivities. Under steady states, single trajectory simulation
+would be enough for accurate measurement. This would be advantageous
+compared to transient state simulations, which need multiple
+trajectories to produce reliable average results.
+
+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.
+an important role in the interfacial thermal transport
+process. Moderately lower coverages allow higher contact between
+capping agent and solvent, and thus could further enhance the heat
+transfer process.
 
 Our measurement results, particularly of the UA models, agree with
 available experimental data. This indicates that our force field
@@ -849,7 +964,13 @@ modelings.
 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.
+modelings. 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 is
+limited for accurate computation of derivatives data.
 
 Vlugt {\it et al.} has investigated the surface thiol structures for
 nanocrystal gold and pointed out that they differs from those of the
@@ -859,7 +980,6 @@ simulating structures with spherical symmetry.
 and measure the corresponding thermal gradient is desirable for
 simulating structures with spherical symmetry.
 
-
 \section{Acknowledgments}
 Support for this project was provided by the National Science
 Foundation under grant CHE-0848243. Computational time was provided by
