group/interfacial/interfacial.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{multirow}
%\usepackage{booktabs}
%\usepackage{bibentry}
%\usepackage{mathrsfs}
%\usepackage[ref]{overcite}
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

%\renewcommand\citemid{\ } % no comma in optional reference note
\bibpunct{[}{]}{,}{s}{}{;}
\bibliographystyle{aip}

\begin{document}

\title{Simulating interfacial thermal conductance at metal-solvent
  interfaces: the role of chemical capping agents}

\author{Shenyu Kuang and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}

With the Non-Isotropic Velocity Scaling algorithm (NIVS) we have
developed, an unphysical thermal flux can be effectively set up even
for non-homogeneous systems like interfaces in non-equilibrium
molecular dynamics simulations. In this work, this algorithm is
applied for simulating thermal conductance at metal / organic solvent
interfaces with various coverages of butanethiol capping
agents. Different solvents and force field models were tested. Our
results suggest that the United-Atom models are able to provide an
estimate of the interfacial thermal conductivity comparable to
experiments in our simulations with satisfactory computational
efficiency. From our results, the acoustic impedance mismatch between
metal and liquid phase is effectively reduced by the capping
agents, and thus leads to interfacial thermal conductance
enhancement. Furthermore, this effect is closely related to the
capping agent coverage on the metal surfaces and the type of solvent
molecules, and is affected by the models used in the simulations.

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
[BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM]
Interfacial thermal conductance is extensively studied both
experimentally and computationally, and systems with interfaces
present are generally heterogeneous. Although interfaces are commonly
barriers to heat transfer, it has been
reported\cite{doi:10.1021/la904855s} that under specific circustances,
e.g. with certain capping agents present on the surface, interfacial
conductance can be significantly enhanced. However, heat conductance
of molecular and nano-scale interfaces will be affected by the
chemical details of the surface and is challenging to
experimentalist. The lower thermal flux through interfaces is even
more difficult to measure with EMD and forward NEMD simulation
methods. Therefore, developing good simulation methods will be
desirable in order to investigate thermal transport across interfaces.

Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
retains the desirable features of RNEMD (conservation of linear
momentum and total energy, compatibility with periodic boundary
conditions) while establishing true thermal distributions in each of
the two slabs. Furthermore, it allows more effective thermal exchange
between particles of different identities, and thus enables extensive
study of interfacial conductance.

\section{Methodology}
\subsection{Algorithm}
[BACKGROUND FOR MD METHODS]
There have been many algorithms for computing thermal conductivity
using molecular dynamics simulations. However, interfacial conductance
is at least an order of magnitude smaller. This would make the
calculation even more difficult for those slowly-converging
equilibrium methods. Imposed-flux non-equilibrium
methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
the response of temperature or momentum gradients are easier to
measure than the flux, if unknown, and thus, is a preferable way to
the forward NEMD methods. Although the momentum swapping approach for
flux-imposing can be used for exchanging energy between particles of
different identity, the kinetic energy transfer efficiency is affected
by the mass difference between the particles, which limits its
application on heterogeneous interfacial systems.

The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach in
non-equilibrium MD simulations is able to impose relatively large
kinetic energy flux without obvious perturbation to the velocity
distribution of the simulated systems. Furthermore, this approach has
the advantage in heterogeneous interfaces in that kinetic energy flux
can be applied between regions of particles of arbitary identity, and
the flux quantity is not restricted by particle mass difference.

The 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. 

\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:
\begin{equation}
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.

When the interfacial conductance is {\it not} small, two ways can be
used to define $G$.

One way is to assume the temperature is discretely different on two
sides of the interface, $G$ can be calculated with the thermal flux
applied $J$ and the maximum temperature difference measured along the
thermal gradient max($\Delta T$), which occurs at the interface, as:
\begin{equation}
G=\frac{J}{\Delta T}
\label{discreteG}
\end{equation}

The other approach is to assume a continuous temperature profile along
the thermal gradient axis (e.g. $z$) and define $G$ at the point where
the magnitude of thermal conductivity $\lambda$ change reach its
maximum, given that $\lambda$ is well-defined throughout the space:
\begin{equation}
G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
         = \Big|\frac{\partial}{\partial z}\left(-J_z\Big/
           \left(\frac{\partial T}{\partial z}\right)\right)\Big|
         = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
         \Big/\left(\frac{\partial T}{\partial z}\right)^2
\label{derivativeG}
\end{equation}

With the temperature profile obtained from simulations, one is able to
approximate the first and second derivatives of $T$ with finite
difference method and thus calculate $G^\prime$.

In what follows, both definitions are used for calculation and comparison.

[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}.

\begin{figure}
\includegraphics[width=\linewidth]{demoPic}
\caption{A sample showing how a metal slab has its (111) surface
  covered by capping agent molecules and solvated by hexane.}
\label{demoPic}
\end{figure}

