--- interfacial/interfacial.tex	2011/07/27 03:27:28	3753
+++ interfacial/interfacial.tex	2011/07/29 21:06:30	3759
@@ -83,7 +83,12 @@ thermal conductance. Wang {\it et al.} studied heat tr
 traditional methods developed for homogeneous systems. 
 
 Experimentally, various interfaces have been investigated for their
-thermal conductance. Wang {\it et al.} studied heat transport through
+thermal conductance. Cahill and coworkers studied nanoscale thermal
+transport from metal nanoparticle/fluid interfaces, to epitaxial
+TiN/single crystal oxides interfaces, to hydrophilic and hydrophobic
+interfaces between water and solids with different self-assembled
+monolayers.\cite{Wilson:2002uq,PhysRevB.67.054302,doi:10.1021/jp048375k,PhysRevLett.96.186101}
+Wang {\it et al.} studied heat transport through
 long-chain hydrocarbon monolayers on gold substrate at individual
 molecular level,\cite{Wang10082007} Schmidt {\it et al.} studied the
 role of CTAB on thermal transport between gold nanorods and
@@ -106,7 +111,8 @@ difficult to measure with Equilibrium MD or forward NE
 measurements for heat conductance of interfaces between the capping
 monolayer on Au and a solvent phase have yet to be studied with their
 approach. The comparatively low thermal flux through interfaces is
-difficult to measure with Equilibrium MD or forward NEMD simulation
+difficult to measure with Equilibrium
+MD\cite{doi:10.1080/0026897031000068578} or forward NEMD simulation
 methods. Therefore, the Reverse NEMD (RNEMD)
 methods\cite{MullerPlathe:1997xw,kuang:164101} would have the
 advantage of applying this difficult to measure flux (while measuring
@@ -187,14 +193,18 @@ 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.
+temperature of the two separated phases.  For an applied flux $J_z$
+operating over a simulation time $t$ on a periodically-replicated slab
+of dimensions $L_x \times L_y$, $E_{total} = J_z *(t)*(2 L_x L_y)$.
 
 When the interfacial conductance is {\it not} small, there are two
 ways to define $G$. One 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}):
+occurs at the Gibbs deviding surface (Figure \ref{demoPic}). This is
+known as the Kapitza conductance, which is the inverse of the Kapitza
+resistance.
 \begin{equation}
   G=\frac{J}{\Delta T}
 \label{discreteG}
@@ -350,7 +360,9 @@ 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,TraPPE-UA.thiols} (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}
 
@@ -470,88 +482,229 @@ our simulations.
   \end{minipage}
 \end{table*}
 
-\subsection{Vibrational Power Spectrum}
 
-To investigate the mechanism of interfacial thermal conductance, the
-vibrational power spectrum was computed. Power spectra were taken for
-individual components in different simulations. To obtain these
-spectra, simulations were run after equilibration, in the 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: 
-\begin{equation}
-C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
-\label{vCorr} 
-\end{equation}
-The power spectrum is constructed via a Fourier transform of the
-symmetrized velocity autocorrelation function,
-\begin{equation}
-  \hat{f}(\omega) =
-  \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
-\label{fourier} 
-\end{equation}
-
-\section{Results and Discussions}
-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.
+\section{Results}
+There are many factors contributing to the measured interfacial
+conductance; some of these factors are physically motivated
+(e.g. coverage of the surface by the capping agent coverage and
+solvent identity), while some are governed by parameters of the
+methodology (e.g. applied flux and the formulas used to obtain the
+conductance). In this section we discuss the major physical and
+calculational effects on the computed conductivity.
 
-\subsection{How Simulation Parameters Affects $G$}
-We have varied our protocol or other parameters of the simulations in
-order to investigate how these factors would affect the measurement of
-$G$'s. It turned out that while some of these parameters would not
-affect the results substantially, some other changes to the
-simulations would have a significant impact on the measurement
-results.
+\subsection{Effects due to capping agent coverage}
 
-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
-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.
+A series of different initial conditions with a range of surface
+coverages was prepared and solvated with various with both of the
+solvent molecules. These systems were then equilibrated and their
+interfacial thermal conductivity was measured with the NIVS
+algorithm. Figure \ref{coverage} demonstrates the trend of conductance
+with respect to surface coverage. 
 
