--- interfacial/interfacial.tex 2011/06/29 18:14:23 3729 +++ interfacial/interfacial.tex 2011/07/29 15:45:14 3755 @@ -23,8 +23,8 @@ \setlength{\belowcaptionskip}{30 pt} %\renewcommand\citemid{\ } % no comma in optional reference note -\bibpunct{[}{]}{,}{s}{}{;} -\bibliographystyle{aip} +\bibpunct{[}{]}{,}{n}{}{;} +\bibliographystyle{achemso} \begin{document} @@ -45,20 +45,22 @@ We have developed a Non-Isotropic Velocity Scaling alg \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. +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} @@ -71,98 +73,158 @@ leads to higher interfacial thermal transfer efficienc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \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. +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. -Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS) +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 more effective thermal exchange -between particles of different identities, and thus enables extensive -study of interfacial conductance. +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{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. +\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 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 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 quantity is not restricted by particle mass difference. +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 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. +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 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: +\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 - + 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. +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, 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: +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} + 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 +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| @@ -173,101 +235,103 @@ With the temperature profile obtained from simulations \label{derivativeG} \end{equation} -With the temperature profile obtained from simulations, one is able to +With temperature profiles obtained from simulation, one is able to approximate the first and second derivatives of $T$ with finite -difference method and thus calculate $G^\prime$. +difference methods and calculate $G^\prime$. In what follows, both +definitions have been used, and are compared in the results. -In what follows, both definitions are used for calculation and comparison. +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}. -[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}. +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]{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.} +\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{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. +\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. -[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] +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 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. +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 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. +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. -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 +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 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| +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 @@ -276,176 +340,329 @@ G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^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 various components. Therefore, force field -parameter descriptions are needed for interactions both between the -same type of particles and between particles of different species. +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}. +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 +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, 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] +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. -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] +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. -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. +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. -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} +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. -[TABULATED FORCE FIELD PARAMETERS NEEDED] +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. -[SURFACE RECONSTRUCTION PREVENTS SIMULATION TEMP TO GO HIGHER] +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} -[REARRANGEMENT NEEDED] -\subsection{Toluene Solvent} +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. -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. +\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 conductivity ($G$ and - $G^\prime$) values for the Au/butanethiol/toluene interface at - different temperatures using a range of energy fluxes.} - \begin{tabular}{cccc} + \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 - $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\ - (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\ + Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\ + (or bare surface) & model & (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 \\ + 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{AuThiolToluene} + \label{modelTest} \end{center} \end{minipage} \end{table*} -\subsection{Hexane Solvent} +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. -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. +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. -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] +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. -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. +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. -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. +\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} -%MAY ADD MORE DATA TO TABLE +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 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.} + $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}{cccccc} + \begin{tabular}{ccccccc} \hline\hline - Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\ - coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) & + $\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 - 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 \\ + 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} @@ -453,81 +670,224 @@ For the all-atom model, the liquid hexane phase was no \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. +\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 the Au/butanethiol/hexane interface - with all-atom model and different capping agent coverage at - 200K using a range of energy fluxes.} + $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}{cccc} + \begin{tabular}{ccccc} \hline\hline - Thiol & $J_z$ & $G$ & $G^\prime$ \\ - coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\ + $\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 - 0.0 & 0.95 & 28.5 & 27.2 \\ - & 1.88 & 30.3 & 28.9 \\ - 100.0 & 2.87 & 551 & 294 \\ - & 3.81 & 494 & 193 \\ + 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{AuThiolHexaneAA} + \label{AuThiolToluene} \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. +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 (upper panel) and for Au/thiol/hexane simulation in - all-atom model (lower panel).} -\label{vibration} + environments.} +\label{specAu} \end{figure} -% 600dpi, letter size. too large? +\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 +Dame. +\newpage \bibliography{interfacial}