With a simulation cell setup following the above manner, one is able
to equilibrate the system and impose an unphysical thermal flux
between the liquid and the metal phase with the NIVS algorithm. Under
a stablized thermal gradient induced by periodically applying the
unphysical flux, one is able to obtain a temperature profile and the
physical thermal flux corresponding to it, which equals to the
unphysical flux applied by NIVS. These data enables the evaluation of
the interfacial thermal conductance of a surface. Figure \ref{gradT}
is an example how those stablized thermal gradient can be used to
obtain the 1st and 2nd derivatives of the temperature profile.

\begin{figure}
\includegraphics[width=\linewidth]{gradT}
\caption{The 1st and 2nd derivatives of temperature profile can be
  obtained with finite difference approximation.}
\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]

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 spacing 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
solvent region and lead to undesired phenomena such as solvent boiling
or freezing when a thermal flux is applied. Conversely, too few
solvent molecules would change the normal behavior of the liquid
phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
these extreme cases did not happen to our simulations. And the
corresponding spacing is usually $35 \sim 60$\AA.

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.

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:
\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
         = |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big|
         \Big/\left(\frac{\partial T}{\partial n}\right)^2
\label{derivativeG2}
\end{equation}

\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.

The Au-Au interactions in metal lattice slab is described by the
quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
potentials include zero-point quantum corrections and are
reparametrized for accurate surface energies compared to the
Sutton-Chen potentials\cite{Chen90}. 

Figure [REF] 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]

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]

The capping agent in our simulations, the butanethiol molecules can
either use UA or AA model. The TraPPE-UA force fields 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 [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: 


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
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]

However, the Lennard-Jones parameters between Au and other types of
particles in our simulations are not yet well-established. For these
interactions, we attempt to derive their parameters using the Mixing
Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters
for Au is first extracted from the Au-CH$_x$ parameters by applying
the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
parameters in our simulations.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      \caption{Lennard-Jones parameters for Au-non-Metal
        interactions in our simulations.}
      
      \begin{tabular}{ccc}
        \hline\hline
        Non-metal atom   & $\sigma$ & $\epsilon$ \\
        (or pseudo-atom) & \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  \\
        \hline\hline
      \end{tabular}
      \label{MnM}
    \end{center}
  \end{minipage}
\end{table*}


\section{Results and Discussions}
[MAY HAVE A BRIEF SUMMARY]
\subsection{How Simulation Parameters Affects $G$}
[MAY NOT PUT AT FIRST]
We have varied our protocol or other parameters of the simulations in
order to investigate how these factors would affect the measurement of
$G$'s. It turned out that while some of these parameters would not
affect the results substantially, some other changes to the
simulations would have a significant impact on the measurement
results.

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. 

As stated in our computational details, the spacing filled with
solvent molecules can be chosen within a range. This allows some
change of solvent molecule numbers for the same Au-butanethiol
surfaces. We did this study on our Au-butanethiol/hexane
simulations. Nevertheless, the results obtained from systems of
different $N_{hexane}$ did not indicate that the measurement of $G$ is
susceptible to this parameter. For computational efficiency concern,
smaller system size would be preferable, given that the liquid phase
structure is not affected.

Our NIVS algorithm allows change of unphysical thermal flux both in
direction and in quantity. This feature extends our investigation of
interfacial thermal conductance. However, the magnitude of this
thermal flux is not arbitary if one aims to obtain a stable and
reliable thermal gradient. A temperature profile would be
substantially affected by noise when $|J_z|$ has a much too low
magnitude; while an excessively large $|J_z|$ that overwhelms the
conductance capacity of the interface would prevent a thermal gradient
to reach a stablized steady state. NIVS has the advantage of allowing
$J$ to vary in a wide range such that the optimal flux range for $G$
measurement can generally be simulated by the algorithm. Within the
optimal range, we were able to study how $G$ would change according to
the thermal flux across the interface. For our simulations, we denote
$J_z$ to be positive when the physical thermal flux is from the liquid
to metal, and negative vice versa. The $G$'s measured under different
$J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These
results do not suggest that $G$ is dependent on $J_z$ within this flux
range. The linear response of flux to thermal gradient simplifies our
investigations in that we can rely on $G$ measurement with only a
couple $J_z$'s and do not need to test a large series of fluxes.

%ADD MORE TO TABLE
\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      \caption{Computed interfacial thermal conductivity ($G$ and
        $G^\prime$) values for the 100\% covered Au-butanethiol/hexane
        interfaces with UA model and different hexane molecule numbers
        at different temperatures using a range of energy fluxes.}
      
      \begin{tabular}{cccccccc}
        \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$) &
        \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() \\
        \hline\hline
      \end{tabular}
      \label{AuThiolHexaneUA}
    \end{center}
  \end{minipage}
\end{table*}