-As stated in our computational details, the spacing filled with
-solvent molecules can be chosen within a range. This allows some
-change of solvent molecule numbers for the same Au-butanethiol
-surfaces. We did this study on our Au-butanethiol/hexane
-simulations. Nevertheless, the results obtained from systems of
-different $N_{hexane}$ did not indicate that the measurement of $G$ is
-susceptible to this parameter. For computational efficiency concern,
-smaller system size would be preferable, given that the liquid phase
-structure is not affected.
+\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.}
+\label{coverage}
+\end{figure}
 
-Our NIVS algorithm allows change of unphysical thermal flux both in
-direction and in quantity. This feature extends our investigation of
-interfacial thermal conductance. However, the magnitude of this
-thermal flux is not arbitary if one aims to obtain a stable and
-reliable thermal gradient. A temperature profile would be
-substantially affected by noise when $|J_z|$ has a much too low
-magnitude; while an excessively large $|J_z|$ that overwhelms the
-conductance capacity of the interface would prevent a thermal gradient
-to reach a stablized steady state. NIVS has the advantage of allowing
-$J$ to vary in a wide range such that the optimal flux range for $G$
-measurement can generally be simulated by the algorithm. Within the
-optimal range, we were able to study how $G$ would change according to
-the thermal flux across the interface. For our simulations, we denote
-$J_z$ to be positive when the physical thermal flux is from the liquid
-to metal, and negative vice versa. The $G$'s measured under different
-$J_z$ is listed in Table \ref{AuThiolHexaneUA} and
-\ref{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.
+In partially covered surfaces, the derivative definition for
+$G^\prime$ (Eq. \ref{derivativeG}) becomes difficult to apply, as the
+location of maximum change of $\lambda$ becomes washed out.  The
+discrete definition (Eq. \ref{discreteG}) is easier to apply, as the
+Gibbs dividing surface is still well-defined. Therefore, $G$ (not
+$G^\prime$) was used in this section.
+
+From Figure \ref{coverage}, one can see the significance of the
+presence of capping agents. When even a small fraction of the Au(111)
+surface sites are covered with butanethiols, the conductivity exhibits
+an enhancement by at least a factor of 3.  Cappping agents are clearly
+playing a major role in thermal transport at metal / organic solvent
+surfaces.
+
+We note a non-monotonic behavior in the interfacial conductance as a
+function of surface coverage. The maximum conductance (largest $G$)
+happens when the surfaces are about 75\% covered with butanethiol
+caps.  The reason for this behavior is not entirely clear.  One
+explanation is that incomplete butanethiol coverage allows small gaps
+between butanethiols to form. These gaps can be filled by transient
+solvent molecules.  These solvent molecules couple very strongly with
+the hot capping agent molecules near the surface, and can then carry
+away (diffusively) the excess thermal energy from the surface.
+
+There appears to be a competition between the conduction of the
+thermal energy away from the surface by the capping agents (enhanced
+by greater coverage) and the coupling of the capping agents with the
+solvent (enhanced by interdigitation at lower coverages).  This
+competition would lead to the non-monotonic coverage behavior observed
+here.
+
+Results for rigid body toluene solvent, as well as the UA hexane, are
+within the ranges expected from prior experimental
+work.\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406} This suggests
+that explicit hydrogen atoms might not be required for modeling
+thermal transport in these systems.  C-H vibrational modes do not see
+significant excited state population at low temperatures, and are not
+likely to carry lower frequency excitations from the solid layer into
+the bulk liquid.
+
+The toluene solvent does not exhibit the same behavior as hexane in
+that $G$ remains at approximately the same magnitude when the capping
+coverage increases from 25\% to 75\%.  Toluene, as a rigid planar
+molecule, cannot occupy the relatively small gaps between the capping
+agents as easily as the chain-like {\it n}-hexane.  The effect of
+solvent coupling to the capping agent is therefore weaker in toluene
+except at the very lowest coverage levels.  This effect counters the
+coverage-dependent conduction of heat away from the metal surface,
+leading to a much flatter $G$ vs. coverage trend than is observed in
+{\it n}-hexane.
+
+\subsection{Effects due to Solvent \& Solvent Models}
+In addition to UA solvent and capping agent models, AA models have
+also been included in our simulations.  In most of this work, the same
+(UA or AA) model for solvent and capping agent was used, but it is
+also possible to utilize different models for different components.
+We have also included isotopic substitutions (Hydrogen to Deuterium)
+to decrease the explicit vibrational overlap between solvent and
+capping agent. Table \ref{modelTest} summarizes the results of these
+studies.
+
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      
+      \caption{Computed interfacial thermal conductance ($G$ and
+        $G^\prime$) values for interfaces using various models for
+        solvent and capping agent (or without capping agent) at
+        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
+        or capping agent molecules; ``Avg.'' denotes results that are
+        averages of simulations under different applied thermal flux
+        values $(J_z)$. Error estimates are indicated in
+        parentheses.)}
+      
+      \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    & 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}
+    \end{center}
+  \end{minipage}
+\end{table*}
 
