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{[}{]}{,}{n}{}{;}
\bibliographystyle{achemso}

\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}
Interfacial thermal conductance is extensively studied both
experimentally and computationally\cite{cahill:793}, due to its
importance in nanoscale science and technology. Reliability of
nanoscale devices depends on their thermal transport
properties. Unlike bulk homogeneous materials, nanoscale materials
features significant presence of interfaces, and these interfaces
could dominate the heat transfer behavior of these
materials. Furthermore, these materials are generally heterogeneous,
which challenges traditional research methods for homogeneous
systems.

Heat conductance of molecular and nano-scale interfaces will be
affected by the chemical details of the surface. Experimentally,
various interfaces have been investigated for their thermal
conductance properties. Wang {\it et al.} studied heat transport
through long-chain hydrocarbon monolayers on gold substrate at
individual molecular level\cite{Wang10082007}; Schmidt {\it et al.}
studied the role of CTAB on thermal transport between gold nanorods
and solvent\cite{doi:10.1021/jp8051888}; Juv\'e {\it et al.} studied
the cooling dynamics, which is controlled by thermal interface
resistence of glass-embedded metal
nanoparticles\cite{PhysRevB.80.195406}. Although interfaces are
commonly barriers for heat transport, Alper {\it et al.} suggested
that specific ligands (capping agents) could completely eliminate this
barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}.

Theoretical and computational models have also been used to study the
interfacial thermal transport in order to gain an understanding of
this phenomena at the molecular level. Recently, Hase and coworkers
employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
study thermal transport from hot Au(111) substrate to a self-assembled
monolayer of alkylthiol with relatively long chain (8-20 carbon
atoms)\cite{hase:2010,hase:2011}. However, ensemble averaged
measurements for heat conductance of interfaces between the capping
monolayer on Au and a solvent phase has yet to be studied. 
The comparatively low thermal flux through interfaces is
difficult to measure with Equilibrium MD or forward NEMD simulation
methods. Therefore, the Reverse NEMD (RNEMD) methods would have the
advantage of having this difficult to measure flux known when studying
the thermal transport across interfaces, given that the simulation
methods being able to effectively apply an unphysical flux in
non-homogeneous systems.

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 effective thermal exchange
between particles of different identities, and thus makes the study of
interfacial conductance much simpler.

The work presented here deals with the Au(111) surface covered to
varying degrees by butanethiol, a capping agent with short carbon
chain, and solvated with organic solvents of different molecular
properties. Different models were used for both the capping agent and
the solvent force field parameters. Using the NIVS algorithm, the
thermal transport across these interfaces was studied and the
underlying mechanism for the phenomena was investigated.

[MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]

\section{Methodology}
\subsection{Imposd-Flux Methods in MD Simulations}
Steady state MD simulations has the advantage that not many
trajectories are needed to study the relationship between thermal flux
and thermal gradients. For systems including low conductance
interfaces one must have a method capable of generating or measuring
relatively small fluxes, compared to those required for bulk
conductivity. This requirement makes the calculation even more
difficult for those slowly-converging equilibrium
methods\cite{Viscardy:2007lq}. Forward methods may impose gradient,
but in interfacial conditions it is not clear what behavior to impose
at the interfacial boundaries. Imposed-flux reverse non-equilibrium
methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
the thermal response becomes easier to measure than the flux. Although
M\"{u}ller-Plathe's original momentum swapping approach can be used
for exchanging energy between particles of different identity, the
kinetic energy transfer efficiency is affected by the mass difference
between the particles, which limits its application on heterogeneous
interfacial systems.

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

The NIVS algorithm scales the velocity vectors in two separate regions
of a simulation system with respective diagonal scaling matricies. To
determine these scaling factors in the matricies, a set of equations
including linear momentum conservation and kinetic energy conservation
constraints and target energy flux satisfaction is solved. With the
scaling operation applied to the system in a set frequency, bulk
temperature gradients can be easily established, and these can be used
for computing thermal conductivities. The NIVS algorithm conserves
momenta and energy and does not depend on an external thermostat.

\subsection{Defining Interfacial Thermal Conductivity $G$}
Given a system with thermal gradients and the corresponding thermal
flux, for interfaces with a relatively low interfacial conductance,
the bulk regions on either side of an interface rapidly come to a
state in which the two phases have relatively homogeneous (but
distinct) temperatures. The interfacial thermal conductivity $G$ can
therefore be approximated as:
\begin{equation}
G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
    \langle T_\mathrm{cold}\rangle \right)}
