group/interfacial/interfacial.tex

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

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

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

\begin{document}

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

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

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}

We have developed a Non-Isotropic Velocity Scaling algorithm for
setting up and maintaining stable thermal gradients in non-equilibrium
molecular dynamics simulations. This approach effectively imposes
unphysical thermal flux even between particles of different
identities, conserves linear momentum and kinetic energy, and
minimally perturbs the velocity profile of a system when compared with
previous RNEMD methods. We have used this method to simulate thermal
conductance at metal / organic solvent interfaces both with and
without the presence of thiol-based capping agents.  We obtained
values comparable with experimental values, and observed significant
conductance enhancement with the presence of capping agents. Computed
power spectra indicate the acoustic impedance mismatch between metal
and liquid phase is greatly reduced by the capping agents and thus
leads to higher interfacial thermal transfer efficiency.

\end{abstract}

\newpage

%\narrowtext

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

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

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

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

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

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

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

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

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

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

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

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

[IMPOSE G DEFINITION INTO OUR SYSTEMS]
To facilitate the use of the above definitions in calculating $G$ and
$G^\prime$, we have a metal slab with its (111) surfaces perpendicular
to the $z$-axis of our simulation cells. With or withour capping
agents on the surfaces, the metal slab is solvated with organic
solvents, as illustrated in Figure \ref{demoPic}.

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

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

\begin{figure}
\includegraphics[width=\linewidth]{gradT}
\caption{The 1st and 2nd derivatives of temperature profile can be
  obtained with finite difference approximation.}
\label{gradT}
\end{figure}

\section{Computational Details}
\subsection{System Geometry}
In our simulations, Au is used to construct a metal slab with bare
(111) surface perpendicular to the $z$-axis. Different slab thickness
(layer numbers of Au) are simulated. This metal slab is first
equilibrated under normal pressure (1 atm) and a desired
temperature. After equilibration, butanethiol is used as the capping
agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
atoms in the butanethiol molecules would occupy the three-fold sites
of the surfaces, and the maximal butanethiol capacity on Au surface is
$1/3$ of the total number of surface Au atoms[CITATION]. A series of
different coverage surfaces is investigated in order to study the
relation between coverage and conductance.

[COVERAGE DISCRIPTION] However, since the interactions between surface
Au and butanethiol is non-bonded, the capping agent molecules are
allowed to migrate to an empty neighbor three-fold site during a
simulation. Therefore, the initial configuration would not severely
affect the sampling of a variety of configurations of the same
coverage, and the final conductance measurement would be an average
effect of these configurations explored in the simulations. [MAY NEED FIGURES]

After the modified Au-butanethiol surface systems are equilibrated
under canonical ensemble, Packmol\cite{packmol} is used to pack
organic solvent molecules in the previously vacuum part of the
simulation cells, which guarantees that short range repulsive
interactions do not disrupt the simulations. Two solvents are
investigated, one which has little vibrational overlap with the
alkanethiol and plane-like shape (toluene), and one which has similar
vibrational frequencies and chain-like shape ({\it n}-hexane). The
spacing filled by solvent molecules, i.e. the gap between periodically
repeated Au-butanethiol surfaces should be carefully chosen so that it
would not be too short to affect the liquid phase structure, nor too
long, leading to over cooling (freezing) or heating (boiling) when a
thermal flux is applied. In our simulations, this spacing is usually
$35 \sim 60$\AA.

The initial configurations generated by Packmol are further
equilibrated with the $x$ and $y$ dimensions fixed, only allowing
length scale change in $z$ dimension. This is to ensure that the
equilibration of liquid phase does not affect the metal crystal
structure in $x$ and $y$ dimensions. Further equilibration are run
under NVT and then NVE ensembles.

After the systems reach equilibrium, NIVS is implemented to impose a
periodic unphysical thermal flux between the metal and the liquid
phase. Most of our simulations are under an average temperature of
$\sim$200K. Therefore, this flux usually comes from the metal to the
liquid so that the liquid has a higher temperature and would not
freeze due to excessively low temperature. This induced temperature
gradient is stablized and the simulation cell is devided evenly into
N slabs along the $z$-axis and the temperatures of each slab are
recorded. When the slab width $d$ of each slab is the same, the
derivatives of $T$ with respect to slab number $n$ can be directly
used for $G^\prime$ calculations:
\begin{equation}
G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
         \Big/\left(\frac{\partial T}{\partial z}\right)^2
         = |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
         \Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
         = |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big|
         \Big/\left(\frac{\partial T}{\partial n}\right)^2
\label{derivativeG2}
\end{equation}

\subsection{Force Field Parameters}
Our simulations include various components. Therefore, force field
parameter descriptions are needed for interactions both between the
same type of particles and between particles of different species.

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