+To facilitate direct comparison between force fields, systems with the
+same capping agent and solvent were prepared with the same length
+scales for the simulation cells.
+
+On bare metal / solvent surfaces, different force field models for
+hexane yield similar results for both $G$ and $G^\prime$, and these
+two definitions agree with each other very well. This is primarily an
+indicator of weak interactions between the metal and the solvent, and
+is a typical case for acoustic impedance mismatch between these two
+phases.  
+
+For the fully-covered surfaces, the choice of force field for the
+capping agent and solvent has a large impact on the calulated values
+of $G$ and $G^\prime$.  The AA thiol to AA solvent conductivities are
+much larger than their UA to UA counterparts, and these values exceed
+the experimental estimates by a large measure.  The AA force field
+allows significant energy to go into C-H (or C-D) stretching modes,
+and since these modes are high frequency, this non-quantum behavior is
+likely responsible for the overestimate of the conductivity.  Compared
+to the AA model, the UA model yields more reasonable conductivity
+values with much higher computational efficiency.
+
+\subsubsection{Are electronic excitations in the metal important?}
+Because they lack electronic excitations, the QSC and related embedded
+atom method (EAM) models for gold are known to predict unreasonably
+low values for bulk conductivity
+($\lambda$).\cite{kuang:164101,ISI:000207079300006,Clancy:1992} If the
+conductance between the phases ($G$) is governed primarily by phonon
+excitation (and not electronic degrees of freedom), one would expect a
+classical model to capture most of the interfacial thermal
+conductance.  Our results for $G$ and $G^\prime$ indicate that this is
+indeed the case, and suggest that the modeling of interfacial thermal
+transport depends primarily on the description of the interactions
+between the various components at the interface.  When the metal is
+chemically capped, the primary barrier to thermal conductivity appears
+to be the interface between the capping agent and the surrounding
+solvent, so the excitations in the metal have little impact on the
+value of $G$.
+
+\subsection{Effects due to methodology and simulation parameters}
+
+We have varied the parameters of the simulations in order to
+investigate how these factors would affect the computation of $G$.  Of
+particular interest are: 1) the length scale for the applied thermal
+gradient (modified by increasing the amount of solvent in the system),
+2) the sign and magnitude of the applied thermal flux, 3) the average
+temperature of the simulation (which alters the solvent density during
+equilibration), and 4) the definition of the interfacial conductance
+(Eqs. (\ref{discreteG}) or (\ref{derivativeG})) used in the
+calculation.
+
+Systems of different lengths were prepared by altering the number of
+solvent molecules and extending the length of the box along the $z$
+axis to accomodate the extra solvent.  Equilibration at the same
+temperature and pressure conditions led to nearly identical surface
+areas ($L_x$ and $L_y$) available to the metal and capping agent,
+while the extra solvent served mainly to lengthen the axis that was
+used to apply the thermal flux.  For a given value of the applied
+flux, the different $z$ length scale has only a weak effect on the
+computed conductivities (Table \ref{AuThiolHexaneUA}).
+
+\subsubsection{Effects of applied flux}
+The NIVS algorithm allows changes in both the sign and magnitude of
+the applied flux.  It is possible to reverse the direction of heat
+flow simply by changing the sign of the flux, and thermal gradients
+which would be difficult to obtain experimentally ($5$ K/\AA) can be
+easily simulated.  However, the magnitude of the applied flux is not
+arbitary if one aims to obtain a stable and reliable thermal gradient.
+A temperature gradient can be lost in the noise if $|J_z|$ is too
+small, and excessive $|J_z|$ values can cause phase transitions if the
+extremes of the simulation cell become widely separated in
+temperature.  Also, if $|J_z|$ is too large for the bulk conductivity
+of the materials, the thermal gradient will never reach a stable
+state.  
+
+Within a reasonable range of $J_z$ values, we were able to study how
+$G$ changes as a function of this flux.  In what follows, we use
+positive $J_z$ values to denote the case where energy is being
+transferred by the method from the metal phase and into the liquid.
+The resulting gradient therefore has a higher temperature in the
+liquid phase.  Negative flux values reverse this transfer, and result
+in higher temperature metal phases.  The conductance measured under
+different applied $J_z$ values is listed in Tables
+\ref{AuThiolHexaneUA} and \ref{AuThiolToluene}. These results do not
+indicate that $G$ depends strongly on $J_z$ within this flux
+range. The linear response of flux to thermal gradient simplifies our
+investigations in that we can rely on $G$ measurement with only a
+small number $J_z$ values.  
+
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
@@ -563,29 +716,24 @@ a large series of fluxes.
       
       \begin{tabular}{ccccccc}
         \hline\hline