\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, there are two
ways to define $G$.

One way is to assume the temperature is discrete on the two sides of
the interface. $G$ can be calculated using the applied thermal flux
$J$ and the maximum temperature difference measured along the thermal
gradient max($\Delta T$), which occurs at the Gibbs deviding surface
(Figure \ref{demoPic}):
\begin{equation}
G=\frac{J}{\Delta T}
\label{discreteG}
\end{equation}

\begin{figure}
\includegraphics[width=\linewidth]{method}
\caption{Interfacial conductance can be calculated by applying an
  (unphysical) kinetic energy flux between two slabs, one located
  within the metal and another on the edge of the periodic box.  The
  system responds by forming a thermal response or a gradient.  In
  bulk liquids, this gradient typically has a single slope, but in
  interfacial systems, there are distinct thermal conductivity
  domains.  The interfacial conductance, $G$ is found by measuring the
  temperature gap at the Gibbs dividing surface, or by using second
  derivatives of the thermal profile.}
\label{demoPic}
\end{figure}

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

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

In what follows, both definitions have been used for calculation and
are compared in the results.

To compare the above definitions ($G$ and $G^\prime$), we have modeled
a metal slab with its (111) surfaces perpendicular to the $z$-axis of
our simulation cells. Both with and without capping agents on the
surfaces, the metal slab is solvated with simple organic solvents, as
illustrated in Figure \ref{gradT}.

With the simulation cell described above, we are able to equilibrate
the system and impose an unphysical thermal flux between the liquid
and the metal phase using the NIVS algorithm. By periodically applying
the unphysical flux, we are able to obtain a temperature profile and
its spatial derivatives. These quantities enable the evaluation of the
interfacial thermal conductance of a surface. Figure \ref{gradT} is an
example of how an applied thermal flux can be used to obtain the 1st
and 2nd derivatives of the temperature profile.

\begin{figure}
\includegraphics[width=\linewidth]{gradT}
\caption{A sample of Au-butanethiol/hexane interfacial system and the
  temperature profile after a kinetic energy flux is imposed to
  it. The 1st and 2nd derivatives of the temperature profile can be
  obtained with finite difference approximation (lower panel).}
\label{gradT}
\end{figure}

\section{Computational Details}
\subsection{Simulation Protocol}
The NIVS algorithm has been implemented in our MD simulation code,
OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
simulations. Different metal slab thickness (layer numbers of Au) was
simulated. Metal slabs were first equilibrated under atmospheric
pressure (1 atm) and a desired temperature (e.g. 200K). After
equilibration, butanethiol capping agents were placed at three-fold
hollow sites on the Au(111) surfaces. These sites could be either a
{\it fcc} or {\it hcp} site. However, Hase {\it et al.} found that
they are equivalent in a heat transfer process\cite{hase:2010}, so
they are not distinguished in our study. The maximum butanethiol
capacity on Au surface is $1/3$ of the total number of surface Au
atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
series of different coverages was derived by evenly eliminating
butanethiols on the surfaces, and was investigated in order to study
the relation between coverage and interfacial conductance.

The capping agent molecules were allowed to migrate during the
simulations. They distributed themselves uniformly and sampled a
number of three-fold sites throughout out study. Therefore, the
initial configuration would not noticeably affect the sampling of a
variety of configurations of the same coverage, and the final
conductance measurement would be an average effect of these
configurations explored in the simulations. [MAY NEED SNAPSHOTS]

After the modified Au-butanethiol surface systems were equilibrated
under canonical ensemble, organic solvent molecules were packed in the
previously empty part of the simulation cells\cite{packmol}. Two
solvents were investigated, one which has little vibrational overlap
with the alkanethiol and a planar shape (toluene), and one which has
similar vibrational frequencies and chain-like shape ({\it n}-hexane).

The space filled by solvent molecules, i.e. the gap between
periodically repeated Au-butanethiol surfaces should be carefully
chosen. A very long length scale for the thermal gradient axis ($z$)
may cause excessively hot or cold temperatures in the middle of the
solvent region and lead to undesired phenomena such as solvent boiling
or freezing when a thermal flux is applied. Conversely, too few
solvent molecules would change the normal behavior of the liquid
phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
these extreme cases did not happen to our simulations. And the
corresponding spacing is usually $35 \sim 75$\AA.

