--- interfacial/interfacial.tex	2011/07/25 19:11:33	3749
+++ interfacial/interfacial.tex	2011/07/26 23:01:51	3752
@@ -73,31 +73,28 @@ 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 in nanotechnology, interfacial
+thermal conductance has been studied extensively both experimentally
+and computationally.\cite{cahill:793} Unlike bulk materials, nanoscale
+materials have a significant fraction of their atoms at interfaces,
+and the chemical details of these interfaces govern the heat transfer
+behavior. Furthermore, the interfaces are
+heterogeneous (e.g. solid - liquid), which provides a challenge to
+traditional methods developed for homogeneous 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
+Experimentally, various interfaces have been investigated for their
+thermal conductance. 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} and 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}.
+nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
+normally considered barriers for heat transport, Alper {\it et al.}
+suggested that specific ligands (capping agents) could completely
+eliminate this barrier
+($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}
 
 Theoretical and computational models have also been used to study the
 interfacial thermal transport in order to gain an understanding of
@@ -105,18 +102,22 @@ atoms)\cite{hase:2010,hase:2011}. However, ensemble av
 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
+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 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.
+methods. Therefore, the Reverse NEMD (RNEMD)
+methods\cite{MullerPlathe:1997xw,kuang:164101} would have the
+advantage of applying this difficult to measure flux (while measuring
+the resulting gradient), given that the simulation methods being able
+to effectively apply an unphysical flux in non-homogeneous systems.
+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 chemistry details of
+interfaces, e.g. hydrophobic and hydrophilic aqueous interfaces.
 
-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
@@ -133,30 +134,28 @@ underlying mechanism for the phenomena was investigate
 thermal transport across these interfaces was studied and the
 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}
-Steady state MD simulations has the advantage that not many
+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 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 thermal response becomes easier to measure than the flux. Although
+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
@@ -173,33 +172,31 @@ 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$}
-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:
+\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.
 
 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
-(Figure \ref{demoPic}):
+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 deviding surface (Figure \ref{demoPic}):
 \begin{equation}
-G=\frac{J}{\Delta T}
+  G=\frac{J}{\Delta T}
 \label{discreteG}
 \end{equation}
 
@@ -219,7 +216,7 @@ the magnitude of thermal conductivity $\lambda$ change
 
 The other approach is to assume a continuous temperature profile along
 the thermal gradient axis (e.g. $z$) and define $G$ at the point where
-the magnitude of thermal conductivity $\lambda$ change reach its
+the magnitude of thermal conductivity ($\lambda$) change reaches its
 maximum, given that $\lambda$ is well-defined throughout the space:
 \begin{equation}
 G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
@@ -230,27 +227,25 @@ With the temperature profile obtained from simulations
 \label{derivativeG}
 \end{equation}
 
-With the temperature profile obtained from simulations, one is able to
+With temperature profiles obtained from simulation, one is able to
 approximate the first and second derivatives of $T$ with finite
-difference methods and thus calculate $G^\prime$.
+difference methods and calculate $G^\prime$. In what follows, both
+definitions have been used, and are compared in the results.
 
-In what follows, both definitions 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 without capping agents on the
-surfaces, the metal slab is solvated with simple organic solvents, as
+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 of how an applied thermal flux 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}
@@ -264,74 +259,71 @@ OpenMD\cite{Meineke:2005gd,openmd}, and was used for o
 \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 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
+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 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.
+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 SNAPSHOTS]
+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 75$\AA.
+these extreme cases did not happen to our simulations. The spacing
+between periodic images of the gold interfaces is $45 \sim 75$\AA.
 
 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.
+$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. 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
+freeze due to lowered temperatures. After this induced temperature
+gradient had stablized, the temperature profile of the simulation cell
+was recorded. To do this, the simulation cell is devided 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|
+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
@@ -340,15 +332,15 @@ All of the above simulation procedures use a time step
 \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.
+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. 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
+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.
 
@@ -358,104 +350,90 @@ particles of different species.
   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}.}
+  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes,TraPPE-UA.thiols} (UA) and \protect\cite{OPLSAA} (AA). Cross-interactions with the Au atoms are given in Table \ref{MnM}.}
 \label{demoMol}
 \end{figure}
 
 The Au-Au interactions in metal lattice slab is described by the
-quantum Sutton-Chen (QSC) formulation\cite{PhysRevB.59.3527}. The QSC
+quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
 potentials include zero-point quantum corrections and are
 reparametrized for accurate surface energies compared to the
-Sutton-Chen potentials\cite{Chen90}. 
+Sutton-Chen potentials.\cite{Chen90}
 
-For both solvent molecules, straight chain {\it n}-hexane and aromatic
-toluene, United-Atom (UA) and All-Atom (AA) models are used
-respectively. The TraPPE-UA
+For the two solvent molecules, {\it n}-hexane and toluene, two
+different atomistic models were utilized. Both solvents were modeled
+using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
 parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
 for our UA solvent molecules. In these models, 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]
+sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones
+potentials are used.
 
-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.
+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 slighly
+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 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.
+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. Additional explicit hydrogen sites were
+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, 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]
+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}[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:
+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. 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, 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.
+The potentials developed from {\it ab initio} calculations by Leng
+{\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
+interactions between Au and aromatic C/H atoms in toluene. However,
+the Lennard-Jones parameters between Au and other types of particles,
+(e.g. AA alkanes) have not yet been established. For these
+interactions, the Lorentz-Berthelot mixing rule can be used to derive
+effective single-atom LJ parameters for the metal using the fit values
+for toluene. These are then used to construct reasonable mixing
+parameters for the interactions between the gold and other atoms.
+Table \ref{MnM} summarizes the ``metal/non-metal'' parameters used in
+our simulations.
 
 \begin{table*}
   \begin{minipage}{\linewidth}
@@ -492,23 +470,26 @@ parameters in our simulations.
   \end{minipage}
 \end{table*}
 
-\subsection{Vibrational Spectrum}
+\subsection{Vibrational Power 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
+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 NVE
+ensemble, and without a thermal gradient. 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:
+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:
+\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}
+  \hat{f}(\omega) =
+  \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
+\label{fourier} 
 \end{equation}
 
 \section{Results and Discussions}
@@ -748,7 +729,7 @@ 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. [NEED
+case, $G$ decrease could not be offset but instead accelerated. [MAY NEED
 SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
 
 A comparison of the results obtained from differenet organic solvents
@@ -974,11 +955,11 @@ Au(111) surface\cite{vlugt:cpc2007154}. This differenc
 
 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.
+Au(111) surface\cite{landman:1998,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.
 
 \section{Acknowledgments}
 Support for this project was provided by the National Science