-        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
+        $\langle T\rangle$ & $N_{hexane}$  & $\rho_{hexane}$ &
         $J_z$ & $G$ & $G^\prime$ \\
-        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
+        (K) &  & (g/cm$^3$) & (GW/m$^2$) &
         \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        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) \\
+        200 & 266 &  0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
+            & 200 &  0.688 &  0.96 & 125(16)   & 90.2(15)  \\
+            &     &        &  1.91 & 139(10)   & 101(10)   \\
+            &     &        &  2.83 & 141(6)    & 89.9(9.8) \\
+            & 166 &  0.681 &  0.97 & 141(30)   & 78(22)    \\
+            &     &        &  1.92 & 138(4)    & 98.9(9.5) \\
         \hline
-        250 & 200 & 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) \\
+        250 & 200 &  0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
+            &     &        & -0.95 & 49.4(0.3) & 45.7(2.1) \\
+            & 166 &  0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
+            &     &        &  1.44 & 76.2(5.0) & 64.8(3.8) \\
+            &     &        & -0.95 & 56.4(2.5) & 54.4(1.1) \\
+            &     &        & -1.85 & 47.8(1.1) & 53.5(1.5) \\
         \hline\hline
       \end{tabular}
       \label{AuThiolHexaneUA}
@@ -593,30 +741,42 @@ Furthermore, we also attempted to increase system aver
   \end{minipage}
 \end{table*}
 
-Furthermore, we also attempted to increase system average temperatures
-to above 200K. These simulations are first equilibrated in the NPT
-ensemble under normal pressure. As stated above, the TraPPE-UA model
-for hexane tends to predict a lower boiling point. In our simulations,
-hexane had diffculty to remain in liquid phase when NPT equilibration
-temperature is higher than 250K. Additionally, the equilibrated liquid
-hexane density under 250K becomes lower than experimental value. This
-expanded liquid phase leads to lower contact between hexane and
-butanethiol as well.[MAY NEED SLAB DENSITY FIGURE] 
-And this reduced contact would
-probably be accountable for a lower interfacial thermal conductance,
-as shown in Table \ref{AuThiolHexaneUA}. 
+The sign of $J_z$ is a different matter, however, as this can alter
+the temperature on the two sides of the interface. The average
+temperature values reported are for the entire system, and not for the
+liquid phase, so at a given $\langle T \rangle$, the system with
+positive $J_z$ has a warmer liquid phase.  This means that if the
+liquid carries thermal energy via convective transport, {\it positive}
+$J_z$ values will result in increased molecular motion on the liquid
+side of the interface, and this will increase the measured
+conductivity.
 
-A similar study for TraPPE-UA toluene agrees with the above result as
-well. Having a higher boiling point, toluene tends to remain liquid in
-our simulations even equilibrated under 300K in NPT
-ensembles. Furthermore, the expansion of the toluene liquid phase is
-not as significant as that of the hexane. This prevents severe
-decrease of liquid-capping agent contact and the results (Table
-\ref{AuThiolToluene}) show only a slightly decreased interface
-conductance. Therefore, solvent-capping agent contact should play an
-important role in the thermal transport process across the interface
-in that higher degree of contact could yield increased conductance.
+\subsubsection{Effects due to average temperature}
 
+We also studied the effect of average system temperature on the
+interfacial conductance.  The simulations are first equilibrated in
+the NPT ensemble at 1 atm.  The TraPPE-UA model for hexane tends to
+predict a lower boiling point (and liquid state density) than
+experiments.  This lower-density liquid phase leads to reduced contact
+between the hexane and butanethiol, and this accounts for our
+observation of lower conductance at higher temperatures as shown in
+Table \ref{AuThiolHexaneUA}.  In raising the average temperature from
+200K to 250K, the density drop of ~20\% in the solvent phase leads to
+a ~65\% drop in the conductance.
+
+Similar behavior is observed in the TraPPE-UA model for toluene,
+although this model has better agreement with the experimental
+densities of toluene.  The expansion of the toluene liquid phase is
+not as significant as that of the hexane (8.3\% over 100K), and this
+limits the effect to ~20\% drop in thermal conductivity  (Table
+\ref{AuThiolToluene}).
+
+Although we have not mapped out the behavior at a large number of
+temperatures, is clear that there will be a strong temperature
+dependence in the interfacial conductance when the physical properties
+of one side of the interface (notably the density) change rapidly as a
+function of temperature.
+
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
@@ -643,329 +803,203 @@ Besides lower interfacial thermal conductance, surface
   \end{minipage}
 \end{table*}
 
-Besides lower interfacial thermal conductance, surfaces in relatively
-high temperatures are susceptible to reconstructions, when
-butanethiols have a full coverage on the Au(111) surface. These
-reconstructions include surface Au atoms migrated outward to the S
-atom layer, and butanethiol molecules embedded into the original
-surface Au layer. The driving force for this behavior is the strong
-Au-S interactions in our simulations. And these reconstructions lead
-to higher ratio of Au-S attraction and thus is energetically
-favorable. Furthermore, this phenomenon agrees with experimental
-results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. Vlugt
-{\it et al.} had kept their Au(111) slab rigid so that their
-simulations can reach 300K without surface reconstructions. Without
-this practice, simulating 100\% thiol covered interfaces under higher
-temperatures could hardly avoid surface reconstructions. However, our
-measurement is based on assuming homogeneity on $x$ and $y$ dimensions
-so that measurement of $T$ at particular $z$ would be an effective
-average of the particles of the same type. Since surface
-reconstructions could eliminate the original $x$ and $y$ dimensional
-homogeneity, measurement of $G$ is more difficult to conduct under
-higher temperatures. Therefore, most of our measurements are
-undertaken at $\langle T\rangle\sim$200K.
+Besides the lower interfacial thermal conductance, surfaces at
+relatively high temperatures are susceptible to reconstructions,
+particularly when butanethiols fully cover the Au(111) surface. These
+reconstructions include surface Au atoms which migrate outward to the
+S atom layer, and butanethiol molecules which embed into the surface
+Au layer. The driving force for this behavior is the strong Au-S
+interactions which are modeled here with a deep Lennard-Jones
+potential. This phenomenon agrees with reconstructions that have beeen
+experimentally
+observed.\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}.  Vlugt
+{\it et al.} kept their Au(111) slab rigid so that their simulations
+could reach 300K without surface
+reconstructions.\cite{vlugt:cpc2007154} Since surface reconstructions
+blur the interface, the measurement of $G$ becomes more difficult to
+conduct at higher temperatures.  For this reason, most of our
+measurements are undertaken at $\langle T\rangle\sim$200K where
+reconstruction is minimized.
 
 However, when the surface is not completely covered by butanethiols,
-the simulated system is more resistent to the reconstruction
-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.
+the simulated system appears to be more resistent to the
+reconstruction. O ur Au / butanethiol / toluene system had the Au(111)
+surfaces 90\% covered by butanethiols, but did not see this above
+phenomena even at $\langle T\rangle\sim$300K.  That said, we did
+observe butanethiols migrating to neighboring three-fold sites during
+a simulation.  Since the interface persisted in these simulations,
+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. 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.
+\section{Discussion}
 
-\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}
+The primary result of this work is that the capping agent acts as an
+efficient thermal coupler between solid and solvent phases.  One of
+the ways the capping agent can carry out this role is to down-shift
+between the phonon vibrations in the solid (which carry the heat from
+the gold) and the molecular vibrations in the liquid (which carry some
+of the heat in the solvent).
 
-It turned out that with partial covered butanethiol on the Au(111)
-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.
+To investigate the mechanism of interfacial thermal conductance, the
+vibrational power spectrum was computed. Power spectra were taken for
+individual components in different simulations. To obtain these
+spectra, simulations were run after equilibration in the
+microcanonical (NVE) ensemble and without a thermal
+gradient. Snapshots of configurations were collected at a frequency
+that is higher than that of the fastest vibrations 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}
+The power spectrum is constructed via a Fourier transform of the
+symmetrized velocity autocorrelation function,
+\begin{equation}
+  \hat{f}(\omega) =
+  \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
+\label{fourier} 
+\end{equation}
 
-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.
+\subsection{The role of specific vibrations}
+The vibrational spectra for gold slabs in different environments are
+shown as in Figure \ref{specAu}. Regardless of the presence of
+solvent, the gold surfaces which are covered by butanethiol molecules
+exhibit an additional peak observed at a frequency of
+$\sim$165cm$^{-1}$.  We attribute this peak to the S-Au bonding
+vibration. This vibration enables efficient thermal coupling of the
+surface Au layer to the capping agents. Therefore, in our simulations,
+the Au / S interfaces do not appear to be the primary barrier to
+thermal transport when compared with the butanethiol / solvent
+interfaces.
 
-Interestingly, as one could observe from our results, the maximum
-conductance enhancement (largest $G$) happens while the surfaces are
-about 75\% covered with butanethiols. This again indicates that
-solvent-capping agent contact has an important role of the thermal
-transport process. Slightly lower butanethiol coverage allows small
-gaps between butanethiols to form. And these gaps could be filled with
-solvent molecules, which acts like ``heat conductors'' on the
-surface. The higher degree of interaction between these solvent
-molecules and capping agents increases the enhancement effect and thus
-produces a higher $G$ than densely packed butanethiol arrays. However,
-once this maximum conductance enhancement is reached, $G$ decreases
-when butanethiol coverage continues to decrease. Each capping agent
-molecule reaches its maximum capacity for thermal
-conductance. Therefore, even higher solvent-capping agent contact
-would not offset this effect. Eventually, when butanethiol coverage
-continues to decrease, solvent-capping agent contact actually
-decreases with the disappearing of butanethiol molecules. In this
-case, $G$ decrease could not be offset but instead accelerated. [MAY 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
-transport process. The deuterated hexane (UA) results do not appear to
-be much different from those of normal hexane (UA), given that
-butanethiol (UA) is non-deuterated for both solvents. These UA model
-studies, even though eliminating C-H vibration samplings, still have
-C-C vibrational frequencies different from each other. However, these
-differences in the infrared range do not seem to produce an observable
-difference for the results of $G$ (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}
+\includegraphics[width=\linewidth]{vibration}
+\caption{Vibrational power spectra for gold in different solvent
+  environments.  The presence of the butanethiol capping molecules
+  adds a vibrational peak at $\sim$165cm$^{-1}$. The butanethiol
+  spectra exhibit a corresponding peak.}
+\label{specAu}
 \end{figure}
 
-Furthermore, results for rigid body toluene solvent, as well as other
-UA-hexane solvents, are reasonable within the general experimental
-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.
+Also in this figure, we show the vibrational power spectrum for the
+bound butanethiol molecules, which also exhibits the same
+$\sim$165cm$^{-1}$ peak.
 
-However, results for Au-butanethiol/toluene do not show an identical
-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 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
-increasing until the capping agent coverage reaches a relatively low
-level. This becomes an offset for decreasing butanethiol molecules on
-its effect to the process of interfacial thermal transport. Thus, one
-can see a plateau of $G$ vs. butanethiol coverage in our results.
+\subsection{Overlap of power spectra}
+A comparison of the results obtained from the two different organic
+solvents can also provide useful information of the interfacial
+thermal transport process.  In particular, the vibrational overlap
+between the butanethiol and the organic solvents suggests a highly
+efficient thermal exchange between these components.  Very high
+thermal conductivity was observed when AA models were used and C-H
+vibrations were treated classically.  The presence of extra degrees of
+freedom in the AA force field yields higher heat exchange rates
+between the two phases and results in a much higher conductivity than
+in the UA force field.
 
-\subsection{Influence of Chosen Molecule Model on $G$}
-In addition to UA solvent/capping agent models, AA models are included
-in our simulations as well. Besides simulations of the same (UA or AA)
-model for solvent and capping agent, different models can be applied
-to different components. Furthermore, regardless of models chosen,
-either the solvent or the capping agent can be deuterated, similar to
-the previous section. Table \ref{modelTest} summarizes the results of
-these studies.
+The similarity in the vibrational modes available to solvent and
+capping agent can be reduced by deuterating one of the two components
+(Fig. \ref{aahxntln}).  Once either the hexanes or the butanethiols
+are deuterated, one can observe a significantly lower $G$ and
+$G^\prime$ values (Table \ref{modelTest}).
 
-\begin{table*}
-  \begin{minipage}{\linewidth}
-    \begin{center}
-      
-      \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. (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}{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    & 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}
-    \end{center}
-  \end{minipage}
-\end{table*}
-
-To facilitate direct comparison, the same system with differnt models
-for different components uses the same length scale for their
-simulation cells. Without the presence of capping agent, using
-different models for hexane yields similar results for both $G$ and
-$G^\prime$, and these two definitions agree with eath other very
-well. This indicates very weak interaction between the metal and the
-solvent, and is a typical case for acoustic impedance mismatch between
-these two phases.
-
-As for Au(111) surfaces completely covered by butanethiols, the choice
-of models for capping agent and solvent could impact the measurement
-of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
-interfaces, using AA model for both butanethiol and hexane yields
-substantially higher conductivity values than using UA model for at
-least one component of the solvent and capping agent, which exceeds
-the 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
+\caption{Spectra obtained for all-atom (AA) Au / butanethiol / 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).}
+  vibrational overlap with hexane decreases significantly.  Since
+  aromatic molecules and the butanethiol are vibrationally dissimilar,
+  the change is not as dramatic when toluene is the solvent (right).}
 \label{aahxntln}
 \end{figure}
 
-However, for Au-butanethiol/toluene interfaces, having the AA
+For the Au / butanethiol / toluene interfaces, having the AA
 butanethiol deuterated did not yield a significant change in the
-measurement results. Compared to the C-H vibrational overlap between
-hexane and butanethiol, both of which have alkyl chains, that overlap
-between toluene and butanethiol is not so significant and thus does
-not have as much contribution to the 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.
+measured conductance. Compared to the C-H vibrational overlap between
+hexane and butanethiol, both of which have alkyl chains, the overlap
+between toluene and butanethiol is not as significant and thus does
+not contribute as much to the heat exchange process. 
 
-Although the QSC model for Au is known to predict an overly low value
-for bulk metal gold conductivity\cite{kuang:164101}, our computational
-results for $G$ and $G^\prime$ do not seem to be affected by this
-drawback of the model for metal. Instead, our results suggest that the
-modeling of interfacial thermal transport behavior relies mainly on
-the accuracy of the interaction descriptions between components
-occupying the interfaces.
+Previous observations of Zhang {\it et al.}\cite{hase:2010} indicate
+that the {\it intra}molecular heat transport due to alkylthiols is
+highly efficient.  Combining our observations with those of Zhang {\it
+  et al.}, it appears that butanethiol acts as a channel to expedite
+heat flow from the gold surface and into the alkyl chain.  The
+acoustic impedance mismatch between the metal and the liquid phase can
+therefore be effectively reduced with the presence of suitable capping
+agents.
 
-\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.
+Deuterated models in the UA force field did not decouple the thermal
+transport as well as in the AA force field.  The UA models, even
+though they have eliminated the high frequency C-H vibrational
+overlap, still have significant overlap in the lower-frequency
+portions of the infrared spectra (Figure \ref{uahxnua}).  Deuterating
+the UA models did not decouple the low frequency region enough to
+produce an observable difference for the results of $G$ (Table
+\ref{modelTest}). 
 