The initial configurations generated are further equilibrated with the
$x$ and $y$ dimensions fixed, only allowing length scale change in $z$
dimension. This is to ensure that the equilibration of liquid phase
does not affect the metal crystal structure in $x$ and $y$ dimensions.
To investigate this effect, comparisons were made with simulations
that allow changes of $L_x$ and $L_y$ during NPT equilibration, and
the results are shown in later sections. After ensuring the liquid
phase reaches equilibrium at atmospheric pressure (1 atm), further
equilibration are followed under NVT and then NVE ensembles.

After the systems reach equilibrium, NIVS is implemented to impose a
periodic unphysical thermal flux between the metal and the liquid
phase. Most of our simulations are under an average temperature of
$\sim$200K. Therefore, this flux usually comes from the metal to the
liquid so that the liquid has a higher temperature and would not
freeze due to excessively low temperature. After this induced
temperature gradient is stablized, the temperature profile of the
simulation cell is recorded. To do this, the simulation cell is
devided evenly into $N$ slabs along the $z$-axis and $N$ is maximized
for highest possible spatial resolution but not too many to have some
slabs empty most of the time. The average temperatures of each slab
are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
the same, the derivatives of $T$ with respect to slab number $n$ can
be directly used for $G^\prime$ calculations:
\begin{equation}
G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
         \Big/\left(\frac{\partial T}{\partial z}\right)^2
         = |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}

All of the above simulation procedures use a time step of 1 fs. And
each equilibration / stabilization step usually takes 100 ps, or
longer, if necessary.

\subsection{Force Field Parameters}
Our simulations include various components. Figure \ref{demoMol}
demonstrates the sites defined for both United-Atom and All-Atom
models of the organic solvent and capping agent molecules in our
simulations. Force field parameter descriptions are needed for
interactions both between the same type of particles and between
particles of different species.

\begin{figure}
\includegraphics[width=\linewidth]{structures}
\caption{Structures of the capping agent and solvents utilized in
  these simulations. The chemically-distinct sites (a-e) are expanded
  in terms of constituent atoms for both United Atom (UA) and All Atom
  (AA) force fields.  Most parameters are from
  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} (UA) and
  \protect\cite{OPLSAA} (AA).  Cross-interactions with the Au atoms are given
  in Table \ref{MnM}.}
\label{demoMol}
\end{figure}

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

For both solvent molecules, straight chain {\it n}-hexane and aromatic
toluene, United-Atom (UA) and All-Atom (AA) models are used
respectively. The TraPPE-UA
parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
for our UA solvent molecules. In these models, sites are located at
the carbon centers for alkyl groups. Bonding interactions, including
bond stretches and bends and torsions, were used for intra-molecular
sites not separated by more than 3 bonds. Otherwise, for non-bonded
interactions, Lennard-Jones potentials are used. [CHECK CITATION]

By eliminating explicit hydrogen atoms, these models are simple and
computationally efficient, while maintains good accuracy. However, the
TraPPE-UA for alkanes is known to predict a lower boiling point than
experimental values. Considering that after an unphysical thermal flux
is applied to a system, the temperature of ``hot'' area in the liquid
phase would be significantly higher than the average of the system, to
prevent over heating and boiling of the liquid phase, the average
temperature in our simulations should be much lower than the liquid
boiling point.

For UA-toluene model, the non-bonded potentials between
inter-molecular sites have a similar Lennard-Jones formulation. For
intra-molecular interactions, considering the stiffness of the benzene
ring, rigid body constraints are applied for further computational
efficiency. All bonds in the benzene ring and between the ring and the
methyl group remain rigid during the progress of simulations.

Besides the TraPPE-UA models, AA models for both organic solvents are
included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
force field is used. Additional explicit hydrogen sites were
included. Besides bonding and non-bonded site-site interactions,
partial charges and the electrostatic interactions were added to each
CT and HC site. For toluene, the United Force Field developed by
Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} is
adopted. Without the rigid body constraints, bonding interactions were
included. For the aromatic ring, improper torsions (inversions) were
added as an extra potential for maintaining the planar shape.
[CHECK CITATION]

The capping agent in our simulations, the butanethiol molecules can
either use UA or AA model. The TraPPE-UA force fields includes
parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
UA butanethiol model in our simulations. The OPLS-AA also provides
parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
surfaces do not have the hydrogen atom bonded to sulfur. To adapt this
change and derive suitable parameters for butanethiol adsorbed on
Au(111) surfaces, we adopt the S parameters from Luedtke and
Landman\cite{landman:1998}[CHECK CITATION]
 and modify parameters for its neighbor C
