ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/interfacial/interfacial.tex
(Generate patch)

Comparing interfacial/interfacial.tex (file contents):
Revision 3749 by skuang, Mon Jul 25 19:11:33 2011 UTC vs.
Revision 3755 by skuang, Fri Jul 29 15:45:14 2011 UTC

# Line 73 | Line 73 | Interfacial thermal conductance is extensively studied
73   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74  
75   \section{Introduction}
76 < Interfacial thermal conductance is extensively studied both
77 < experimentally and computationally\cite{cahill:793}, due to its
78 < importance in nanoscale science and technology. Reliability of
79 < nanoscale devices depends on their thermal transport
80 < properties. Unlike bulk homogeneous materials, nanoscale materials
81 < features significant presence of interfaces, and these interfaces
82 < could dominate the heat transfer behavior of these
83 < materials. Furthermore, these materials are generally heterogeneous,
84 < which challenges traditional research methods for homogeneous
85 < systems.
76 > Due to the importance of heat flow in nanotechnology, interfacial
77 > thermal conductance has been studied extensively both experimentally
78 > and computationally.\cite{cahill:793} Unlike bulk materials, nanoscale
79 > materials have a significant fraction of their atoms at interfaces,
80 > and the chemical details of these interfaces govern the heat transfer
81 > behavior. Furthermore, the interfaces are
82 > heterogeneous (e.g. solid - liquid), which provides a challenge to
83 > traditional methods developed for homogeneous systems.
84  
85 < Heat conductance of molecular and nano-scale interfaces will be
86 < affected by the chemical details of the surface. Experimentally,
87 < various interfaces have been investigated for their thermal
88 < conductance properties. Wang {\it et al.} studied heat transport
89 < through long-chain hydrocarbon monolayers on gold substrate at
90 < individual molecular level\cite{Wang10082007}; Schmidt {\it et al.}
91 < studied the role of CTAB on thermal transport between gold nanorods
92 < and solvent\cite{doi:10.1021/jp8051888}; Juv\'e {\it et al.} studied
85 > Experimentally, various interfaces have been investigated for their
86 > thermal conductance. Cahill and coworkers studied nanoscale thermal
87 > transport from metal nanoparticle/fluid interfaces, to epitaxial
88 > TiN/single crystal oxides interfaces, to hydrophilic and hydrophobic
89 > interfaces between water and solids with different self-assembled
90 > monolayers.\cite{Wilson:2002uq,PhysRevB.67.054302,doi:10.1021/jp048375k,PhysRevLett.96.186101}
91 > Wang {\it et al.} studied heat transport through
92 > long-chain hydrocarbon monolayers on gold substrate at individual
93 > molecular level,\cite{Wang10082007} Schmidt {\it et al.} studied the
94 > role of CTAB on thermal transport between gold nanorods and
95 > solvent,\cite{doi:10.1021/jp8051888} and Juv\'e {\it et al.} studied
96   the cooling dynamics, which is controlled by thermal interface
97   resistence of glass-embedded metal
98 < nanoparticles\cite{PhysRevB.80.195406}. Although interfaces are
99 < commonly barriers for heat transport, Alper {\it et al.} suggested
100 < that specific ligands (capping agents) could completely eliminate this
101 < barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}.
98 > nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
99 > normally considered barriers for heat transport, Alper {\it et al.}
100 > suggested that specific ligands (capping agents) could completely
101 > eliminate this barrier
102 > ($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}
103  
104   Theoretical and computational models have also been used to study the
105   interfacial thermal transport in order to gain an understanding of
# Line 105 | Line 107 | atoms)\cite{hase:2010,hase:2011}. However, ensemble av
107   employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
108   study thermal transport from hot Au(111) substrate to a self-assembled
109   monolayer of alkylthiol with relatively long chain (8-20 carbon
110 < atoms)\cite{hase:2010,hase:2011}. However, ensemble averaged
110 > atoms).\cite{hase:2010,hase:2011} However, ensemble averaged
111   measurements for heat conductance of interfaces between the capping
112 < monolayer on Au and a solvent phase has yet to be studied.
113 < The comparatively low thermal flux through interfaces is
114 < difficult to measure with Equilibrium MD or forward NEMD simulation
115 < methods. Therefore, the Reverse NEMD (RNEMD) methods would have the
116 < advantage of having this difficult to measure flux known when studying
117 < the thermal transport across interfaces, given that the simulation
118 < methods being able to effectively apply an unphysical flux in
119 < non-homogeneous systems.
112 > monolayer on Au and a solvent phase have yet to be studied with their
113 > approach. The comparatively low thermal flux through interfaces is
114 > difficult to measure with Equilibrium
115 > MD\cite{doi:10.1080/0026897031000068578} or forward NEMD simulation
116 > methods. Therefore, the Reverse NEMD (RNEMD)
117 > methods\cite{MullerPlathe:1997xw,kuang:164101} would have the
118 > advantage of applying this difficult to measure flux (while measuring
119 > the resulting gradient), given that the simulation methods being able
120 > to effectively apply an unphysical flux in non-homogeneous systems.
121 > Garde and coworkers\cite{garde:nl2005,garde:PhysRevLett2009} applied
122 > this approach to various liquid interfaces and studied how thermal
123 > conductance (or resistance) is dependent on chemistry details of
124 > interfaces, e.g. hydrophobic and hydrophilic aqueous interfaces.
125  
126 < Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
126 > Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS)
127   algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
128   retains the desirable features of RNEMD (conservation of linear
129   momentum and total energy, compatibility with periodic boundary
# Line 133 | Line 140 | underlying mechanism for the phenomena was investigate
140   thermal transport across these interfaces was studied and the
141   underlying mechanism for the phenomena was investigated.
142  
136 [MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]
137
143   \section{Methodology}
144   \subsection{Imposd-Flux Methods in MD Simulations}
145 < Steady state MD simulations has the advantage that not many
145 > Steady state MD simulations have an advantage in that not many
146   trajectories are needed to study the relationship between thermal flux
147 < and thermal gradients. For systems including low conductance
148 < interfaces one must have a method capable of generating or measuring
149 < relatively small fluxes, compared to those required for bulk
150 < conductivity. This requirement makes the calculation even more
151 < difficult for those slowly-converging equilibrium
152 < methods\cite{Viscardy:2007lq}. Forward methods may impose gradient,
153 < but in interfacial conditions it is not clear what behavior to impose
154 < at the interfacial boundaries. Imposed-flux reverse non-equilibrium
155 < methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
156 < the thermal response becomes easier to measure than the flux. Although
147 > and thermal gradients. For systems with low interfacial conductance,
148 > one must have a method capable of generating or measuring relatively
149 > small fluxes, compared to those required for bulk conductivity. This
150 > requirement makes the calculation even more difficult for
151 > slowly-converging equilibrium methods.\cite{Viscardy:2007lq} Forward
152 > NEMD methods impose a gradient (and measure a flux), but at interfaces
153 > it is not clear what behavior should be imposed at the boundaries
154 > between materials.  Imposed-flux reverse non-equilibrium
155 > methods\cite{MullerPlathe:1997xw} set the flux {\it a priori} and
156 > the thermal response becomes an easy-to-measure quantity.  Although
157   M\"{u}ller-Plathe's original momentum swapping approach can be used
158   for exchanging energy between particles of different identity, the
159   kinetic energy transfer efficiency is affected by the mass difference
160   between the particles, which limits its application on heterogeneous
161   interfacial systems.
162  
163 < The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach to
164 < non-equilibrium MD simulations is able to impose a wide range of
163 > The non-isotropic velocity scaling (NIVS) \cite{kuang:164101} approach
164 > to non-equilibrium MD simulations is able to impose a wide range of
165   kinetic energy fluxes without obvious perturbation to the velocity
166   distributions of the simulated systems. Furthermore, this approach has
167   the advantage in heterogeneous interfaces in that kinetic energy flux
# Line 173 | Line 178 | momenta and energy and does not depend on an external
178   for computing thermal conductivities. The NIVS algorithm conserves
179   momenta and energy and does not depend on an external thermostat.
180  
181 < \subsection{Defining Interfacial Thermal Conductivity $G$}
182 < Given a system with thermal gradients and the corresponding thermal
183 < flux, for interfaces with a relatively low interfacial conductance,
184 < the bulk regions on either side of an interface rapidly come to a
185 < state in which the two phases have relatively homogeneous (but
186 < distinct) temperatures. The interfacial thermal conductivity $G$ can
187 < therefore be approximated as:
181 > \subsection{Defining Interfacial Thermal Conductivity ($G$)}
182 >
183 > For an interface with relatively low interfacial conductance, and a
184 > thermal flux between two distinct bulk regions, the regions on either
185 > side of the interface rapidly come to a state in which the two phases
186 > have relatively homogeneous (but distinct) temperatures. The
187 > interfacial thermal conductivity $G$ can therefore be approximated as:
188   \begin{equation}
189 < G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
189 >  G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
190      \langle T_\mathrm{cold}\rangle \right)}
191   \label{lowG}
192   \end{equation}
193 < where ${E_{total}}$ is the imposed non-physical kinetic energy
194 < transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle
195 <  T_\mathrm{cold}\rangle}$ are the average observed temperature of the
196 < two separated phases.
193 > where ${E_{total}}$ is the total imposed non-physical kinetic energy
194 > transfer during the simulation and ${\langle T_\mathrm{hot}\rangle}$
195 > and ${\langle T_\mathrm{cold}\rangle}$ are the average observed
196 > temperature of the two separated phases.
197  
198   When the interfacial conductance is {\it not} small, there are two
199 < ways to define $G$.
200 <
201 < One way is to assume the temperature is discrete on the two sides of
202 < the interface. $G$ can be calculated using the applied thermal flux
203 < $J$ and the maximum temperature difference measured along the thermal
204 < gradient max($\Delta T$), which occurs at the Gibbs deviding surface
205 < (Figure \ref{demoPic}):
199 > ways to define $G$. One common way is to assume the temperature is
200 > discrete on the two sides of the interface. $G$ can be calculated
201 > using the applied thermal flux $J$ and the maximum temperature
202 > difference measured along the thermal gradient max($\Delta T$), which
203 > occurs at the Gibbs deviding surface (Figure \ref{demoPic}). This is
204 > known as the Kapitza conductance, which is the inverse of the Kapitza
205 > resistance.
206   \begin{equation}
207 < G=\frac{J}{\Delta T}
207 >  G=\frac{J}{\Delta T}
208   \label{discreteG}
209   \end{equation}
210  
# Line 219 | Line 224 | the magnitude of thermal conductivity $\lambda$ change
224  
225   The other approach is to assume a continuous temperature profile along
226   the thermal gradient axis (e.g. $z$) and define $G$ at the point where
227 < the magnitude of thermal conductivity $\lambda$ change reach its
227 > the magnitude of thermal conductivity ($\lambda$) change reaches its
228   maximum, given that $\lambda$ is well-defined throughout the space:
229   \begin{equation}
230   G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
# Line 230 | Line 235 | With the temperature profile obtained from simulations
235   \label{derivativeG}
236   \end{equation}
237  
238 < With the temperature profile obtained from simulations, one is able to
238 > With temperature profiles obtained from simulation, one is able to
239   approximate the first and second derivatives of $T$ with finite
240 < difference methods and thus calculate $G^\prime$.
240 > difference methods and calculate $G^\prime$. In what follows, both
241 > definitions have been used, and are compared in the results.
242  
243 < In what follows, both definitions have been used for calculation and
244 < are compared in the results.
245 <
246 < To compare the above definitions ($G$ and $G^\prime$), we have modeled
247 < a metal slab with its (111) surfaces perpendicular to the $z$-axis of
242 < our simulation cells. Both with and without capping agents on the
243 < surfaces, the metal slab is solvated with simple organic solvents, as
243 > To investigate the interfacial conductivity at metal / solvent
244 > interfaces, we have modeled a metal slab with its (111) surfaces
245 > perpendicular to the $z$-axis of our simulation cells. The metal slab
246 > has been prepared both with and without capping agents on the exposed
247 > surface, and has been solvated with simple organic solvents, as
248   illustrated in Figure \ref{gradT}.
249  
250   With the simulation cell described above, we are able to equilibrate
251   the system and impose an unphysical thermal flux between the liquid
252   and the metal phase using the NIVS algorithm. By periodically applying
253 < the unphysical flux, we are able to obtain a temperature profile and
254 < its spatial derivatives. These quantities enable the evaluation of the
255 < interfacial thermal conductance of a surface. Figure \ref{gradT} is an
256 < example of how an applied thermal flux can be used to obtain the 1st
253 < and 2nd derivatives of the temperature profile.
253 > the unphysical flux, we obtained a temperature profile and its spatial
254 > derivatives. Figure \ref{gradT} shows how an applied thermal flux can
255 > be used to obtain the 1st and 2nd derivatives of the temperature
256 > profile.
257  
258   \begin{figure}
259   \includegraphics[width=\linewidth]{gradT}
# Line 264 | Line 267 | OpenMD\cite{Meineke:2005gd,openmd}, and was used for o
267   \section{Computational Details}
268   \subsection{Simulation Protocol}
269   The NIVS algorithm has been implemented in our MD simulation code,
270 < OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
271 < simulations. Different metal slab thickness (layer numbers of Au) was
272 < simulated. Metal slabs were first equilibrated under atmospheric
273 < pressure (1 atm) and a desired temperature (e.g. 200K). After
274 < equilibration, butanethiol capping agents were placed at three-fold
275 < hollow sites on the Au(111) surfaces. These sites could be either a
276 < {\it fcc} or {\it hcp} site. However, Hase {\it et al.} found that
277 < they are equivalent in a heat transfer process\cite{hase:2010}, so
275 < they are not distinguished in our study. The maximum butanethiol
270 > OpenMD\cite{Meineke:2005gd,openmd}, and was used for our simulations.
271 > Metal slabs of 6 or 11 layers of Au atoms were first equilibrated
272 > under atmospheric pressure (1 atm) and 200K. After equilibration,
273 > butanethiol capping agents were placed at three-fold hollow sites on
274 > the Au(111) surfaces. These sites are either {\it fcc} or {\it
275 >  hcp} sites, although Hase {\it et al.} found that they are
276 > equivalent in a heat transfer process,\cite{hase:2010} so we did not
277 > distinguish between these sites in our study. The maximum butanethiol
278   capacity on Au surface is $1/3$ of the total number of surface Au
279   atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
280   structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
281 < series of different coverages was derived by evenly eliminating
282 < butanethiols on the surfaces, and was investigated in order to study
283 < the relation between coverage and interfacial conductance.
281 > series of lower coverages was also prepared by eliminating
282 > butanethiols from the higher coverage surface in a regular manner. The
283 > lower coverages were prepared in order to study the relation between
284 > coverage and interfacial conductance.
285  
286   The capping agent molecules were allowed to migrate during the
287   simulations. They distributed themselves uniformly and sampled a
288   number of three-fold sites throughout out study. Therefore, the
289 < initial configuration would not noticeably affect the sampling of a
289 > initial configuration does not noticeably affect the sampling of a
290   variety of configurations of the same coverage, and the final
291   conductance measurement would be an average effect of these
292 < configurations explored in the simulations. [MAY NEED SNAPSHOTS]
292 > configurations explored in the simulations.
293  
294 < After the modified Au-butanethiol surface systems were equilibrated
295 < under canonical ensemble, organic solvent molecules were packed in the
296 < previously empty part of the simulation cells\cite{packmol}. Two
294 > After the modified Au-butanethiol surface systems were equilibrated in
295 > the canonical (NVT) ensemble, organic solvent molecules were packed in
296 > the previously empty part of the simulation cells.\cite{packmol} Two
297   solvents were investigated, one which has little vibrational overlap
298 < with the alkanethiol and a planar shape (toluene), and one which has
299 < similar vibrational frequencies and chain-like shape ({\it n}-hexane).
298 > with the alkanethiol and which has a planar shape (toluene), and one
299 > which has similar vibrational frequencies to the capping agent and
300 > chain-like shape ({\it n}-hexane).
301  
302 < The space filled by solvent molecules, i.e. the gap between
303 < periodically repeated Au-butanethiol surfaces should be carefully
304 < chosen. A very long length scale for the thermal gradient axis ($z$)
301 < may cause excessively hot or cold temperatures in the middle of the
302 > The simulation cells were not particularly extensive along the
303 > $z$-axis, as a very long length scale for the thermal gradient may
304 > cause excessively hot or cold temperatures in the middle of the
305   solvent region and lead to undesired phenomena such as solvent boiling
306   or freezing when a thermal flux is applied. Conversely, too few
307   solvent molecules would change the normal behavior of the liquid
308   phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
309 < these extreme cases did not happen to our simulations. And the
310 < corresponding spacing is usually $35 \sim 75$\AA.
309 > these extreme cases did not happen to our simulations. The spacing
310 > between periodic images of the gold interfaces is $45 \sim 75$\AA.
311  
312   The initial configurations generated are further equilibrated with the
313 < $x$ and $y$ dimensions fixed, only allowing length scale change in $z$
314 < dimension. This is to ensure that the equilibration of liquid phase
315 < does not affect the metal crystal structure in $x$ and $y$ dimensions.
316 < To investigate this effect, comparisons were made with simulations
317 < that allow changes of $L_x$ and $L_y$ during NPT equilibration, and
318 < the results are shown in later sections. After ensuring the liquid
319 < phase reaches equilibrium at atmospheric pressure (1 atm), further
320 < equilibration are followed under NVT and then NVE ensembles.
313 > $x$ and $y$ dimensions fixed, only allowing the $z$-length scale to
314 > change. This is to ensure that the equilibration of liquid phase does
315 > not affect the metal's crystalline structure. Comparisons were made
316 > with simulations that allowed changes of $L_x$ and $L_y$ during NPT
317 > equilibration. No substantial changes in the box geometry were noticed
318 > in these simulations. After ensuring the liquid phase reaches
319 > equilibrium at atmospheric pressure (1 atm), further equilibration was
320 > carried out under canonical (NVT) and microcanonical (NVE) ensembles.
321  
322 < After the systems reach equilibrium, NIVS is implemented to impose a
323 < periodic unphysical thermal flux between the metal and the liquid
324 < phase. Most of our simulations are under an average temperature of
325 < $\sim$200K. Therefore, this flux usually comes from the metal to the
322 > After the systems reach equilibrium, NIVS was used to impose an
323 > unphysical thermal flux between the metal and the liquid phases. Most
324 > of our simulations were done under an average temperature of
325 > $\sim$200K. Therefore, thermal flux usually came from the metal to the
326   liquid so that the liquid has a higher temperature and would not
327 < freeze due to excessively low temperature. After this induced
328 < temperature gradient is stablized, the temperature profile of the
329 < simulation cell is recorded. To do this, the simulation cell is
330 < devided evenly into $N$ slabs along the $z$-axis and $N$ is maximized
328 < for highest possible spatial resolution but not too many to have some
329 < slabs empty most of the time. The average temperatures of each slab
327 > freeze due to lowered temperatures. After this induced temperature
328 > gradient had stablized, the temperature profile of the simulation cell
329 > was recorded. To do this, the simulation cell is devided evenly into
330 > $N$ slabs along the $z$-axis. The average temperatures of each slab
331   are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
332   the same, the derivatives of $T$ with respect to slab number $n$ can
333 < be directly used for $G^\prime$ calculations:
334 < \begin{equation}
334 < G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
333 > be directly used for $G^\prime$ calculations: \begin{equation}
334 >  G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
335           \Big/\left(\frac{\partial T}{\partial z}\right)^2
336           = |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
337           \Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
# Line 340 | Line 340 | All of the above simulation procedures use a time step
340   \label{derivativeG2}
341   \end{equation}
342  
343 < All of the above simulation procedures use a time step of 1 fs. And
344 < each equilibration / stabilization step usually takes 100 ps, or
345 < longer, if necessary.
343 > All of the above simulation procedures use a time step of 1 fs. Each
344 > equilibration stage took a minimum of 100 ps, although in some cases,
345 > longer equilibration stages were utilized.
346  
347   \subsection{Force Field Parameters}
348 < Our simulations include various components. Figure \ref{demoMol}
349 < demonstrates the sites defined for both United-Atom and All-Atom
350 < models of the organic solvent and capping agent molecules in our
351 < simulations. Force field parameter descriptions are needed for
348 > Our simulations include a number of chemically distinct components.
349 > Figure \ref{demoMol} demonstrates the sites defined for both
350 > United-Atom and All-Atom models of the organic solvent and capping
351 > agents in our simulations. Force field parameters are needed for
352   interactions both between the same type of particles and between
353   particles of different species.
354  
# Line 358 | Line 358 | particles of different species.
358    these simulations. The chemically-distinct sites (a-e) are expanded
359    in terms of constituent atoms for both United Atom (UA) and All Atom
360    (AA) force fields.  Most parameters are from
361 <  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} (UA) and
362 <  \protect\cite{OPLSAA} (AA).  Cross-interactions with the Au atoms are given
363 <  in Table \ref{MnM}.}
361 >  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes,TraPPE-UA.thiols}
362 >  (UA) and \protect\cite{OPLSAA} (AA). Cross-interactions with the Au
363 >  atoms are given in Table \ref{MnM}.}
364   \label{demoMol}
365   \end{figure}
366  
367   The Au-Au interactions in metal lattice slab is described by the
368 < quantum Sutton-Chen (QSC) formulation\cite{PhysRevB.59.3527}. The QSC
368 > quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
369   potentials include zero-point quantum corrections and are
370   reparametrized for accurate surface energies compared to the
371 < Sutton-Chen potentials\cite{Chen90}.
371 > Sutton-Chen potentials.\cite{Chen90}
372  
373 < For both solvent molecules, straight chain {\it n}-hexane and aromatic
374 < toluene, United-Atom (UA) and All-Atom (AA) models are used
375 < respectively. The TraPPE-UA
373 > For the two solvent molecules, {\it n}-hexane and toluene, two
374 > different atomistic models were utilized. Both solvents were modeled
375 > using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
376   parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
377   for our UA solvent molecules. In these models, sites are located at
378   the carbon centers for alkyl groups. Bonding interactions, including
379   bond stretches and bends and torsions, were used for intra-molecular
380 < sites not separated by more than 3 bonds. Otherwise, for non-bonded
381 < interactions, Lennard-Jones potentials are used. [CHECK CITATION]
380 > sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones
381 > potentials are used.
382  
383 < By eliminating explicit hydrogen atoms, these models are simple and
384 < computationally efficient, while maintains good accuracy. However, the
385 < TraPPE-UA for alkanes is known to predict a lower boiling point than
386 < experimental values. Considering that after an unphysical thermal flux
387 < is applied to a system, the temperature of ``hot'' area in the liquid
388 < phase would be significantly higher than the average of the system, to
389 < prevent over heating and boiling of the liquid phase, the average
390 < temperature in our simulations should be much lower than the liquid
391 < boiling point.
383 > By eliminating explicit hydrogen atoms, the TraPPE-UA models are
384 > simple and computationally efficient, while maintaining good accuracy.
385 > However, the TraPPE-UA model for alkanes is known to predict a slighly
386 > lower boiling point than experimental values. This is one of the
387 > reasons we used a lower average temperature (200K) for our
388 > simulations. If heat is transferred to the liquid phase during the
389 > NIVS simulation, the liquid in the hot slab can actually be
390 > substantially warmer than the mean temperature in the simulation. The
391 > lower mean temperatures therefore prevent solvent boiling.
392  
393 < For UA-toluene model, the non-bonded potentials between
394 < inter-molecular sites have a similar Lennard-Jones formulation. For
395 < intra-molecular interactions, considering the stiffness of the benzene
396 < ring, rigid body constraints are applied for further computational
397 < efficiency. All bonds in the benzene ring and between the ring and the
398 < methyl group remain rigid during the progress of simulations.
393 > For UA-toluene, the non-bonded potentials between intermolecular sites
394 > have a similar Lennard-Jones formulation. The toluene molecules were
395 > treated as a single rigid body, so there was no need for
396 > intramolecular interactions (including bonds, bends, or torsions) in
397 > this solvent model.
398  
399   Besides the TraPPE-UA models, AA models for both organic solvents are
400 < included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
401 < force field is used. Additional explicit hydrogen sites were
400 > included in our studies as well. The OPLS-AA\cite{OPLSAA} force fields
401 > were used. For hexane, additional explicit hydrogen sites were
402   included. Besides bonding and non-bonded site-site interactions,
403   partial charges and the electrostatic interactions were added to each
404 < CT and HC site. For toluene, the United Force Field developed by
405 < Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} is
406 < adopted. Without the rigid body constraints, bonding interactions were
408 < included. For the aromatic ring, improper torsions (inversions) were
409 < added as an extra potential for maintaining the planar shape.
410 < [CHECK CITATION]
404 > CT and HC site. For toluene, a flexible model for the toluene molecule
405 > was utilized which included bond, bend, torsion, and inversion
406 > potentials to enforce ring planarity.
407  
408 < The capping agent in our simulations, the butanethiol molecules can
409 < either use UA or AA model. The TraPPE-UA force fields includes
408 > The butanethiol capping agent in our simulations, were also modeled
409 > with both UA and AA model. The TraPPE-UA force field includes
410   parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
411   UA butanethiol model in our simulations. The OPLS-AA also provides
412   parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
413 < surfaces do not have the hydrogen atom bonded to sulfur. To adapt this
414 < change and derive suitable parameters for butanethiol adsorbed on
415 < Au(111) surfaces, we adopt the S parameters from Luedtke and
416 < Landman\cite{landman:1998}[CHECK CITATION]
417 < and modify parameters for its neighbor C
418 < atom for charge balance in the molecule. Note that the model choice
419 < (UA or AA) of capping agent can be different from the
420 < solvent. Regardless of model choice, the force field parameters for
425 < interactions between capping agent and solvent can be derived using
426 < Lorentz-Berthelot Mixing Rule:
413 > surfaces do not have the hydrogen atom bonded to sulfur. To derive
414 > suitable parameters for butanethiol adsorbed on Au(111) surfaces, we
415 > adopt the S parameters from Luedtke and Landman\cite{landman:1998} and
416 > modify the parameters for the CTS atom to maintain charge neutrality
417 > in the molecule.  Note that the model choice (UA or AA) for the capping
418 > agent can be different from the solvent. Regardless of model choice,
419 > the force field parameters for interactions between capping agent and
420 > solvent can be derived using Lorentz-Berthelot Mixing Rule:
421   \begin{eqnarray}
422 < \sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
423 < \epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
422 >  \sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
423 >  \epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
424   \end{eqnarray}
425  
426 < To describe the interactions between metal Au and non-metal capping
427 < agent and solvent particles, we refer to an adsorption study of alkyl
428 < thiols on gold surfaces by Vlugt {\it et
429 <  al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones
430 < form of potential parameters for the interaction between Au and
431 < pseudo-atoms CH$_x$ and S based on a well-established and widely-used
432 < effective potential of Hautman and Klein\cite{hautman:4994} for the
433 < Au(111) surface. As our simulations require the gold lattice slab to
434 < be non-rigid so that it could accommodate kinetic energy for thermal
441 < transport study purpose, the pair-wise form of potentials is
442 < preferred.
426 > To describe the interactions between metal (Au) and non-metal atoms,
427 > we refer to an adsorption study of alkyl thiols on gold surfaces by
428 > Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
429 > Lennard-Jones form of potential parameters for the interaction between
430 > Au and pseudo-atoms CH$_x$ and S based on a well-established and
431 > widely-used effective potential of Hautman and Klein for the Au(111)
432 > surface.\cite{hautman:4994} As our simulations require the gold slab
433 > to be flexible to accommodate thermal excitation, the pair-wise form
434 > of potentials they developed was used for our study.
435  
436 < Besides, the potentials developed from {\it ab initio} calculations by
437 < Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
438 < interactions between Au and aromatic C/H atoms in toluene. A set of
439 < pseudo Lennard-Jones parameters were provided for Au in their force
440 < fields. By using the Mixing Rule, this can be used to derive pair-wise
441 < potentials for non-bonded interactions between Au and non-metal sites.
436 > The potentials developed from {\it ab initio} calculations by Leng
437 > {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
438 > interactions between Au and aromatic C/H atoms in toluene. However,
439 > the Lennard-Jones parameters between Au and other types of particles,
440 > (e.g. AA alkanes) have not yet been established. For these
441 > interactions, the Lorentz-Berthelot mixing rule can be used to derive
442 > effective single-atom LJ parameters for the metal using the fit values
443 > for toluene. These are then used to construct reasonable mixing
444 > parameters for the interactions between the gold and other atoms.
445 > Table \ref{MnM} summarizes the ``metal/non-metal'' parameters used in
446 > our simulations.
447  
451 However, the Lennard-Jones parameters between Au and other types of
452 particles, such as All-Atom normal alkanes in our simulations are not
453 yet well-established. For these interactions, we attempt to derive
454 their parameters using the Mixing Rule. To do this, Au pseudo
455 Lennard-Jones parameters for ``Metal-non-Metal'' (MnM) interactions
456 were first extracted from the Au-CH$_x$ parameters by applying the
457 Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
458 parameters in our simulations.
459
448   \begin{table*}
449    \begin{minipage}{\linewidth}
450      \begin{center}
# Line 492 | Line 480 | parameters in our simulations.
480    \end{minipage}
481   \end{table*}
482  
495 \subsection{Vibrational Spectrum}
496 To investigate the mechanism of interfacial thermal conductance, the
497 vibrational spectrum is utilized as a complementary tool. Vibrational
498 spectra were taken for individual components in different
499 simulations. To obtain these spectra, simulations were run after
500 equilibration, in the NVE ensemble. Snapshots of configurations were
501 collected at a frequency that is higher than that of the fastest
502 vibrations occuring in the simulations. With these configurations, the
503 velocity auto-correlation functions can be computed:
504 \begin{equation}
505 C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
506 \label{vCorr}
507 \end{equation}
508 Followed by Fourier transforms, the power spectrum can be constructed:
509 \begin{equation}
510 \hat{f}(\omega) = \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
511 \label{fourier}
512 \end{equation}
483  
484 < \section{Results and Discussions}
485 < In what follows, how the parameters and protocol of simulations would
486 < affect the measurement of $G$'s is first discussed. With a reliable
487 < protocol and set of parameters, the influence of capping agent
488 < coverage on thermal conductance is investigated. Besides, different
489 < force field models for both solvents and selected deuterated models
490 < were tested and compared. Finally, a summary of the role of capping
491 < agent in the interfacial thermal transport process is given.
484 > \section{Results}
485 > There are many factors contributing to the measured interfacial
486 > conductance; some of these factors are physically motivated
487 > (e.g. coverage of the surface by the capping agent coverage and
488 > solvent identity), while some are governed by parameters of the
489 > methodology (e.g. applied flux and the formulas used to obtain the
490 > conductance). In this section we discuss the major physical and
491 > calculational effects on the computed conductivity.
492  
493 < \subsection{How Simulation Parameters Affects $G$}
524 < We have varied our protocol or other parameters of the simulations in
525 < order to investigate how these factors would affect the measurement of
526 < $G$'s. It turned out that while some of these parameters would not
527 < affect the results substantially, some other changes to the
528 < simulations would have a significant impact on the measurement
529 < results.
493 > \subsection{Effects due to capping agent coverage}
494  
495 < In some of our simulations, we allowed $L_x$ and $L_y$ to change
496 < during equilibrating the liquid phase. Due to the stiffness of the
497 < crystalline Au structure, $L_x$ and $L_y$ would not change noticeably
498 < after equilibration. Although $L_z$ could fluctuate $\sim$1\% after a
499 < system is fully equilibrated in the NPT ensemble, this fluctuation, as
500 < well as those of $L_x$ and $L_y$ (which is significantly smaller),
537 < would not be magnified on the calculated $G$'s, as shown in Table
538 < \ref{AuThiolHexaneUA}. This insensivity to $L_i$ fluctuations allows
539 < reliable measurement of $G$'s without the necessity of extremely
540 < cautious equilibration process.
495 > A series of different initial conditions with a range of surface
496 > coverages was prepared and solvated with various with both of the
497 > solvent molecules. These systems were then equilibrated and their
498 > interfacial thermal conductivity was measured with the NIVS
499 > algorithm. Figure \ref{coverage} demonstrates the trend of conductance
500 > with respect to surface coverage.
501  
502 + \begin{figure}
503 + \includegraphics[width=\linewidth]{coverage}
504 + \caption{Comparison of interfacial thermal conductivity ($G$) values
505 +  for the Au-butanethiol/solvent interface with various UA models and
506 +  different capping agent coverages at $\langle T\rangle\sim$200K.}
507 + \label{coverage}
508 + \end{figure}
509 +
510 + 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.
511 +
512 + 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.
513 +
514 + 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.
515 +
516 + 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.
517 +
518 + 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.
519 +
520 + 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.
521 +
522 + \subsection{Effects due to Solvent \& Solvent Models}
523 + 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.
524 +
525 + \begin{table*}
526 +  \begin{minipage}{\linewidth}
527 +    \begin{center}
528 +      
529 +      \caption{Computed interfacial thermal conductance ($G$ and
530 +        $G^\prime$) values for interfaces using various models for
531 +        solvent and capping agent (or without capping agent) at
532 +        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
533 +        or capping agent molecules; ``Avg.'' denotes results that are
534 +        averages of simulations under different applied thermal flux values $(J_z)$. Error
535 +        estimates are indicated in parentheses.)}
536 +      
537 +      \begin{tabular}{llccc}
538 +        \hline\hline
539 +        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
540 +        (or bare surface) & model & (GW/m$^2$) &
541 +        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
542 +        \hline
543 +        UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
544 +              & UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
545 +              & AA hexane    & Avg. & 131(6)    & 122(10)   \\
546 +              & UA toluene   & 1.96 & 187(16)   & 151(11)   \\
547 +              & AA toluene   & 1.89 & 200(36)   & 149(53)   \\
548 +        \hline
549 +        AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
550 +              & AA hexane    & Avg. & 442(14)   & 356(31)   \\
551 +              & AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
552 +              & UA toluene   & 1.98 & 125(25)   & 97(60)    \\
553 +              & AA toluene   & 3.79 & 487(56)   & 290(42)   \\
554 +        \hline
555 +        AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
556 +              & AA hexane    & 1.92 & 243(29)   & 191(11)   \\
557 +              & AA toluene   & 1.93 & 364(36)   & 322(67)   \\
558 +        \hline
559 +        bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
560 +              & UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
561 +              & AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
562 +              & UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
563 +        \hline\hline
564 +      \end{tabular}
565 +      \label{modelTest}
566 +    \end{center}
567 +  \end{minipage}
568 + \end{table*}
569 +
570 + 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.  
571 +
572 + 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.
573 +
574 + 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.
575 +
576 + 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.
577 +
578 + \begin{figure}
579 + \includegraphics[width=\linewidth]{aahxntln}
580 + \caption{Spectra obtained for all-atom (AA) Au / butanethiol / solvent
581 +  systems. When butanethiol is deuterated (lower left), its
582 +  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).}
583 + \label{aahxntln}
584 + \end{figure}
585 +
586 + 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.
587 +
588 + \subsubsection{Are electronic excitations in the metal important?}
589 + 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$.
590 +
591 + \subsection{Effects due to methodology and simulation parameters}
592 +
593 + START HERE
594 +
595 + We have varied our protocol or other parameters of the simulations in order to investigate how these factors would affect the computation of $G$.
596 +
597 + 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.
598 +
599   As stated in our computational details, the spacing filled with
600   solvent molecules can be chosen within a range. This allows some
601   change of solvent molecule numbers for the same Au-butanethiol
# Line 549 | Line 606 | Our NIVS algorithm allows change of unphysical thermal
606   smaller system size would be preferable, given that the liquid phase
607   structure is not affected.
608  
609 + \subsubsection{Effects of applied flux}
610   Our NIVS algorithm allows change of unphysical thermal flux both in
611   direction and in quantity. This feature extends our investigation of
612   interfacial thermal conductance. However, the magnitude of this
# Line 612 | Line 670 | Furthermore, we also attempted to increase system aver
670    \end{minipage}
671   \end{table*}
672  
673 + \subsubsection{Effects due to average temperature}
674 +
675   Furthermore, we also attempted to increase system average temperatures
676   to above 200K. These simulations are first equilibrated in the NPT
677   ensemble under normal pressure. As stated above, the TraPPE-UA model
# Line 696 | Line 756 | affected by surface reconstructions.
756   interfaces even at a relatively high temperature without being
757   affected by surface reconstructions.
758  
699 \subsection{Influence of Capping Agent Coverage on $G$}
700 To investigate the influence of butanethiol coverage on interfacial
701 thermal conductance, a series of different coverage Au-butanethiol
702 surfaces is prepared and solvated with various organic
703 molecules. These systems are then equilibrated and their interfacial
704 thermal conductivity are measured with our NIVS algorithm. Figure
705 \ref{coverage} demonstrates the trend of conductance change with
706 respect to different coverages of butanethiol. To study the isotope
707 effect in interfacial thermal conductance, deuterated UA-hexane is
708 included as well.
759  
760 < \begin{figure}
711 < \includegraphics[width=\linewidth]{coverage}
712 < \caption{Comparison of interfacial thermal conductivity ($G$) values
713 <  for the Au-butanethiol/solvent interface with various UA models and
714 <  different capping agent coverages at $\langle T\rangle\sim$200K
715 <  using certain energy flux respectively.}
716 < \label{coverage}
717 < \end{figure}
760 > \section{Discussion}
761  
762 < It turned out that with partial covered butanethiol on the Au(111)
763 < surface, the derivative definition for $G^\prime$
764 < (Eq. \ref{derivativeG}) was difficult to apply, due to the difficulty
765 < in locating the maximum of change of $\lambda$. Instead, the discrete
766 < definition (Eq. \ref{discreteG}) is easier to apply, as the Gibbs
767 < deviding surface can still be well-defined. Therefore, $G$ (not
768 < $G^\prime$) was used for this section.
762 > \subsection{Capping agent acts as a vibrational coupler between solid
763 >  and solvent phases}
764 > To investigate the mechanism of interfacial thermal conductance, the
765 > vibrational power spectrum was computed. Power spectra were taken for
766 > individual components in different simulations. To obtain these
767 > spectra, simulations were run after equilibration, in the NVE
768 > ensemble, and without a thermal gradient. Snapshots of configurations
769 > were collected at a frequency that is higher than that of the fastest
770 > vibrations occuring in the simulations. With these configurations, the
771 > velocity auto-correlation functions can be computed:
772 > \begin{equation}
773 > C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
774 > \label{vCorr}
775 > \end{equation}
776 > The power spectrum is constructed via a Fourier transform of the
777 > symmetrized velocity autocorrelation function,
778 > \begin{equation}
779 >  \hat{f}(\omega) =
780 >  \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
781 > \label{fourier}
782 > \end{equation}
783  
727 From Figure \ref{coverage}, one can see the significance of the
728 presence of capping agents. Even when a fraction of the Au(111)
729 surface sites are covered with butanethiols, the conductivity would
730 see an enhancement by at least a factor of 3. This indicates the
731 important role cappping agent is playing for thermal transport
732 phenomena on metal / organic solvent surfaces.
784  
785 < Interestingly, as one could observe from our results, the maximum
735 < conductance enhancement (largest $G$) happens while the surfaces are
736 < about 75\% covered with butanethiols. This again indicates that
737 < solvent-capping agent contact has an important role of the thermal
738 < transport process. Slightly lower butanethiol coverage allows small
739 < gaps between butanethiols to form. And these gaps could be filled with
740 < solvent molecules, which acts like ``heat conductors'' on the
741 < surface. The higher degree of interaction between these solvent
742 < molecules and capping agents increases the enhancement effect and thus
743 < produces a higher $G$ than densely packed butanethiol arrays. However,
744 < once this maximum conductance enhancement is reached, $G$ decreases
745 < when butanethiol coverage continues to decrease. Each capping agent
746 < molecule reaches its maximum capacity for thermal
747 < conductance. Therefore, even higher solvent-capping agent contact
748 < would not offset this effect. Eventually, when butanethiol coverage
749 < continues to decrease, solvent-capping agent contact actually
750 < decreases with the disappearing of butanethiol molecules. In this
751 < case, $G$ decrease could not be offset but instead accelerated. [NEED
752 < SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
753 <
754 < A comparison of the results obtained from differenet organic solvents
755 < can also provide useful information of the interfacial thermal
756 < transport process. The deuterated hexane (UA) results do not appear to
757 < be much different from those of normal hexane (UA), given that
758 < butanethiol (UA) is non-deuterated for both solvents. These UA model
759 < studies, even though eliminating C-H vibration samplings, still have
760 < C-C vibrational frequencies different from each other. However, these
761 < differences in the infrared range do not seem to produce an observable
762 < difference for the results of $G$ (Figure \ref{uahxnua}).
763 <
764 < \begin{figure}
765 < \includegraphics[width=\linewidth]{uahxnua}
766 < \caption{Vibrational spectra obtained for normal (upper) and
767 <  deuterated (lower) hexane in Au-butanethiol/hexane
768 <  systems. Butanethiol spectra are shown as reference. Both hexane and
769 <  butanethiol were using United-Atom models.}
770 < \label{uahxnua}
771 < \end{figure}
772 <
773 < Furthermore, results for rigid body toluene solvent, as well as other
774 < UA-hexane solvents, are reasonable within the general experimental
775 < ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
776 < suggests that explicit hydrogen might not be a required factor for
777 < modeling thermal transport phenomena of systems such as
778 < Au-thiol/organic solvent.
779 <
780 < However, results for Au-butanethiol/toluene do not show an identical
781 < trend with those for Au-butanethiol/hexane in that $G$ remains at
782 < approximately the same magnitue when butanethiol coverage differs from
783 < 25\% to 75\%. This might be rooted in the molecule shape difference
784 < for planar toluene and chain-like {\it n}-hexane. Due to this
785 < difference, toluene molecules have more difficulty in occupying
786 < relatively small gaps among capping agents when their coverage is not
787 < too low. Therefore, the solvent-capping agent contact may keep
788 < increasing until the capping agent coverage reaches a relatively low
789 < level. This becomes an offset for decreasing butanethiol molecules on
790 < its effect to the process of interfacial thermal transport. Thus, one
791 < can see a plateau of $G$ vs. butanethiol coverage in our results.
792 <
793 < \subsection{Influence of Chosen Molecule Model on $G$}
794 < In addition to UA solvent/capping agent models, AA models are included
795 < in our simulations as well. Besides simulations of the same (UA or AA)
796 < model for solvent and capping agent, different models can be applied
797 < to different components. Furthermore, regardless of models chosen,
798 < either the solvent or the capping agent can be deuterated, similar to
799 < the previous section. Table \ref{modelTest} summarizes the results of
800 < these studies.
801 <
802 < \begin{table*}
803 <  \begin{minipage}{\linewidth}
804 <    \begin{center}
805 <      
806 <      \caption{Computed interfacial thermal conductivity ($G$ and
807 <        $G^\prime$) values for interfaces using various models for
808 <        solvent and capping agent (or without capping agent) at
809 <        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
810 <        or capping agent molecules; ``Avg.'' denotes results that are
811 <        averages of simulations under different $J_z$'s. Error
812 <        estimates indicated in parenthesis.)}
813 <      
814 <      \begin{tabular}{llccc}
815 <        \hline\hline
816 <        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
817 <        (or bare surface) & model & (GW/m$^2$) &
818 <        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
819 <        \hline
820 <        UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
821 <              & UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
822 <              & AA hexane    & Avg. & 131(6)    & 122(10)   \\
823 <              & UA toluene   & 1.96 & 187(16)   & 151(11)   \\
824 <              & AA toluene   & 1.89 & 200(36)   & 149(53)   \\
825 <        \hline
826 <        AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
827 <              & AA hexane    & Avg. & 442(14)   & 356(31)   \\
828 <              & AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
829 <              & UA toluene   & 1.98 & 125(25)   & 97(60)    \\
830 <              & AA toluene   & 3.79 & 487(56)   & 290(42)   \\
831 <        \hline
832 <        AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
833 <              & AA hexane    & 1.92 & 243(29)   & 191(11)   \\
834 <              & AA toluene   & 1.93 & 364(36)   & 322(67)   \\
835 <        \hline
836 <        bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
837 <              & UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
838 <              & AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
839 <              & UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
840 <        \hline\hline
841 <      \end{tabular}
842 <      \label{modelTest}
843 <    \end{center}
844 <  \end{minipage}
845 < \end{table*}
846 <
847 < To facilitate direct comparison, the same system with differnt models
848 < for different components uses the same length scale for their
849 < simulation cells. Without the presence of capping agent, using
850 < different models for hexane yields similar results for both $G$ and
851 < $G^\prime$, and these two definitions agree with eath other very
852 < well. This indicates very weak interaction between the metal and the
853 < solvent, and is a typical case for acoustic impedance mismatch between
854 < these two phases.
855 <
856 < As for Au(111) surfaces completely covered by butanethiols, the choice
857 < of models for capping agent and solvent could impact the measurement
858 < of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
859 < interfaces, using AA model for both butanethiol and hexane yields
860 < substantially higher conductivity values than using UA model for at
861 < least one component of the solvent and capping agent, which exceeds
862 < the general range of experimental measurement results. This is
863 < probably due to the classically treated C-H vibrations in the AA
864 < model, which should not be appreciably populated at normal
865 < temperatures. In comparison, once either the hexanes or the
866 < butanethiols are deuterated, one can see a significantly lower $G$ and
867 < $G^\prime$. In either of these cases, the C-H(D) vibrational overlap
868 < between the solvent and the capping agent is removed (Figure
869 < \ref{aahxntln}). Conclusively, the improperly treated C-H vibration in
870 < the AA model produced over-predicted results accordingly. Compared to
871 < the AA model, the UA model yields more reasonable results with higher
872 < computational efficiency.
873 <
874 < \begin{figure}
875 < \includegraphics[width=\linewidth]{aahxntln}
876 < \caption{Spectra obtained for All-Atom model Au-butanethil/solvent
877 <  systems. When butanethiol is deuterated (lower left), its
878 <  vibrational overlap with hexane would decrease significantly,
879 <  compared with normal butanethiol (upper left). However, this
880 <  dramatic change does not apply to toluene as much (right).}
881 < \label{aahxntln}
882 < \end{figure}
883 <
884 < However, for Au-butanethiol/toluene interfaces, having the AA
885 < butanethiol deuterated did not yield a significant change in the
886 < measurement results. Compared to the C-H vibrational overlap between
887 < hexane and butanethiol, both of which have alkyl chains, that overlap
888 < between toluene and butanethiol is not so significant and thus does
889 < not have as much contribution to the heat exchange
890 < process. Conversely, extra degrees of freedom such as the C-H
891 < vibrations could yield higher heat exchange rate between these two
892 < phases and result in a much higher conductivity.
893 <
894 < Although the QSC model for Au is known to predict an overly low value
895 < for bulk metal gold conductivity\cite{kuang:164101}, our computational
896 < results for $G$ and $G^\prime$ do not seem to be affected by this
897 < drawback of the model for metal. Instead, our results suggest that the
898 < modeling of interfacial thermal transport behavior relies mainly on
899 < the accuracy of the interaction descriptions between components
900 < occupying the interfaces.
901 <
902 < \subsection{Role of Capping Agent in Interfacial Thermal Conductance}
785 > \subsubsection{The role of specific vibrations}
786   The vibrational spectra for gold slabs in different environments are
787   shown as in Figure \ref{specAu}. Regardless of the presence of
788   solvent, the gold surfaces covered by butanethiol molecules, compared
# Line 910 | Line 793 | Simultaneously, the vibrational overlap between butane
793   simulations, the Au/S interfaces do not appear major heat barriers
794   compared to the butanethiol / solvent interfaces.
795  
796 + \subsubsection{Overlap of power spectrum}
797   Simultaneously, the vibrational overlap between butanethiol and
798   organic solvents suggests higher thermal exchange efficiency between
799   these two components. Even exessively high heat transport was observed
# Line 933 | Line 817 | capping agents.
817   \label{specAu}
818   \end{figure}
819  
820 < [MAY ADD COMPARISON OF AU SLAB WIDTHS]
820 > \subsubsection{Isotopic substitution and vibrational overlap}
821 > A comparison of the results obtained from the two different organic
822 > solvents can also provide useful information of the interfacial
823 > thermal transport process. The deuterated hexane (UA) results do not
824 > appear to be substantially different from those of normal hexane (UA),
825 > given that butanethiol (UA) is non-deuterated for both solvents. The
826 > UA models, even though they have eliminated C-H vibrational overlap,
827 > still have significant overlap in the infrared spectra.  Because
828 > differences in the infrared range do not seem to produce an observable
829 > difference for the results of $G$ (Figure \ref{uahxnua}).
830  
831 + \begin{figure}
832 + \includegraphics[width=\linewidth]{uahxnua}
833 + \caption{Vibrational spectra obtained for normal (upper) and
834 +  deuterated (lower) hexane in Au-butanethiol/hexane
835 +  systems. Butanethiol spectra are shown as reference. Both hexane and
836 +  butanethiol were using United-Atom models.}
837 + \label{uahxnua}
838 + \end{figure}
839 +
840   \section{Conclusions}
841   The NIVS algorithm we developed has been applied to simulations of
842   Au-butanethiol surfaces with organic solvents. This algorithm allows
# Line 974 | Line 876 | Au(111) surface\cite{vlugt:cpc2007154}. This differenc
876  
877   Vlugt {\it et al.} has investigated the surface thiol structures for
878   nanocrystal gold and pointed out that they differs from those of the
879 < Au(111) surface\cite{vlugt:cpc2007154}. This difference might lead to
880 < change of interfacial thermal transport behavior as well. To
881 < investigate this problem, an effective means to introduce thermal flux
882 < and measure the corresponding thermal gradient is desirable for
883 < simulating structures with spherical symmetry.
879 > Au(111) surface\cite{landman:1998,vlugt:cpc2007154}. This difference
880 > might lead to change of interfacial thermal transport behavior as
881 > well. To investigate this problem, an effective means to introduce
882 > thermal flux and measure the corresponding thermal gradient is
883 > desirable for simulating structures with spherical symmetry.
884  
885   \section{Acknowledgments}
886   Support for this project was provided by the National Science
887   Foundation under grant CHE-0848243. Computational time was provided by
888   the Center for Research Computing (CRC) at the University of Notre
889 < Dame. \newpage
889 > Dame.
890 > \newpage
891  
892   \bibliography{interfacial}
893  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines