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}
Due to the importance of heat flow in nanotechnology, interfacial
thermal conductance has been studied extensively both experimentally
and computationally.\cite{cahill:793} Unlike bulk materials, nanoscale
materials have a significant fraction of their atoms at interfaces,
and the chemical details of these interfaces govern the heat transfer
behavior. Furthermore, the interfaces are
heterogeneous (e.g. solid - liquid), which provides a challenge to
traditional methods developed for homogeneous systems. 

Experimentally, various interfaces have been investigated for their
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
solvent,\cite{doi:10.1021/jp8051888} and Juv\'e {\it et al.} studied
the cooling dynamics, which is controlled by thermal interface
resistence of glass-embedded metal
nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
normally considered barriers for heat transport, Alper {\it et al.}
suggested that specific ligands (capping agents) could completely
eliminate this barrier
($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}

Theoretical and computational models have also been used to study the
interfacial thermal transport in order to gain an understanding of
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 have yet to be studied with their
approach. The comparatively low thermal flux through interfaces is
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
the resulting gradient), given that the simulation methods being able
to effectively apply an unphysical flux in non-homogeneous systems.
Garde and coworkers\cite{garde:nl2005,garde:PhysRevLett2009} applied
this approach to various liquid interfaces and studied how thermal
conductance (or resistance) is dependent on chemistry details of
interfaces, e.g. hydrophobic and hydrophilic aqueous interfaces.

Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS)
algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
retains the desirable features of RNEMD (conservation of linear
momentum and total energy, compatibility with periodic boundary
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.

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

The non-isotropic velocity scaling (NIVS) \cite{kuang:164101} approach
to non-equilibrium MD simulations is able to impose a wide range of
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$)}

For an interface with relatively low interfacial conductance, and a
thermal flux between two distinct bulk regions, the regions on either
side of the interface rapidly come to a state in which the two phases
have relatively homogeneous (but distinct) temperatures. The
interfacial thermal conductivity $G$ can therefore be approximated as:
\begin{equation}
  G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
    \langle T_\mathrm{cold}\rangle \right)}
\label{lowG}
\end{equation}
where ${E_{total}}$ is the total imposed non-physical kinetic energy
transfer during the simulation and ${\langle T_\mathrm{hot}\rangle}$
and ${\langle T_\mathrm{cold}\rangle}$ are the average observed
temperature of the two separated phases.

When the interfacial conductance is {\it not} small, there are two
ways to define $G$. One 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}). This is
known as the Kapitza conductance, which is the inverse of the Kapitza
resistance.
\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 reaches 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 temperature profiles obtained from simulation, one is able to
approximate the first and second derivatives of $T$ with finite
difference methods and calculate $G^\prime$. In what follows, both
definitions have been used, and are compared in the results.

To investigate the interfacial conductivity at metal / solvent
interfaces, we have modeled a metal slab with its (111) surfaces
perpendicular to the $z$-axis of our simulation cells. The metal slab
has been prepared both with and without capping agents on the exposed
surface, and has been solvated with simple organic solvents, as
illustrated in Figure \ref{gradT}.

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

\begin{figure}
\includegraphics[width=\linewidth]{gradT}
\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.
Metal slabs of 6 or 11 layers of Au atoms were first equilibrated
under atmospheric pressure (1 atm) and 200K. After equilibration,
butanethiol capping agents were placed at three-fold hollow sites on
the Au(111) surfaces. These sites are either {\it fcc} or {\it
  hcp} sites, although Hase {\it et al.} found that they are
equivalent in a heat transfer process,\cite{hase:2010} so we did not
distinguish between these sites in our study. The maximum butanethiol
capacity on Au surface is $1/3$ of the total number of surface Au
atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
series of lower coverages was also prepared by eliminating
butanethiols from the higher coverage surface in a regular manner. The
lower coverages were prepared in order to study the relation between
coverage and interfacial conductance.

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

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

The simulation cells were not particularly extensive along the
$z$-axis, as a very long length scale for the thermal gradient may
cause excessively hot or cold temperatures in the middle of the
solvent region and lead to undesired phenomena such as solvent boiling
or freezing when a thermal flux is applied. Conversely, too few
solvent molecules would change the normal behavior of the liquid
phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
these extreme cases did not happen to our simulations. The spacing
between periodic images of the gold interfaces is $45 \sim 75$\AA.

The initial configurations generated are further equilibrated with the
$x$ and $y$ dimensions fixed, only allowing the $z$-length scale to
change. This is to ensure that the equilibration of liquid phase does
not affect the metal's crystalline structure. Comparisons were made
with simulations that allowed changes of $L_x$ and $L_y$ during NPT
equilibration. No substantial changes in the box geometry were noticed
in these simulations. After ensuring the liquid phase reaches
equilibrium at atmospheric pressure (1 atm), further equilibration was
carried out under canonical (NVT) and microcanonical (NVE) ensembles.

After the systems reach equilibrium, NIVS was used to impose an
unphysical thermal flux between the metal and the liquid phases. Most
of our simulations were done under an average temperature of
$\sim$200K. Therefore, thermal flux usually came from the metal to the
liquid so that the liquid has a higher temperature and would not
freeze due to lowered temperatures. After this induced temperature
gradient had stablized, the temperature profile of the simulation cell
was recorded. To do this, the simulation cell is devided evenly into
$N$ slabs along the $z$-axis. The average temperatures of each slab
are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
the same, the derivatives of $T$ with respect to slab number $n$ can
be directly used for $G^\prime$ calculations: \begin{equation}
  G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
         \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. Each