Furthermore, we also attempted to increase system average temperatures
to above 200K. These simulations are first equilibrated in the NPT
ensemble under normal pressure. As stated above, the TraPPE-UA model
for hexane tends to predict a lower boiling point. In our simulations,
hexane had diffculty to remain in liquid phase when NPT equilibration
temperature is higher than 250K. Additionally, the equilibrated liquid
hexane density under 250K becomes lower than experimental value. This
expanded liquid phase leads to lower contact between hexane and
butanethiol as well.[MAY NEED FIGURE] And this reduced contact would
probably be accountable for a lower interfacial thermal conductance,
as shown in Table \ref{AuThiolHexaneUA}. 

A similar study for TraPPE-UA toluene agrees with the above result as
well. Having a higher boiling point, toluene tends to remain liquid in
our simulations even equilibrated under 300K in NPT
ensembles. Furthermore, the expansion of the toluene liquid phase is
not as significant as that of the hexane. This prevents severe
decrease of liquid-capping agent contact and the results (Table
\ref{AuThiolToluene}) show only a slightly decreased interface
conductance. Therefore, solvent-capping agent contact should play an
important role in the thermal transport process across the interface
in that higher degree of contact could yield increased conductance.

[ADD 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.}
      
      \begin{tabular}{cccc}
        \hline\hline
        $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
        (K) & (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() \\
        \hline\hline
      \end{tabular}
      \label{AuThiolToluene}
    \end{center}
  \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.

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 $<T>\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.

\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.

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. 

From Table \ref{tlnUhxnUhxnD}, one can see the significance of the
presence of capping agents. Even when a fraction of the Au(111)
surface sites are covered with butanethiols, the conductivity would
see an enhancement by at least a factor of 3. This indicates the
important role cappping agent is playing for thermal transport
phenomena on metal/organic solvent surfaces.

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.

A comparison of the results obtained from differenet organic solvents
can also provide useful information of the interfacial thermal
transport process. The deuterated hexane (UA) results do not appear to
be much different from those of normal hexane (UA), given that
butanethiol (UA) is non-deuterated for both solvents. These UA model
studies, even though eliminating C-H vibration samplings, still have
C-C vibrational frequencies different from each other. However, these
differences in the infrared range do not seem to produce an observable
difference for the results of $G$. [MAY NEED FIGURE]

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.

However, results for Au-butanethiol/toluene do not show an identical
trend with those for Au-butanethiol/hexane in that $G$'s remain at
approximately the same magnitue when butanethiol coverage differs from
25\% to 75\%. This might be rooted in the molecule shape difference
for plane-like toluene and chain-like {\it n}-hexane. Due to this
difference, toluene molecules have more difficulty in occupying
relatively small gaps among capping agents when their coverage is not
too low. Therefore, the solvent-capping agent contact may keep
increasing until the capping agent coverage reaches a relatively low
level. This becomes an offset for decreasing butanethiol molecules on
its effect to the process of interfacial thermal transport. Thus, one
can see a plateau of $G$ vs. butanethiol coverage in our results.

[NEED ERROR ESTIMATE, 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
to different components. Furthermore, regardless of models chosen,
either the solvent or the capping agent can be deuterated, similar to
the previous section. Table \ref{modelTest} summarizes the results of
these studies.

[MORE DATA; ERROR ESTIMATE]
\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed interfacial thermal conductivity ($G$ and
        $G^\prime$) values for interfaces using various models for
        solvent and capping agent (or without capping agent) at
        $\langle T\rangle\sim$200K.}
      
      \begin{tabular}{ccccc}
        \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() \\
        \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 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.

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.


Although the QSC model for Au is known to predict an overly low value
for bulk metal gold conductivity[CITE NIVSRNEMD], 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.

\subsection{Mechanism of Interfacial Thermal Conductance Enhancement
  by Capping Agent}
%OR\subsection{Vibrational spectrum study on conductance mechanism}

[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL, EQN'S]

To investigate the mechanism of this interfacial thermal conductance,
the vibrational spectra of various gold systems were obtained and are
shown as in the upper panel of Fig. \ref{vibration}. To obtain these
spectra, one first runs a simulation in the NVE ensemble and collects
snapshots of configurations; these configurations are used to compute
the velocity auto-correlation functions, which is used to construct a
power spectrum via a Fourier transform.

 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}
\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.

\section{Conclusions}
The NIVS algorithm we developed has been applied to simulations of
Au-butanethiol surfaces with organic solvents. This algorithm allows
effective unphysical thermal flux transferred between the metal and
the liquid phase. With the flux applied, we were able to measure the
corresponding thermal gradient and to obtain interfacial thermal
conductivities. Our simulations have seen significant conductance
enhancement with the presence of capping agent, compared to the bare
gold/liquid interfaces. The acoustic impedance mismatch between the
metal and the liquid phase is effectively eliminated by proper capping
agent. Furthermore, the coverage precentage of the capping agent plays
an important role in the interfacial thermal transport process.

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
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.

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.


\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

\bibliography{interfacial}

\end{doublespace}
\end{document}

Revision:	3732
Committed:	Fri Jul 8 17:07:00 2011 UTC (13 years ago) by skuang
Content type:	application/x-tex
File size:	40320 byte(s)
Log Message:	more results and discussions, added conclusions and abstract