atom for charge balance in the molecule. Note that the model choice
(UA or AA) of capping agent can be different from the
solvent. Regardless of model choice, the force field parameters for
interactions between capping agent and solvent can be derived using
Lorentz-Berthelot Mixing Rule:
\begin{eqnarray}
\sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
\epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
\end{eqnarray}

To describe the interactions between metal Au and non-metal capping
agent and solvent particles, we refer to an adsorption study of alkyl
thiols on gold surfaces by Vlugt {\it et
  al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones
form of potential parameters for the interaction between Au and
pseudo-atoms CH$_x$ and S based on a well-established and widely-used
effective potential of Hautman and Klein\cite{hautman:4994} for the
Au(111) surface. As our simulations require the gold lattice slab to
be non-rigid so that it could accommodate kinetic energy for thermal
transport study purpose, the pair-wise form of potentials is
preferred.

Besides, the potentials developed from {\it ab initio} calculations by
Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
interactions between Au and aromatic C/H atoms in toluene. A set of
pseudo Lennard-Jones parameters were provided for Au in their force
fields. By using the Mixing Rule, this can be used to derive pair-wise
potentials for non-bonded interactions between Au and non-metal sites.

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

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      \caption{Non-bonded interaction parameters (including cross
        interactions with Au atoms) for both force fields used in this
        work.}      
      \begin{tabular}{lllllll}
        \hline\hline
        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
        \hline
        United Atom (UA)
        &CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
        &CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
        &CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
        &CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
        \hline
        All Atom (AA) 
        &CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
        &CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
        &CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
        &HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
        &CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
        &HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
        \hline
        Both UA and AA 
        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
        \hline\hline
      \end{tabular}
      \label{MnM}
    \end{center}
  \end{minipage}
\end{table*}

\subsection{Vibrational Spectrum}
To investigate the mechanism of interfacial thermal conductance, the
vibrational spectrum is utilized as a complementary tool. Vibrational
spectra were taken for individual components in different
simulations. To obtain these spectra, simulations were run after
equilibration, in the NVE ensemble. Snapshots of configurations were
collected at a frequency that is higher than that of the fastest
vibrations occuring in the simulations. With these configurations, the
velocity auto-correlation functions can be computed:
\begin{equation}
C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
\label{vCorr}
\end{equation}
Followed by Fourier transforms, the power spectrum can be constructed:
\begin{equation}
\hat{f}(\omega) = \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
\label{fourier}
\end{equation}

\section{Results and Discussions}
In what follows, how the parameters and protocol of simulations would
affect the measurement of $G$'s is first discussed. With a reliable
protocol and set of parameters, the influence of capping agent
coverage on thermal conductance is investigated. Besides, different
force field models for both solvents and selected deuterated models
were tested and compared. Finally, a summary of the role of capping
agent in the interfacial thermal transport process is given.

\subsection{How Simulation Parameters Affects $G$}
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
crystalline Au structure, $L_x$ and $L_y$ would not change noticeably
after equilibration. Although $L_z$ could fluctuate $\sim$1\% after a
system is fully equilibrated in the NPT ensemble, this fluctuation, as
well as those of $L_x$ and $L_y$ (which is significantly smaller),
would not be magnified on the calculated $G$'s, as shown in Table
\ref{AuThiolHexaneUA}. This insensivity to $L_i$ fluctuations allows
reliable measurement of $G$'s without the necessity of extremely
cautious equilibration process.

As stated in our computational details, the spacing filled with
solvent molecules can be chosen within a range. This allows some
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{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.

\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. Error estimates indicated in parenthesis.}
      
      \begin{tabular}{ccccccc}
        \hline\hline
        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
        $J_z$ & $G$ & $G^\prime$ \\
        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
            & 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
            &     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
            &     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
            &     &     &       &  1.91 & 139(10)   & 101(10)   \\
            &     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
            & 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
            &     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
            &     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
            &     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
        \hline
        250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
            &     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
            & 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
            &     & No  & 0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
            &     &     &       &  1.44 & 76.2(5.0) & 64.8(3.8) \\
            &     &     &       & -0.95 & 56.4(2.5) & 54.4(1.1) \\
            &     &     &       & -1.85 & 47.8(1.1) & 53.5(1.5) \\
        \hline\hline
      \end{tabular}
      \label{AuThiolHexaneUA}
    \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 SLAB DENSITY 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.

\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. Error estimates indicated in parenthesis.}
      
      \begin{tabular}{ccccc}
        \hline\hline
        $\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
        (K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        200 & 0.933 &  2.15 & 204(12) & 113(12) \\
            &       & -1.86 & 180(3)  & 135(21) \\
            &       & -3.93 & 176(5)  & 113(12) \\
        \hline
        300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
            &       & -4.19 & 135(9)  & 113(12) \\
        \hline\hline
      \end{tabular}
      \label{AuThiolToluene}
    \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 had the Au(111) surfaces 90\%
covered by butanethiols, but did not see this above phenomena even at
$\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
capping agents could help prevent surface reconstruction in that they
provide other means of capping agent relaxation. It is observed that
butanethiols can migrate to their neighbor empty sites during a
simulation. Therefore, we were able to obtain $G$'s for these
interfaces even at a relatively high temperature without being
affected by surface reconstructions.

\subsection{Influence of Capping Agent Coverage on $G$}
To investigate the influence of butanethiol coverage on interfacial
thermal conductance, a series of different coverage Au-butanethiol
surfaces is prepared and solvated with various organic
molecules. These systems are then equilibrated and their interfacial
thermal conductivity are measured with our NIVS algorithm. Figure
\ref{coverage} demonstrates the trend of conductance change with
respect to different coverages of butanethiol. To study the isotope
effect in interfacial thermal conductance, deuterated UA-hexane is
included as well.

\begin{figure}
\includegraphics[width=\linewidth]{coverage}
\caption{Comparison of interfacial thermal conductivity ($G$) values
  for the Au-butanethiol/solvent interface with various UA models and
  different capping agent coverages at $\langle T\rangle\sim$200K
  using certain energy flux respectively.}
\label{coverage}
\end{figure}

It turned out that with partial covered butanethiol on the Au(111)
surface, the derivative definition for $G^\prime$
(Eq. \ref{derivativeG}) was difficult to apply, due to the difficulty
in locating the maximum of change of $\lambda$. Instead, the discrete
definition (Eq. \ref{discreteG}) is easier to apply, as the Gibbs
deviding surface can still be well-defined. Therefore, $G$ (not
$G^\prime$) was used for this section.

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

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

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

\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
  systems. When butanethiol is deuterated (lower left), its
  vibrational overlap with hexane would decrease significantly,
  compared with normal butanethiol (upper left). However, this
  dramatic change does not apply to toluene as much (right).}