equilibration stage took a minimum of 100 ps, although in some cases,
longer equilibration stages were utilized.

\subsection{Force Field Parameters}
Our simulations include a number of chemically distinct components.
Figure \ref{demoMol} demonstrates the sites defined for both
United-Atom and All-Atom models of the organic solvent and capping
agents in our simulations. Force field parameters are needed for
interactions both between the same type of particles and between
particles of different species.

\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,TraPPE-UA.thiols}
  (UA) and \protect\cite{OPLSAA} (AA). Cross-interactions with the Au
  atoms are given in Table \ref{MnM}.}
\label{demoMol}
\end{figure}

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

For the two solvent molecules, {\it n}-hexane and toluene, two
different atomistic models were utilized. Both solvents were modeled
using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
for our UA solvent molecules. In these models, sites are located at
the carbon centers for alkyl groups. Bonding interactions, including
bond stretches and bends and torsions, were used for intra-molecular
sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones
potentials are used.

By eliminating explicit hydrogen atoms, the TraPPE-UA models are
simple and computationally efficient, while maintaining good accuracy.
However, the TraPPE-UA model for alkanes is known to predict a slighly
lower boiling point than experimental values. This is one of the
reasons we used a lower average temperature (200K) for our
simulations. If heat is transferred to the liquid phase during the
NIVS simulation, the liquid in the hot slab can actually be
substantially warmer than the mean temperature in the simulation. The
lower mean temperatures therefore prevent solvent boiling.

For UA-toluene, the non-bonded potentials between intermolecular sites
have a similar Lennard-Jones formulation. The toluene molecules were
treated as a single rigid body, so there was no need for
intramolecular interactions (including bonds, bends, or torsions) in
this solvent model.

Besides the TraPPE-UA models, AA models for both organic solvents are
included in our studies as well. The OPLS-AA\cite{OPLSAA} force fields
were used. For hexane, additional explicit hydrogen sites were
included. Besides bonding and non-bonded site-site interactions,
partial charges and the electrostatic interactions were added to each
CT and HC site. For toluene, a flexible model for the toluene molecule
was utilized which included bond, bend, torsion, and inversion
potentials to enforce ring planarity. 

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

To describe the interactions between metal (Au) and non-metal atoms,
we refer to an adsorption study of alkyl thiols on gold surfaces by
Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
Lennard-Jones form of potential parameters for the interaction between
Au and pseudo-atoms CH$_x$ and S based on a well-established and
widely-used effective potential of Hautman and Klein for the Au(111)
surface.\cite{hautman:4994} As our simulations require the gold slab
to be flexible to accommodate thermal excitation, the pair-wise form
of potentials they developed was used for our study.

The potentials developed from {\it ab initio} calculations by Leng
{\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
interactions between Au and aromatic C/H atoms in toluene. However,
the Lennard-Jones parameters between Au and other types of particles,
(e.g. AA alkanes) have not yet been established. For these
interactions, the Lorentz-Berthelot mixing rule can be used to derive
effective single-atom LJ parameters for the metal using the fit values
for toluene. These are then used to construct reasonable mixing
parameters for the interactions between the gold and other atoms.
Table \ref{MnM} summarizes the ``metal/non-metal'' parameters used in
our simulations.

\begin{table*}
  \begin{minipage}{\linewidth}
    \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*}


\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{Effects due to capping agent coverage}

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. 

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

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.

The similarity in the vibrational modes available to solvent and capping agent can be reduced by deuterating one of the two components.  Once either the hexanes or the butanethiols are deuterated, one can see a significantly lower $G$ and $G^\prime$ (Figure \ref{aahxntln}).  Compared to the AA model, the UA model yields more reasonable conductivity values with much higher computational efficiency.

\begin{figure}
\includegraphics[width=\linewidth]{aahxntln}
\caption{Spectra obtained for all-atom (AA) Au / butanethiol / solvent
  systems. When butanethiol is deuterated (lower left), its
  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}

For the Au / butanethiol / toluene interfaces, having the AA butanethiol deuterated did not yield a significant change in the 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.  The presence of extra degrees of freedom in the AA force field for toluene yields higher heat exchange rates between the two phases and results in a much higher conductivity than in the UA force field.

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

START HERE 

We have varied our protocol or other parameters of the simulations in order to investigate how these factors would affect the computation of $G$. 

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.

\subsubsection{Effects of applied flux}
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*}

\subsubsection{Effects due to average temperature}

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.


\section{Discussion}

\subsection{Capping agent acts as a vibrational coupler between solid
  and solvent phases}
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}


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

\subsubsection{Overlap of power spectrum}
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}

\subsubsection{Isotopic substitution and vibrational overlap}
A comparison of the results obtained from the two different organic
solvents can also provide useful information of the interfacial
thermal transport process. The deuterated hexane (UA) results do not
appear to be substantially different from those of normal hexane (UA),
given that butanethiol (UA) is non-deuterated for both solvents. The
UA models, even though they have eliminated C-H vibrational overlap,
still have significant overlap in the infrared spectra.  Because 
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}

\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{landman:1998,vlugt:cpc2007154}. This difference
might lead to change of interfacial thermal transport behavior as
well. To investigate this problem, an effective means to introduce
thermal flux and measure the corresponding thermal gradient is
desirable for simulating structures with spherical symmetry.

\section{Acknowledgments}
Support for this project was provided by the National Science
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:	3755
Committed:	Fri Jul 29 15:45:14 2011 UTC (12 years, 11 months ago) by skuang
Content type:	application/x-tex
File size:	50087 byte(s)
Log Message:	new edits