--- interfacial/interfacial.tex 2011/07/29 19:26:44 3756 +++ interfacial/interfacial.tex 2011/09/27 21:02:48 3763 @@ -44,24 +44,24 @@ Notre Dame, Indiana 46556} \begin{doublespace} \begin{abstract} - -With the Non-Isotropic Velocity Scaling algorithm (NIVS) we have -developed, an unphysical thermal flux can be effectively set up even -for non-homogeneous systems like interfaces in non-equilibrium -molecular dynamics simulations. In this work, this algorithm is -applied for simulating thermal conductance at metal / organic solvent -interfaces with various coverages of butanethiol capping -agents. Different solvents and force field models were tested. Our -results suggest that the United-Atom models are able to provide an -estimate of the interfacial thermal conductivity comparable to -experiments in our simulations with satisfactory computational -efficiency. From our results, the acoustic impedance mismatch between -metal and liquid phase is effectively reduced by the capping -agents, and thus leads to interfacial thermal conductance -enhancement. Furthermore, this effect is closely related to the -capping agent coverage on the metal surfaces and the type of solvent -molecules, and is affected by the models used in the simulations. - + With the Non-Isotropic Velocity Scaling (NIVS) approach to Reverse + Non-Equilibrium Molecular Dynamics (RNEMD) it is possible to impose + an unphysical thermal flux between different regions of + inhomogeneous systems such as solid / liquid interfaces. We have + applied NIVS to compute the interfacial thermal conductance at a + metal / organic solvent interface that has been chemically capped by + butanethiol molecules. Our calculations suggest that the acoustic + impedance mismatch between the metal and liquid phases is + effectively reduced by the capping agents, leading to a greatly + enhanced conductivity at the interface. Specifically, the chemical + bond between the metal and the capping agent introduces a + vibrational overlap that is not present without the capping agent, + and the overlap between the vibrational spectra (metal to cap, cap + to solvent) provides a mechanism for rapid thermal transport across + the interface. Our calculations also suggest that this is a + non-monotonic function of the fractional coverage of the surface, as + moderate coverages allow convective heat transport of solvent + molecules that have been in close contact with the capping agent. \end{abstract} \newpage @@ -73,28 +73,29 @@ Due to the importance of heat flow in nanotechnology, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} -Due to the importance of heat flow in nanotechnology, interfacial -thermal conductance has been studied extensively both experimentally -and computationally.\cite{cahill:793} Unlike bulk materials, nanoscale -materials have a significant fraction of their atoms at interfaces, -and the chemical details of these interfaces govern the heat transfer -behavior. Furthermore, the interfaces are +Due to the importance of heat flow (and heat removal) in +nanotechnology, interfacial thermal conductance has been studied +extensively both experimentally and computationally.\cite{cahill:793} +Nanoscale materials have a significant fraction of their atoms at +interfaces, and the chemical details of these interfaces govern the +thermal transport properties. Furthermore, the interfaces are often heterogeneous (e.g. solid - liquid), which provides a challenge to -traditional methods developed for homogeneous systems. +computational methods which have been developed for homogeneous or +bulk systems. -Experimentally, various interfaces have been investigated for their -thermal conductance. Cahill and coworkers studied nanoscale thermal +Experimentally, the thermal properties of a number of interfaces have +been investigated. Cahill and coworkers studied nanoscale thermal transport from metal nanoparticle/fluid interfaces, to epitaxial -TiN/single crystal oxides interfaces, to hydrophilic and hydrophobic +TiN/single crystal oxides interfaces, and 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 +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 +cetyltrimethylammonium bromide (CTAB) on the 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 resistance 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 @@ -114,14 +115,21 @@ methods\cite{MullerPlathe:1997xw,kuang:164101} would h 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. +methods\cite{MullerPlathe:1997xw,kuang:164101} would be advantageous +in that they {\it apply} the difficult to measure quantity (flux), +while {\it measuring} the easily-computed quantity (the thermal +gradient). This is particularly true for inhomogeneous interfaces +where it would not be clear how to apply a gradient {\it a priori}. 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. +conductance (or resistance) is dependent on chemical details of a +number of hydrophobic and hydrophilic aqueous interfaces. {\bf And + Luo {\it et al.} studied the thermal conductance of Au-SAM-Au + junctions using the same approach, with comparison to a constant + temperature difference method\cite{Luo20101}. While this latter + approach establishes more thermal distributions compared to the + former RNEMD methods, it does not guarantee momentum or kinetic + energy conservations.} Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS) algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm @@ -135,13 +143,15 @@ properties. Different models were used for both the ca 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. +properties. {\bf To our knowledge, few previous MD inverstigations + have been found to address to these systems yet.} Different models +were used for both the capping agent and the solvent force field +parameters. Using the NIVS algorithm, the thermal transport across +these interfaces was studied and the underlying mechanism for the +phenomena was investigated. \section{Methodology} -\subsection{Imposd-Flux Methods in MD Simulations} +\subsection{Imposed-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, @@ -165,12 +175,12 @@ can be applied between regions of particles of arbitar 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 +can be applied between regions of particles of arbitrary identity, and the flux will not be restricted by difference in particle mass. The NIVS algorithm scales the velocity vectors in two separate regions -of a simulation system with respective diagonal scaling matricies. To -determine these scaling factors in the matricies, a set of equations +of a simulation system with respective diagonal scaling matrices. To +determine these scaling factors in the matrices, a set of equations including linear momentum conservation and kinetic energy conservation constraints and target energy flux satisfaction is solved. With the scaling operation applied to the system in a set frequency, bulk @@ -202,7 +212,7 @@ occurs at the Gibbs deviding surface (Figure \ref{demo 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 +occurs at the Gibbs dividing surface (Figure \ref{demoPic}). This is known as the Kapitza conductance, which is the inverse of the Kapitza resistance. \begin{equation} @@ -215,12 +225,12 @@ resistance. \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.} + system responds by forming a thermal 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} @@ -259,10 +269,14 @@ profile. \begin{figure} \includegraphics[width=\linewidth]{gradT} -\caption{A sample of Au-butanethiol/hexane interfacial system and the - temperature profile after a kinetic energy flux is imposed to - it. The 1st and 2nd derivatives of the temperature profile can be - obtained with finite difference approximation (lower panel).} +\caption{A sample of Au (111) / butanethiol / hexane interfacial + system with the temperature profile after a kinetic energy flux has + been imposed. Note that the largest temperature jump in the thermal + profile (corresponding to the lowest interfacial conductance) is at + the interface between the butanethiol molecules (blue) and the + solvent (grey). First and second derivatives of the temperature + profile are obtained using a finite difference approximation (lower + panel).} \label{gradT} \end{figure} @@ -309,7 +323,8 @@ between periodic images of the gold interfaces is $45 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. +between periodic images of the gold interfaces is $45 \sim 75$\AA in +our simulations. The initial configurations generated are further equilibrated with the $x$ and $y$ dimensions fixed, only allowing the $z$-length scale to @@ -327,8 +342,8 @@ gradient had stablized, the temperature profile of the $\sim$200K. Therefore, thermal flux usually came from the metal to the liquid so that the liquid has a higher temperature and would not freeze due to lowered temperatures. After this induced temperature -gradient had stablized, the temperature profile of the simulation cell -was recorded. To do this, the simulation cell is devided evenly into +gradient had stabilized, the temperature profile of the simulation cell +was recorded. To do this, the simulation cell is divided 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 @@ -359,8 +374,8 @@ particles of different species. \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} + (AA) force fields. Most parameters are from References + \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} @@ -384,7 +399,7 @@ However, the TraPPE-UA model for alkanes is known to p 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 +However, the TraPPE-UA model for alkanes is known to predict a slightly 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 @@ -503,9 +518,10 @@ 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.} +\caption{The interfacial thermal conductivity ($G$) has a + non-monotonic dependence on the degree of surface capping. This + data is for the Au(111) / butanethiol / solvent interface with + various UA force fields at $\langle T\rangle \sim $200K.} \label{coverage} \end{figure} @@ -519,7 +535,7 @@ an enhancement by at least a factor of 3. Cappping ag 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 +an enhancement by at least a factor of 3. Capping agents are clearly playing a major role in thermal transport at metal / organic solvent surfaces. @@ -577,10 +593,11 @@ studies. \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.)} + $\langle T\rangle\sim$200K. Here ``D'' stands for deuterated + solvent or capping agent molecules; ``Avg.'' denotes results + that are averages of simulations under different applied + thermal flux $(J_z)$ values. Error estimates are indicated in + parentheses.} \begin{tabular}{llccc} \hline\hline @@ -627,7 +644,7 @@ capping agent and solvent has a large impact on the ca phases. For the fully-covered surfaces, the choice of force field for the -capping agent and solvent has a large impact on the calulated values +capping agent and solvent has a large impact on the calculated 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 @@ -682,7 +699,7 @@ arbitary if one aims to obtain a stable and reliable t flow simply by changing the sign of the flux, and thermal gradients which would be difficult to obtain experimentally ($5$ K/\AA) can be easily simulated. However, the magnitude of the applied flux is not -arbitary if one aims to obtain a stable and reliable thermal gradient. +arbitrary if one aims to obtain a stable and reliable thermal gradient. A temperature gradient can be lost in the noise if $|J_z|$ is too small, and excessive $|J_z|$ values can cause phase transitions if the extremes of the simulation cell become widely separated in @@ -707,11 +724,12 @@ small number $J_z$ values. \begin{table*} \begin{minipage}{\linewidth} \begin{center} - \caption{Computed interfacial thermal conductivity ($G$ and - $G^\prime$) values for the 100\% covered Au-butanethiol/hexane - interfaces with UA model and different hexane molecule numbers - at different temperatures using a range of energy - fluxes. Error estimates indicated in parenthesis.} + \caption{In the hexane-solvated interfaces, the system size has + little effect on the calculated values for interfacial + conductance ($G$ and $G^\prime$), but the direction of heat + flow (i.e. the sign of $J_z$) can alter the average + temperature of the liquid phase and this can alter the + computed conductivity.} \begin{tabular}{ccccccc} \hline\hline @@ -760,14 +778,14 @@ Table \ref{AuThiolHexaneUA}. In raising the average t between the hexane and butanethiol, and this accounts for our observation of lower conductance at higher temperatures as shown in Table \ref{AuThiolHexaneUA}. In raising the average temperature from -200K to 250K, the density drop of ~20\% in the solvent phase leads to -a ~65\% drop in the conductance. +200K to 250K, the density drop of $\sim$20\% in the solvent phase +leads to a $\sim$40\% drop in the conductance. Similar behavior is observed in the TraPPE-UA model for toluene, although this model has better agreement with the experimental densities of toluene. The expansion of the toluene liquid phase is not as significant as that of the hexane (8.3\% over 100K), and this -limits the effect to ~20\% drop in thermal conductivity (Table +limits the effect to $\sim$20\% drop in thermal conductivity (Table \ref{AuThiolToluene}). Although we have not mapped out the behavior at a large number of @@ -779,10 +797,10 @@ function of temperature. \begin{table*} \begin{minipage}{\linewidth} \begin{center} - \caption{Computed interfacial thermal conductivity ($G$ and - $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene - interface at different temperatures using a range of energy - fluxes. Error estimates indicated in parenthesis.} + \caption{When toluene is the solvent, the interfacial thermal + conductivity is less sensitive to temperature, but again, the + direction of the heat flow can alter the solvent temperature + and can change the computed conductance values.} \begin{tabular}{ccccc} \hline\hline @@ -809,7 +827,7 @@ potential. This phenomenon agrees with reconstructions S atom layer, and butanethiol molecules which embed into the surface Au layer. The driving force for this behavior is the strong Au-S interactions which are modeled here with a deep Lennard-Jones -potential. This phenomenon agrees with reconstructions that have beeen +potential. This phenomenon agrees with reconstructions that have been experimentally observed.\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. Vlugt {\it et al.} kept their Au(111) slab rigid so that their simulations @@ -822,11 +840,11 @@ reconstruction. O ur Au / butanethiol / toluene system However, when the surface is not completely covered by butanethiols, the simulated system appears to be more resistent to the -reconstruction. O ur Au / butanethiol / toluene system had the Au(111) +reconstruction. 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. That said, we did observe butanethiols migrating to neighboring three-fold sites during -a simulation. Since the interface persisted in these simulations, +a simulation. Since the interface persisted in these simulations, we were able to obtain $G$'s for these interfaces even at a relatively high temperature without being affected by surface reconstructions. @@ -845,7 +863,7 @@ that is higher than that of the fastest vibrations occ spectra, simulations were run after equilibration in the microcanonical (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 +that is higher than that of the fastest vibrations occurring in the simulations. With these configurations, the velocity auto-correlation functions can be computed: \begin{equation} @@ -865,24 +883,29 @@ $\sim$170cm$^{-1}$. We attribute this peak to the S-A shown as in Figure \ref{specAu}. Regardless of the presence of solvent, the gold surfaces which are covered by butanethiol molecules exhibit an additional peak observed at a frequency of -$\sim$170cm$^{-1}$. We attribute this peak to the S-Au bonding +$\sim$165cm$^{-1}$. We attribute this peak to the S-Au bonding vibration. This vibration enables efficient thermal coupling of the surface Au layer to the capping agents. Therefore, in our simulations, the Au / S interfaces do not appear to be the primary barrier to thermal transport when compared with the butanethiol / solvent -interfaces. +interfaces. {\bf This confirms the results from Luo {\it et + al.}\cite{Luo20101}, which reported $G$ for Au-SAM junctions + generally twice larger than what we have computed for the + thiol-liquid interfaces.} \begin{figure} \includegraphics[width=\linewidth]{vibration} -\caption{Vibrational power spectra for gold in different solvent - environments. The presence of the butanethiol capping molecules - adds a vibrational peak at $\sim$170cm$^{-1}$.} +\caption{The vibrational power spectrum for thiol-capped gold has an + additional vibrational peak at $\sim $165cm$^{-1}$. Bare gold + surfaces (both with and without a solvent over-layer) are missing + this peak. A similar peak at $\sim $165cm$^{-1}$ also appears in + the vibrational power spectrum for the butanethiol capping agents.} \label{specAu} \end{figure} Also in this figure, we show the vibrational power spectrum for the bound butanethiol molecules, which also exhibits the same -$\sim$170cm$^{-1}$ peak. +$\sim$165cm$^{-1}$ peak. \subsection{Overlap of power spectra} A comparison of the results obtained from the two different organic @@ -894,7 +917,11 @@ in the UA force field. vibrations were treated classically. The presence of extra degrees of freedom in the AA force field yields higher heat exchange rates between the two phases and results in a much higher conductivity than -in the UA force field. +in the UA force field. {\bf Due to the classical models used, this + even includes those high frequency modes which should be unpopulated + at our relatively low temperatures. This artifact causes high + frequency vibrations accountable for thermal transport in classical + MD simulations.} The similarity in the vibrational modes available to solvent and capping agent can be reduced by deuterating one of the two components @@ -939,10 +966,12 @@ produce an observable difference for the results of $G \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.} +\caption{Vibrational power spectra for UA models for the butanethiol + and hexane solvent (upper panel) show the high degree of overlap + between these two molecules, particularly at lower frequencies. + Deuterating a UA model for the solvent (lower panel) does not + decouple the two spectra to the same degree as in the AA force + field (see Fig \ref{aahxntln}).} \label{uahxnua} \end{figure} @@ -960,7 +989,7 @@ agent. Furthermore, the coverage precentage of the cap presence of capping agent, compared with the bare gold / liquid interfaces. The acoustic impedance mismatch between the metal and the liquid phase is effectively eliminated by a chemically-bonded capping -agent. Furthermore, the coverage precentage of the capping agent plays +agent. Furthermore, the coverage percentage of the capping agent plays an important role in the interfacial thermal transport process. Moderately low coverages allow higher contact between capping agent and solvent, and thus could further enhance the heat transfer @@ -997,6 +1026,11 @@ Dame. Foundation under grant CHE-0848243. Computational time was provided by the Center for Research Computing (CRC) at the University of Notre Dame. + +\section{Supporting Information} +This information is available free of charge via the Internet at +http://pubs.acs.org. + \newpage \bibliography{interfacial}