For both solvent molecules, straight chain {\it n}-hexane and aromatic
toluene, United-Atom (UA) and All-Atom (AA) models are used
respectively. The TraPPE-UA
parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
for our UA solvent molecules. In these models, pseudo-atoms are
located at the carbon centers for alkyl groups. By eliminating
explicit hydrogen atoms, these models are simple and computationally
efficient, while maintains good accuracy. However, the TraPPE-UA for
alkanes is known to predict a lower boiling point than experimental
values. Considering that after an unphysical thermal flux is applied
to a system, the temperature of ``hot'' area in the liquid phase would be
significantly higher than the average, to prevent over heating and
boiling of the liquid phase, the average temperature in our
simulations should be much lower than the liquid boiling point. [NEED MORE DISCUSSION]
For UA-toluene model, rigid body constraints are applied, so that the
benzene ring and the methyl-C(aromatic) bond are kept rigid. This
would save computational time.[MORE DETAILS NEEDED]

Besides the TraPPE-UA models, AA models for both organic solvents are
included in our studies as well. For hexane, the OPLS
all-atom\cite{OPLSAA} force field is used. [MORE DETAILS] 
For toluene, the United Force Field developed by Rapp\'{e} {\it et
  al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]

The capping agent in our simulations, the butanethiol molecules can
either use UA or AA model. The TraPPE-UA force fields includes
parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used in
our simulations corresponding to our TraPPE-UA models for solvent. 
 and All-Atom models [NEED CITATIONS]
However, the model choice (UA or AA) of capping agent can be different
from the solvent. Regardless of model choice, the force field
parameters for interactions between capping agent and solvent can be
derived using Lorentz-Berthelot Mixing Rule.

To describe the interactions between metal Au and non-metal capping
agent and solvent, we refer to Vlugt\cite{vlugt:cpc2007154} and derive
other interactions which are not yet finely parametrized. [can add
hautman and klein's paper here and more discussion; need to put
aromatic-metal interaction approximation here]\cite{doi:10.1021/jp034405s}

[TABULATED FORCE FIELD PARAMETERS NEEDED]


[SURFACE RECONSTRUCTION PREVENTS SIMULATION TEMP TO GO HIGHER]


\section{Results}
[REARRANGEMENT NEEDED]
\subsection{Toluene Solvent}

The results (Table \ref{AuThiolToluene}) show a
significant conductance enhancement compared to the gold/water
interface without capping agent and agree with available experimental
data. This indicates that the metal-metal potential, though not
predicting an accurate bulk metal thermal conductivity, does not
greatly interfere with the simulation of the thermal conductance
behavior across a non-metal interface. The solvent model is not
particularly volatile, so the simulation cell does not expand
significantly under higher temperature. We did not observe a
significant conductance decrease when the temperature was increased to
300K. The results show that the two definitions used for $G$ yield
comparable values, though $G^\prime$ tends to be smaller.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      \caption{Computed interfacial thermal conductivity ($G$ and
        $G^\prime$) values for the Au/butanethiol/toluene interface at
        different temperatures using a range of energy fluxes.}
      
      \begin{tabular}{cccc}
        \hline\hline
        $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
        (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        200 & 1.86 & 180 & 135 \\
            & 2.15 & 204 & 113 \\
            & 3.93 & 175 & 114 \\
        300 & 1.91 & 143 & 125 \\
            & 4.19 & 134 & 113 \\
        \hline\hline
      \end{tabular}
      \label{AuThiolToluene}
    \end{center}
  \end{minipage}
\end{table*}

\subsection{Hexane Solvent}

Using the united-atom model, different coverages of capping agent,
temperatures of simulations and numbers of solvent molecules were all
investigated and Table \ref{AuThiolHexaneUA} shows the results of
these computations. The number of hexane molecules in our simulations
does not affect the calculations significantly. However, 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, while too
few solvent molecules would change the normal behavior of the liquid
phase. Our $N_{hexane}$ values were chosen to ensure that these
extreme cases did not happen to our simulations.

Table \ref{AuThiolHexaneUA} enables direct comparison between
different coverages of capping agent, when other system parameters are
held constant. With high coverage of butanethiol on the gold surface,
the interfacial thermal conductance is enhanced
significantly. Interestingly, a slightly lower butanethiol coverage
leads to a moderately higher conductivity. This is probably due to
more solvent/capping agent contact when butanethiol molecules are
not densely packed, which enhances the interactions between the two
phases and lowers the thermal transfer barrier of this interface.
% [COMPARE TO AU/WATER IN PAPER]

It is also noted that the overall simulation temperature is another
factor that affects the interfacial thermal conductance. One
possibility of this effect may be rooted in the decrease in density of
the liquid phase. We observed that when the average temperature
increases from 200K to 250K, the bulk hexane density becomes lower
than experimental value, as the system is equilibrated under NPT
ensemble. This leads to lower contact between solvent and capping
agent, and thus lower conductivity.

Conductivity values are more difficult to obtain under higher
temperatures. This is because the Au surface tends to undergo
reconstructions in relatively high temperatures. Surface Au atoms can
migrate outward to reach higher Au-S contact; and capping agent
molecules can be embedded into the surface Au layer due to the same
driving force. This phenomenon agrees with experimental
results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. A surface
fully covered in capping agent is more susceptible to reconstruction,
possibly because fully coverage prevents other means of capping agent
relaxation, such as migration to an empty neighbor three-fold site.

%MAY ADD MORE DATA TO TABLE
\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      \caption{Computed interfacial thermal conductivity ($G$ and
        $G^\prime$) values for the Au/butanethiol/hexane interface
        with united-atom model and different capping agent coverage
        and solvent molecule numbers at different temperatures using a
        range of energy fluxes.}
      
      \begin{tabular}{cccccc}
        \hline\hline
        Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\
        coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) &
        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        0.0   & 200 & 200 & 0.96 & 43.3 & 42.7 \\
              &     &     & 1.91 & 45.7 & 42.9 \\
              &     & 166 & 0.96 & 43.1 & 53.4 \\
        88.9  & 200 & 166 & 1.94 & 172  & 108  \\
        100.0 & 250 & 200 & 0.96 & 81.8 & 67.0 \\
              &     & 166 & 0.98 & 79.0 & 62.9 \\
              &     &     & 1.44 & 76.2 & 64.8 \\
              & 200 & 200 & 1.92 & 129  & 87.3 \\
              &     &     & 1.93 & 131  & 77.5 \\
              &     & 166 & 0.97 & 115  & 69.3 \\
              &     &     & 1.94 & 125  & 87.1 \\
        \hline\hline
      \end{tabular}
      \label{AuThiolHexaneUA}
    \end{center}
  \end{minipage}
