--- interfacial/interfacial.tex 2011/06/24 16:59:37 3727 +++ interfacial/interfacial.tex 2011/07/12 22:11:11 3737 @@ -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,66 +73,103 @@ 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. +experimentally and computationally\cite{cahill:793}, due to its +importance in nanoscale science and technology. Reliability of +nanoscale devices depends on their thermal transport +properties. Unlike bulk homogeneous materials, nanoscale materials +features significant presence of interfaces, and these interfaces +could dominate the heat transfer behavior of these +materials. Furthermore, these materials are generally heterogeneous, +which challenges traditional research methods for homogeneous +systems. +Heat conductance of molecular and nano-scale interfaces will be +affected by the chemical details of the surface. Experimentally, +various interfaces have been investigated for their thermal +conductance properties. Wang {\it et al.} studied heat transport +through long-chain hydrocarbon monolayers on gold substrate at +individual molecular level\cite{Wang10082007}; Schmidt {\it et al.} +studied the role of CTAB on thermal transport between gold nanorods +and solvent\cite{doi:10.1021/jp8051888}; Juv\'e {\it et al.} studied +the cooling dynamics, which is controlled by thermal interface +resistence of glass-embedded metal +nanoparticles\cite{PhysRevB.80.195406}. Although interfaces are +commonly barriers for heat transport, Alper {\it et al.} suggested +that specific ligands (capping agents) could completely eliminate this +barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}. + +Theoretical and computational models have also been used to study the +interfacial thermal transport in order to gain an understanding of +this phenomena at the molecular level. Recently, Hase and coworkers +employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to +study thermal transport from hot Au(111) substrate to a self-assembled +monolayer of alkylthiolate with relatively long chain (8-20 carbon +atoms)\cite{hase:2010,hase:2011}. However, ensemble averaged +measurements for heat conductance of interfaces between the capping +monolayer on Au and a solvent phase has yet to be studied. +The relatively low thermal flux through interfaces is +difficult to measure with Equilibrium MD or forward NEMD simulation +methods. Therefore, the Reverse NEMD (RNEMD) methods would have the +advantage of having this difficult to measure flux known when studying +the thermal transport across interfaces, given that the simulation +methods being able to effectively apply an unphysical flux in +non-homogeneous systems. + Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS) algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm retains the desirable features of RNEMD (conservation of linear momentum and total energy, compatibility with periodic boundary conditions) while establishing true thermal distributions in each of -the two slabs. Furthermore, it allows 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 this phenomena was investigated. + +[MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG] + \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 +\subsection{Imposd-Flux Methods in MD Simulations} +For systems with low interfacial conductivity one must have a method +capable of generating relatively small fluxes, compared to those +required for bulk conductivity. This requirement makes the calculation +even more difficult for those slowly-converging equilibrium +methods\cite{Viscardy:2007lq}. +Forward methods impose gradient, but in interfacail conditions it is +not clear what behavior to impose at the boundary... + Imposed-flux reverse 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 thermal response becomes easier to +measure than the flux. Although M\"{u}ller-Plathe's original momentum +swapping approach can be used for exchanging energy between particles +of different identity, the kinetic energy transfer efficiency is +affected by the mass difference between the particles, which limits +its application on heterogeneous interfacial systems. -The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach 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 @@ -148,13 +187,14 @@ When the interfacial conductance is {\it not} small, t 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$. +When the interfacial conductance is {\it not} small, there are two +ways 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: +One way is to assume the temperature is discrete on the two sides of +the interface. $G$ can be calculated using the applied thermal flux +$J$ and the maximum temperature difference measured along the thermal +gradient max($\Delta T$), which occurs at the Gibbs deviding surface, +as: \begin{equation} G=\frac{J}{\Delta T} \label{discreteG} @@ -175,16 +215,16 @@ difference method and thus calculate $G^\prime$. 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$. +difference methods and thus calculate $G^\prime$. -In what follows, both definitions are used for calculation and comparison. +In what follows, both definitions have been used for calculation and +are compared in the results. -[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}. +To compare the above definitions ($G$ and $G^\prime$), we have modeled +a metal slab with its (111) surfaces perpendicular to the $z$-axis of +our simulation cells. Both with and withour capping agents on the +surfaces, the metal slab is solvated with simple organic solvents, as +illustrated in Figure \ref{demoPic}. \begin{figure} \includegraphics[width=\linewidth]{demoPic} @@ -193,16 +233,14 @@ With a simulation cell setup following the above manne \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. +With the simulation cell described above, we are able to equilibrate +the system and impose an unphysical thermal flux between the liquid +and the metal phase using the NIVS algorithm. By periodically applying +the unphysical flux, we are able to obtain a temperature profile and +its spatial derivatives. These quantities enable the evaluation of the +interfacial thermal conductance of a surface. Figure \ref{gradT} is an +example how those applied thermal fluxes can be used to obtain the 1st +and 2nd derivatives of the temperature profile. \begin{figure} \includegraphics[width=\linewidth]{gradT} @@ -212,44 +250,56 @@ obtain the 1st and 2nd derivatives of the temperature \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. Different slab thickness (layer numbers of Au) were +simulated. Metal slabs were first equilibrated under atmospheric +pressure (1 atm) and a desired temperature (e.g. 200K). After +equilibration, butanethiol capping agents were placed at three-fold +sites on the Au(111) surfaces. The maximum butanethiol capacity on Au +surface is $1/3$ of the total number of surface Au +atoms\cite{vlugt:cpc2007154}. A series of different coverages was +investigated 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 would not noticeably affect the sampling of a +variety of configurations of the same coverage, and the final +conductance measurement would be an average effect of these +configurations explored in the simulations. [MAY NEED 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 -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 modified Au-butanethiol surface systems were equilibrated +under canonical ensemble, organic solvent molecules were packed in the +previously empty part of the simulation cells\cite{packmol}. Two +solvents were investigated, one which has little vibrational overlap +with the alkanethiol and a planar shape (toluene), and one which has +similar vibrational frequencies and chain-like shape ({\it n}-hexane). +The space filled by solvent molecules, i.e. the gap between +periodically repeated Au-butanethiol surfaces should be carefully +chosen. A very long length scale for the thermal gradient axis ($z$) +may cause excessively hot or cold temperatures in the middle of the +solvent region and lead to undesired phenomena such as solvent boiling +or freezing when a thermal flux is applied. Conversely, too few +solvent molecules would change the normal behavior of the liquid +phase. Therefore, our $N_{solvent}$ values were chosen to ensure that +these extreme cases did not happen to our simulations. And the +corresponding spacing is usually $35 \sim 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 have this flux from the metal to the +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 @@ -267,68 +317,254 @@ G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^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 +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}. -Straight chain {\it n}-hexane and aromatic toluene are respectively -used as solvents. For hexane, both United-Atom\cite{TraPPE-UA.alkanes} -and All-Atom\cite{OPLSAA} force fields are used for comparison; for -toluene, United-Atom\cite{TraPPE-UA.alkylbenzenes} force fields are -used with rigid body constraints applied. (maybe needs more details -about rigid body) +Figure \ref{demoMol} demonstrates how we name our pseudo-atoms of the +organic solvent molecules in our simulations. -Buatnethiol molecules are used as capping agent for some of our -simulations. United-Atom\cite{TraPPE-UA.thiols} and All-Atom models -are respectively used corresponding to the force field type of -solvent. +\begin{figure} +\includegraphics[width=\linewidth]{demoMol} +\caption{Denomination of atoms or pseudo-atoms in our simulations: a) + UA-hexane; b) AA-hexane; c) UA-toluene; d) AA-toluene.} +\label{demoMol} +\end{figure} +For both solvent molecules, straight chain {\it n}-hexane and aromatic +toluene, United-Atom (UA) and All-Atom (AA) models are used +respectively. The TraPPE-UA +parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used +for our UA solvent molecules. In these models, pseudo-atoms are +located at the carbon centers for alkyl groups. By eliminating +explicit hydrogen atoms, these models are simple and computationally +efficient, while maintains good accuracy. However, the TraPPE-UA for +alkanes is known to predict a lower boiling point than experimental +values. Considering that after an unphysical thermal flux is applied +to a system, the temperature of ``hot'' area in the liquid phase would be +significantly higher than the average, to prevent over heating and +boiling of the liquid phase, the average temperature in our +simulations should be much lower than the liquid boiling point. [MORE DISCUSSION] +For UA-toluene model, rigid body constraints are applied, so that the +benzene ring and the methyl-CRar bond are kept rigid. This would save +computational time.[MORE DETAILS] + +Besides the TraPPE-UA models, AA models for both organic solvents are +included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA} +force field is used. [MORE DETAILS] +For toluene, the United Force Field developed by Rapp\'{e} {\it et + al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS] + +The capping agent in our simulations, the butanethiol molecules can +either use UA or AA model. The TraPPE-UA force fields includes +parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for +UA butanethiol model in our simulations. The OPLS-AA also provides +parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111) +surfaces do not have the hydrogen atom bonded to sulfur. To adapt this +change and derive suitable parameters for butanethiol adsorbed on +Au(111) surfaces, we adopt the S parameters from Luedtke and +Landman\cite{landman:1998} and modify parameters for its neighbor C +atom for charge balance in the molecule. Note that the model choice +(UA or AA) of capping agent can be different from the +solvent. Regardless of model choice, the force field parameters for +interactions between capping agent and solvent can be derived using +Lorentz-Berthelot Mixing Rule:[EQN'S] + + 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 parametrized in their work. (can add -hautman and klein's paper here and more discussion; need to put -aromatic-metal interaction approximation here) +agent and solvent particles, we refer to an adsorption study of alkyl +thiols on gold surfaces by Vlugt {\it et + al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones +form of potential parameters for the interaction between Au and +pseudo-atoms CH$_x$ and S based on a well-established and widely-used +effective potential of Hautman and Klein\cite{hautman:4994} for the +Au(111) surface. As our simulations require the gold lattice slab to +be non-rigid so that it could accommodate kinetic energy for thermal +transport study purpose, the pair-wise form of potentials is +preferred. -[TABULATED FORCE FIELD PARAMETERS NEEDED] +Besides, the potentials developed from {\it ab initio} calculations by +Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the +interactions between Au and aromatic C/H atoms in toluene.[MORE DETAILS] -\section{Results} -\subsection{Toluene Solvent} +However, the Lennard-Jones parameters between Au and other types of +particles in our simulations are not yet well-established. For these +interactions, we attempt to derive their parameters using the Mixing +Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters +for Au is first extracted from the Au-CH$_x$ parameters by applying +the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM'' +parameters in our simulations. -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{Lennard-Jones parameters for Au-non-Metal + interactions in our simulations.} + + \begin{tabular}{ccc} + \hline\hline + Non-metal atom & $\sigma$ & $\epsilon$ \\ + (or pseudo-atom) & \AA & kcal/mol \\ + \hline + S & 2.40 & 8.465 \\ + CH3 & 3.54 & 0.2146 \\ + CH2 & 3.54 & 0.1749 \\ + CT3 & 3.365 & 0.1373 \\ + CT2 & 3.365 & 0.1373 \\ + CTT & 3.365 & 0.1373 \\ + HC & 2.865 & 0.09256 \\ + CHar & 3.4625 & 0.1680 \\ + CRar & 3.555 & 0.1604 \\ + CA & 3.173 & 0.0640 \\ + HA & 2.746 & 0.0414 \\ + \hline\hline + \end{tabular} + \label{MnM} + \end{center} + \end{minipage} +\end{table*} + +\section{Results and Discussions} +[MAY HAVE A BRIEF SUMMARY] +\subsection{How Simulation Parameters Affects $G$} +[MAY NOT PUT AT FIRST] +We have varied our protocol or other parameters of the simulations in +order to investigate how these factors would affect the measurement of +$G$'s. It turned out that while some of these parameters would not +affect the results substantially, some other changes to the +simulations would have a significant impact on the measurement +results. + +In some of our simulations, we allowed $L_x$ and $L_y$ to change +during equilibrating the liquid phase. Due to the stiffness of the Au +slab, $L_x$ and $L_y$ would not change noticeably after +equilibration. Although $L_z$ could fluctuate $\sim$1\% after a system +is fully equilibrated in the NPT ensemble, this fluctuation, as well +as those comparably smaller to $L_x$ and $L_y$, would not be magnified +on the calculated $G$'s, as shown in Table \ref{AuThiolHexaneUA}. This +insensivity to $L_i$ fluctuations allows reliable measurement of $G$'s +without the necessity of extremely cautious equilibration process. + +As stated in our computational details, the spacing filled with +solvent molecules can be chosen within a range. This allows some +change of solvent molecule numbers for the same Au-butanethiol +surfaces. We did this study on our Au-butanethiol/hexane +simulations. Nevertheless, the results obtained from systems of +different $N_{hexane}$ did not indicate that the measurement of $G$ is +susceptible to this parameter. For computational efficiency concern, +smaller system size would be preferable, given that the liquid phase +structure is not affected. + +Our NIVS algorithm allows change of unphysical thermal flux both in +direction and in quantity. This feature extends our investigation of +interfacial thermal conductance. However, the magnitude of this +thermal flux is not arbitary if one aims to obtain a stable and +reliable thermal gradient. A temperature profile would be +substantially affected by noise when $|J_z|$ has a much too low +magnitude; while an excessively large $|J_z|$ that overwhelms the +conductance capacity of the interface would prevent a thermal gradient +to reach a stablized steady state. NIVS has the advantage of allowing +$J$ to vary in a wide range such that the optimal flux range for $G$ +measurement can generally be simulated by the algorithm. Within the +optimal range, we were able to study how $G$ would change according to +the thermal flux across the interface. For our simulations, we denote +$J_z$ to be positive when the physical thermal flux is from the liquid +to metal, and negative vice versa. The $G$'s measured under different +$J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These +results do not suggest that $G$ is dependent on $J_z$ within this flux +range. The linear response of flux to thermal gradient simplifies our +investigations in that we can rely on $G$ measurement with only a +couple $J_z$'s and do not need to test a large series of fluxes. + +%ADD MORE TO TABLE \begin{table*} \begin{minipage}{\linewidth} \begin{center} \caption{Computed interfacial thermal conductivity ($G$ and - $G^\prime$) values for the Au/butanethiol/toluene interface 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.} + \begin{tabular}{cccccccc} + \hline\hline + $\langle T\rangle$ & & $L_x$ & $L_y$ & $L_z$ & $J_z$ & + $G$ & $G^\prime$ \\ + (K) & $N_{hexane}$ & \multicolumn{3}{c}{(\AA)} & (GW/m$^2$) & + \multicolumn{2}{c}{(MW/m$^2$/K)} \\ + \hline + 200 & 266 & 29.86 & 25.80 & 113.1 & -0.96 & + 102() & 80.0() \\ + & 200 & 29.84 & 25.81 & 93.9 & 1.92 & + 129() & 87.3() \\ + & & 29.84 & 25.81 & 95.3 & 1.93 & + 131() & 77.5() \\ + & 166 & 29.84 & 25.81 & 85.7 & 0.97 & + 115() & 69.3() \\ + & & & & & 1.94 & + 125() & 87.1() \\ + 250 & 200 & 29.84 & 25.87 & 106.8 & 0.96 & + 81.8() & 67.0() \\ + & 166 & 29.87 & 25.84 & 94.8 & 0.98 & + 79.0() & 62.9() \\ + & & 29.84 & 25.85 & 95.0 & 1.44 & + 76.2() & 64.8() \\ + \hline\hline + \end{tabular} + \label{AuThiolHexaneUA} + \end{center} + \end{minipage} +\end{table*} + +Furthermore, we also attempted to increase system average temperatures +to above 200K. These simulations are first equilibrated in the NPT +ensemble under normal pressure. As stated above, the TraPPE-UA model +for hexane tends to predict a lower boiling point. In our simulations, +hexane had diffculty to remain in liquid phase when NPT equilibration +temperature is higher than 250K. Additionally, the equilibrated liquid +hexane density under 250K becomes lower than experimental value. This +expanded liquid phase leads to lower contact between hexane and +butanethiol as well.[MAY NEED FIGURE] And this reduced contact would +probably be accountable for a lower interfacial thermal conductance, +as shown in Table \ref{AuThiolHexaneUA}. + +A similar study for TraPPE-UA toluene agrees with the above result as +well. Having a higher boiling point, toluene tends to remain liquid in +our simulations even equilibrated under 300K in NPT +ensembles. Furthermore, the expansion of the toluene liquid phase is +not as significant as that of the hexane. This prevents severe +decrease of liquid-capping agent contact and the results (Table +\ref{AuThiolToluene}) show only a slightly decreased interface +conductance. Therefore, solvent-capping agent contact should play an +important role in the thermal transport process across the interface +in that higher degree of contact could yield increased conductance. + +[ADD Lxyz AND ERROR ESTIMATE TO TABLE] +\begin{table*} + \begin{minipage}{\linewidth} + \begin{center} + \caption{Computed interfacial thermal conductivity ($G$ and + $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene + interface at different temperatures using a range of energy + fluxes.} + \begin{tabular}{cccc} \hline\hline $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\ (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\ \hline - 200 & 1.86 & 180 & 135 \\ - & 2.15 & 204 & 113 \\ - & 3.93 & 175 & 114 \\ - 300 & 1.91 & 143 & 125 \\ - & 4.19 & 134 & 113 \\ + 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} @@ -336,145 +572,252 @@ comparable values, though $G^\prime$ tends to be small \end{minipage} \end{table*} -\subsection{Hexane Solvent} +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. -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. +However, when the surface is not completely covered by butanethiols, +the simulated system is more resistent to the reconstruction +above. Our Au-butanethiol/toluene system did not see this phenomena +even at $\langle T\rangle\sim$300K. The Au(111) surfaces have a 90\% coverage of +butanethiols and have empty three-fold sites. These empty sites could +help prevent surface reconstruction in that they provide other means +of capping agent relaxation. It is observed that butanethiols can +migrate to their neighbor empty sites during a simulation. Therefore, +we were able to obtain $G$'s for these interfaces even at a relatively +high temperature without being affected by surface reconstructions. -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] +\subsection{Influence of Capping Agent Coverage on $G$} +To investigate the influence of butanethiol coverage on interfacial +thermal conductance, a series of different coverage Au-butanethiol +surfaces is prepared and solvated with various organic +molecules. These systems are then equilibrated and their interfacial +thermal conductivity are measured with our NIVS algorithm. Table +\ref{tlnUhxnUhxnD} lists these results for direct comparison between +different coverages of butanethiol. To study the isotope effect in +interfacial thermal conductance, deuterated UA-hexane is included as +well. -It 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. +It turned out that with partial covered butanethiol on the Au(111) +surface, the derivative definition for $G$ (Eq. \ref{derivativeG}) has +difficulty to apply, due to the difficulty in locating the maximum of +change of $\lambda$. Instead, the discrete definition +(Eq. \ref{discreteG}) is easier to apply, as max($\Delta T$) can still +be well-defined. Therefore, $G$'s (not $G^\prime$) are used for this +section. -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. +From Table \ref{tlnUhxnUhxnD}, one can see the significance of the +presence of capping agents. Even when a fraction of the Au(111) +surface sites are covered with butanethiols, the conductivity would +see an enhancement by at least a factor of 3. This indicates the +important role cappping agent is playing for thermal transport +phenomena on metal/organic solvent surfaces. -%MAY ADD MORE DATA TO TABLE +Interestingly, as one could observe from our results, the maximum +conductance enhancement (largest $G$) happens while the surfaces are +about 75\% covered with butanethiols. This again indicates that +solvent-capping agent contact has an important role of the thermal +transport process. Slightly lower butanethiol coverage allows small +gaps between butanethiols to form. And these gaps could be filled with +solvent molecules, which acts like ``heat conductors'' on the +surface. The higher degree of interaction between these solvent +molecules and capping agents increases the enhancement effect and thus +produces a higher $G$ than densely packed butanethiol arrays. However, +once this maximum conductance enhancement is reached, $G$ decreases +when butanethiol coverage continues to decrease. Each capping agent +molecule reaches its maximum capacity for thermal +conductance. Therefore, even higher solvent-capping agent contact +would not offset this effect. Eventually, when butanethiol coverage +continues to decrease, solvent-capping agent contact actually +decreases with the disappearing of butanethiol molecules. In this +case, $G$ decrease could not be offset but instead accelerated. + +A comparison of the results obtained from differenet organic solvents +can also provide useful information of the interfacial thermal +transport process. The deuterated hexane (UA) results do not appear to +be much different from those of normal hexane (UA), given that +butanethiol (UA) is non-deuterated for both solvents. These UA model +studies, even though eliminating C-H vibration samplings, still have +C-C vibrational frequencies different from each other. However, these +differences in the infrared range do not seem to produce an observable +difference for the results of $G$. [MAY NEED FIGURE] + +Furthermore, results for rigid body toluene solvent, as well as other +UA-hexane solvents, are reasonable within the general experimental +ranges[CITATIONS]. This suggests that explicit hydrogen might not be a +required factor for modeling thermal transport phenomena of systems +such as Au-thiol/organic solvent. + +However, results for Au-butanethiol/toluene do not show an identical +trend with those for Au-butanethiol/hexane in that $G$'s remain at +approximately the same magnitue when butanethiol coverage differs from +25\% to 75\%. This might be rooted in the molecule shape difference +for plane-like toluene and chain-like {\it n}-hexane. Due to this +difference, toluene molecules have more difficulty in occupying +relatively small gaps among capping agents when their coverage is not +too low. Therefore, the solvent-capping agent contact may keep +increasing until the capping agent coverage reaches a relatively low +level. This becomes an offset for decreasing butanethiol molecules on +its effect to the process of interfacial thermal transport. Thus, one +can see a plateau of $G$ vs. butanethiol coverage in our results. + +[NEED ERROR ESTIMATE, MAY ALSO PUT J HERE] \begin{table*} \begin{minipage}{\linewidth} \begin{center} - \caption{Computed interfacial thermal conductivity ($G$ 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.} + \caption{Computed interfacial thermal conductivity ($G$) values + for the Au-butanethiol/solvent interface with various UA + models and different capping agent coverages at $\langle + T\rangle\sim$200K using certain energy flux respectively.} - \begin{tabular}{cccccc} + \begin{tabular}{cccc} \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)} \\ + Thiol & \multicolumn{3}{c}{$G$(MW/m$^2$/K)} \\ + coverage (\%) & hexane & hexane(D) & toluene \\ \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 \\ + 0.0 & 46.5() & 43.9() & 70.1() \\ + 25.0 & 151() & 153() & 249() \\ + 50.0 & 172() & 182() & 214() \\ + 75.0 & 242() & 229() & 244() \\ + 88.9 & 178() & - & - \\ + 100.0 & 137() & 153() & 187() \\ \hline\hline \end{tabular} - \label{AuThiolHexaneUA} + \label{tlnUhxnUhxnD} \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. +\subsection{Influence of Chosen Molecule Model on $G$} +[MAY COMBINE W MECHANISM STUDY] +In addition to UA solvent/capping agent models, AA models are included +in our simulations as well. Besides simulations of the same (UA or AA) +model for solvent and capping agent, different models can be applied +to different components. Furthermore, regardless of models chosen, +either the solvent or the capping agent can be deuterated, similar to +the previous section. Table \ref{modelTest} summarizes the results of +these studies. + +[MORE DATA; ERROR ESTIMATE] \begin{table*} \begin{minipage}{\linewidth} \begin{center} \caption{Computed interfacial thermal conductivity ($G$ and - $G^\prime$) values for 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 interfaces using various models for + solvent and capping agent (or without capping agent) at + $\langle T\rangle\sim$200K.} - \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)} \\ + Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\ + (or bare surface) & model & (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 \\ + UA & AA hexane & 1.94 & 135() & 129() \\ + & & 2.86 & 126() & 115() \\ + & AA toluene & 1.89 & 200() & 149() \\ + AA & UA hexane & 1.94 & 116() & 129() \\ + & AA hexane & 3.76 & 451() & 378() \\ + & & 4.71 & 432() & 334() \\ + & AA toluene & 3.79 & 487() & 290() \\ + AA(D) & UA hexane & 1.94 & 158() & 172() \\ + bare & AA hexane & 0.96 & 31.0() & 29.4() \\ \hline\hline \end{tabular} - \label{AuThiolHexaneAA} + \label{modelTest} \end{center} \end{minipage} \end{table*} -%subsubsection{Vibrational spectrum study on conductance mechanism} +To facilitate direct comparison, the same system with differnt models +for different components uses the same length scale for their +simulation cells. Without the presence of capping agent, using +different models for hexane yields similar results for both $G$ and +$G^\prime$, and these two definitions agree with eath other very +well. This indicates very weak interaction between the metal and the +solvent, and is a typical case for acoustic impedance mismatch between +these two phases. + +As for Au(111) surfaces completely covered by butanethiols, the choice +of models for capping agent and solvent could impact the measurement +of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane +interfaces, using AA model for both butanethiol and hexane yields +substantially higher conductivity values than using UA model for at +least one component of the solvent and capping agent, which exceeds +the upper bond of experimental value range. This is probably due to +the classically treated C-H vibrations in the AA model, which should +not be appreciably populated at normal temperatures. In comparison, +once either the hexanes or the butanethiols are deuterated, one can +see a significantly lower $G$ and $G^\prime$. In either of these +cases, the C-H(D) vibrational overlap between the solvent and the +capping agent is removed. [MAY NEED FIGURE] Conclusively, the +improperly treated C-H vibration in the AA model produced +over-predicted results accordingly. Compared to the AA model, the UA +model yields more reasonable results with higher computational +efficiency. + +However, for Au-butanethiol/toluene interfaces, having the AA +butanethiol deuterated did not yield a significant change in the +measurement results. +. , so extra degrees of freedom +such as the C-H vibrations could enhance heat exchange between these +two phases and result in a much higher conductivity. + + +Although the QSC model for Au is known to predict an overly low value +for bulk metal gold conductivity[CITE NIVSRNEMD], our computational +results for $G$ and $G^\prime$ do not seem to be affected by this +drawback of the model for metal. Instead, the modeling of interfacial +thermal transport behavior relies mainly on an accurate description of +the interactions between components occupying the interfaces. + +\subsection{Mechanism of Interfacial Thermal Conductance Enhancement + by Capping Agent} +%OR\subsection{Vibrational spectrum study on conductance mechanism} + +[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL, EQN'S] + To investigate the mechanism of this interfacial thermal conductance, the vibrational spectra of various gold systems were obtained and are shown as in the upper panel of Fig. \ref{vibration}. To obtain these spectra, one first runs a simulation in the NVE ensemble and collects snapshots of configurations; these configurations are used to compute the velocity auto-correlation functions, which is used to construct a -power spectrum via a Fourier transform. The gold surfaces covered by -butanethiol molecules 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. +power spectrum via a Fourier transform. + The gold surfaces covered by +butanethiol molecules, compared to bare gold surfaces, exhibit an +additional peak observed at a frequency of $\sim$170cm$^{-1}$, which +is attributed to the vibration of the S-Au bond. This vibration +enables efficient thermal transport from surface Au atoms to the +capping agents. Simultaneously, as shown in the lower panel of +Fig. \ref{vibration}, the large overlap of the vibration spectra of +butanethiol and hexane in the all-atom model, including the C-H +vibration, also suggests high thermal exchange efficiency. The +combination of these two effects produces the drastic interfacial +thermal conductance enhancement in the all-atom model. + +[MAY NEED TO CONVERT TO JPEG] \begin{figure} \includegraphics[width=\linewidth]{vibration} \caption{Vibrational spectra obtained for gold in different @@ -482,14 +825,48 @@ interfacial thermal conductance enhancement in the all all-atom model (lower panel).} \label{vibration} \end{figure} -% 600dpi, letter size. too large? +[COMPARISON OF TWO G'S; AU SLAB WIDTHS; ETC] +% The results show that the two definitions used for $G$ yield +% comparable values, though $G^\prime$ tends to be smaller. +\section{Conclusions} +The NIVS algorithm we developed has been applied to simulations of +Au-butanethiol surfaces with organic solvents. This algorithm allows +effective unphysical thermal flux transferred between the metal and +the liquid phase. With the flux applied, we were able to measure the +corresponding thermal gradient and to obtain interfacial thermal +conductivities. Our simulations have seen significant conductance +enhancement with the presence of capping agent, compared to the bare +gold/liquid interfaces. The acoustic impedance mismatch between the +metal and the liquid phase is effectively eliminated by proper capping +agent. Furthermore, the coverage precentage of the capping agent plays +an important role in the interfacial thermal transport process. + +Our measurement results, particularly of the UA models, agree with +available experimental data. This indicates that our force field +parameters have a nice description of the interactions between the +particles at the interfaces. AA models tend to overestimate the +interfacial thermal conductance in that the classically treated C-H +vibration would be overly sampled. Compared to the AA models, the UA +models have higher computational efficiency with satisfactory +accuracy, and thus are preferable in interfacial thermal transport +modelings. + +Vlugt {\it et al.} has investigated the surface thiol structures for +nanocrystal gold and pointed out that they differs from those of the +Au(111) surface\cite{vlugt:cpc2007154}. This difference might lead to +change of interfacial thermal transport behavior as well. To +investigate this problem, an effective means to introduce thermal flux +and measure the corresponding thermal gradient is desirable for +simulating structures with spherical symmetry. + + \section{Acknowledgments} Support for this project was provided by the National Science Foundation under grant CHE-0848243. Computational time was provided by the Center for Research Computing (CRC) at the University of Notre -Dame. \newpage +Dame. \newpage \bibliography{interfacial}