\label{aahxntln}
\end{figure}

However, for Au-butanethiol/toluene interfaces, having the AA
butanethiol deuterated did not yield a significant change in the
measurement results. Compared to the C-H vibrational overlap between
hexane and butanethiol, both of which have alkyl chains, that overlap
between toluene and butanethiol is not so significant and thus does
not have as much contribution to the heat exchange
process. Conversely, extra degrees of freedom such as the C-H
vibrations could yield higher heat exchange rate between these two
phases and result in a much higher conductivity.

Although the QSC model for Au is known to predict an overly low value
for bulk metal gold conductivity\cite{kuang:164101}, our computational
results for $G$ and $G^\prime$ do not seem to be affected by this
drawback of the model for metal. Instead, our results suggest that the
modeling of interfacial thermal transport behavior relies mainly on
the accuracy of the interaction descriptions between components
occupying the interfaces.

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

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}
\end{figure}

[MAY ADD COMPARISON OF AU SLAB WIDTHS]

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

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. Moderately lower coverages allow higher contact between
capping agent and solvent, and thus could further enhance the heat
transfer process.

Our measurement results, particularly of the UA models, agree with
available experimental data. This indicates that our force field
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. Of the two definitions for $G$, the discrete form
(Eq. \ref{discreteG}) was easier to use and gives out relatively
consistent results, while the derivative form (Eq. \ref{derivativeG})
is not as versatile. Although $G^\prime$ gives out comparable results
and follows similar trend with $G$ when measuring close to fully
covered or bare surfaces, the spatial resolution of $T$ profile is
limited for accurate computation of derivatives data.

Vlugt {\it et al.} has investigated the surface thiol structures for
nanocrystal gold and pointed out that they differs from those of the
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:	3749
Committed:	Mon Jul 25 19:11:33 2011 UTC (12 years, 11 months ago) by skuang
Content type:	application/x-tex
File size:	50279 byte(s)
Log Message:	add new references; some edits.