\end{table*}

For the all-atom model, the liquid hexane phase was not stable under NPT
conditions. Therefore, the simulation length scale parameters are
adopted from previous equilibration results of the united-atom model
at 200K. Table \ref{AuThiolHexaneAA} shows the results of these
simulations. The conductivity values calculated with full capping
agent coverage are substantially larger than observed in the
united-atom model, and is even higher than predicted by
experiments. It is possible that our parameters for metal-non-metal
particle interactions lead to an overestimate of the interfacial
thermal conductivity, although the active C-H vibrations in the
all-atom model (which should not be appreciably populated at normal
temperatures) could also account for this high conductivity. The major
thermal transfer barrier of Au/butanethiol/hexane interface is between
the liquid phase and the capping agent, so extra degrees of freedom
such as the C-H vibrations could enhance heat exchange between these
two phases and result in a much higher conductivity.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed interfacial thermal conductivity ($G$ and
        $G^\prime$) values for the Au/butanethiol/hexane interface
        with all-atom model and different capping agent coverage at
        200K using a range of energy fluxes.}
      
      \begin{tabular}{cccc}
        \hline\hline
        Thiol & $J_z$ & $G$ & $G^\prime$ \\
        coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        0.0   & 0.95 & 28.5 & 27.2 \\
              & 1.88 & 30.3 & 28.9 \\
        100.0 & 2.87 & 551  & 294  \\
              & 3.81 & 494  & 193  \\
        \hline\hline
      \end{tabular}
      \label{AuThiolHexaneAA}
    \end{center}
  \end{minipage}
\end{table*}

%subsubsection{Vibrational spectrum study on conductance mechanism}
To investigate the mechanism of this interfacial thermal conductance,
the vibrational spectra of various gold systems were obtained and are
shown as in the upper panel of Fig. \ref{vibration}. To obtain these
spectra, one first runs a simulation in the NVE ensemble and collects
snapshots of configurations; these configurations are used to compute
the velocity auto-correlation functions, which is used to construct a
power spectrum via a Fourier transform. The gold surfaces covered by
butanethiol molecules exhibit an additional peak observed at a
frequency of $\sim$170cm$^{-1}$, which is attributed to the vibration
of the S-Au bond. This vibration enables efficient thermal transport
from surface Au atoms to the capping agents. Simultaneously, as shown
in the lower panel of Fig. \ref{vibration}, the large overlap of the
vibration spectra of butanethiol and hexane in the all-atom model,
including the C-H vibration, also suggests high thermal exchange
efficiency. The combination of these two effects produces the drastic
interfacial thermal conductance enhancement in the all-atom model.

\begin{figure}
\includegraphics[width=\linewidth]{vibration}
\caption{Vibrational spectra obtained for gold in different
  environments (upper panel) and for Au/thiol/hexane simulation in
  all-atom model (lower panel).}
\label{vibration}
\end{figure}
% 600dpi, letter size. too large?


\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:	3729
Committed:	Wed Jun 29 18:14:23 2011 UTC (13 years ago) by skuang
Content type:	application/x-tex
File size:	25043 byte(s)
Log Message:	add more data, citations, some work on computational details.