-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.
-
-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.
-
 \begin{figure}
-\includegraphics[width=\linewidth]{vibration}
-\caption{Vibrational spectra obtained for gold in different
-  environments.}
-\label{specAu}
+\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}
 
-[MAY ADD COMPARISON OF AU SLAB WIDTHS BUT NOT MUCH TO TALK ABOUT...]
-
 \section{Conclusions}
-The NIVS algorithm we developed has been applied to simulations of
-Au-butanethiol surfaces with organic solvents. This algorithm allows
-effective unphysical thermal flux transferred between the metal and
-the liquid phase. With the flux applied, we were able to measure the
-corresponding thermal gradient and to obtain interfacial thermal
-conductivities. 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.
+The NIVS algorithm has been applied to simulations of
+butanethiol-capped Au(111) surfaces in the presence of organic
+solvents. This algorithm allows the application of unphysical thermal
+flux to transfer heat between the metal and the liquid phase. With the
+flux applied, we were able to measure the corresponding thermal
+gradients and to obtain interfacial thermal conductivities. Under
+steady states, 2-3 ns trajectory simulations are sufficient for
+computation of this quantity.
 
-Our simulations have seen significant conductance enhancement with the
-presence of capping agent, compared to the bare gold / liquid
+Our simulations have seen significant conductance enhancement in the
+presence of capping agent, compared with the bare gold / liquid
 interfaces. The acoustic impedance mismatch between the metal and the
-liquid phase is effectively eliminated by proper capping
+liquid phase is effectively eliminated by a chemically-bonded capping
 agent. Furthermore, the coverage precentage of the capping agent plays
 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.
+process. Moderately low coverages allow higher contact between capping
+agent and solvent, and thus could further enhance the heat transfer
+process, giving a non-monotonic behavior of conductance with
+increasing coverage.
 
-Our measurement results, particularly of the UA models, agree with
-available experimental data. This indicates that our force field
-parameters have a nice description of the interactions between the
-particles at the interfaces. AA models tend to overestimate the
+Our results, particularly using the UA models, agree well with
+available experimental data.  The AA models tend to overestimate the
 interfacial thermal conductance in that the classically treated C-H
-vibration would be overly sampled. Compared to the AA models, the UA
-models have higher computational efficiency with satisfactory
-accuracy, and thus are preferable in interfacial thermal transport
-modelings. Of the two definitions for $G$, the discrete form
+vibrations become too easily populated. Compared to the AA models, the
+UA models have higher computational efficiency with satisfactory
+accuracy, and thus are preferable in modeling interfacial thermal
+transport.
+
+Of the two definitions for $G$, the discrete form
 (Eq. \ref{discreteG}) was easier to use and gives out relatively
 consistent results, while the derivative form (Eq. \ref{derivativeG})
 is not as versatile. Although $G^\prime$ gives out comparable results
 and follows similar trend with $G$ when measuring close to fully
-covered or bare surfaces, the spatial resolution of $T$ profile is
-limited for accurate computation of derivatives data.
+covered or bare surfaces, the spatial resolution of $T$ profile
+required for the use of a derivative form is limited by the number of
+bins and the sampling required to obtain thermal gradient information.
 
-Vlugt {\it et al.} has investigated the surface thiol structures for
-nanocrystal gold and pointed out that they differs from those of the
-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.
+Vlugt {\it et al.} have investigated the surface thiol structures for
+nanocrystalline gold and pointed out that they differ from those of
+the Au(111) surface.\cite{landman:1998,vlugt:cpc2007154} This
+difference could also cause differences in the interfacial thermal
+transport behavior. To investigate this problem, one would need an
+effective method for applying thermal gradients in non-planar
+(i.e. spherical) geometries. 
 
 \section{Acknowledgments}
 Support for this project was provided by the National Science
 Foundation under grant CHE-0848243. Computational time was provided by
 the Center for Research Computing (CRC) at the University of Notre
-Dame. \newpage
+Dame. 
+\newpage
 
 \bibliography{interfacial}
 
