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 3718 by gezelter, Thu Jan 27 16:30:51 2011 UTC vs.
Revision 3751 by gezelter, Tue Jul 26 19:43:10 2011 UTC

# Line 22 | Line 22
22   \setlength{\abovecaptionskip}{20 pt}
23   \setlength{\belowcaptionskip}{30 pt}
24  
25 < %\renewcommand\citemid{\ } % no comma in optional referenc note
26 < \bibpunct{[}{]}{,}{s}{}{;}
27 < \bibliographystyle{aip}
25 > %\renewcommand\citemid{\ } % no comma in optional reference note
26 > \bibpunct{[}{]}{,}{n}{}{;}
27 > \bibliographystyle{achemso}
28  
29   \begin{document}
30  
# Line 44 | Line 44 | The abstract version 2
44   \begin{doublespace}
45  
46   \begin{abstract}
47 < The abstract version 2
47 >
48 > With the Non-Isotropic Velocity Scaling algorithm (NIVS) we have
49 > developed, an unphysical thermal flux can be effectively set up even
50 > for non-homogeneous systems like interfaces in non-equilibrium
51 > molecular dynamics simulations. In this work, this algorithm is
52 > applied for simulating thermal conductance at metal / organic solvent
53 > interfaces with various coverages of butanethiol capping
54 > agents. Different solvents and force field models were tested. Our
55 > results suggest that the United-Atom models are able to provide an
56 > estimate of the interfacial thermal conductivity comparable to
57 > experiments in our simulations with satisfactory computational
58 > efficiency. From our results, the acoustic impedance mismatch between
59 > metal and liquid phase is effectively reduced by the capping
60 > agents, and thus leads to interfacial thermal conductance
61 > enhancement. Furthermore, this effect is closely related to the
62 > capping agent coverage on the metal surfaces and the type of solvent
63 > molecules, and is affected by the models used in the simulations.
64 >
65   \end{abstract}
66  
67   \newpage
# Line 56 | Line 73 | The abstract version 2
73   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74  
75   \section{Introduction}
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 < The intro.
85 > Experimentally, various interfaces have been investigated for their
86 > thermal conductance. Wang {\it et al.} studied heat transport through
87 > long-chain hydrocarbon monolayers on gold substrate at individual
88 > molecular level,\cite{Wang10082007} Schmidt {\it et al.} studied the
89 > role of CTAB on thermal transport between gold nanorods and
90 > solvent,\cite{doi:10.1021/jp8051888} and Juv\'e {\it et al.} studied
91 > the cooling dynamics, which is controlled by thermal interface
92 > resistence of glass-embedded metal
93 > nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
94 > normally considered barriers for heat transport, Alper {\it et al.}
95 > suggested that specific ligands (capping agents) could completely
96 > eliminate this barrier
97 > ($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}
98  
99 + Theoretical and computational models have also been used to study the
100 + interfacial thermal transport in order to gain an understanding of
101 + this phenomena at the molecular level. Recently, Hase and coworkers
102 + employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
103 + study thermal transport from hot Au(111) substrate to a self-assembled
104 + monolayer of alkylthiol with relatively long chain (8-20 carbon
105 + atoms).\cite{hase:2010,hase:2011} However, ensemble averaged
106 + measurements for heat conductance of interfaces between the capping
107 + monolayer on Au and a solvent phase have yet to be studied with their
108 + approach. The comparatively low thermal flux through interfaces is
109 + difficult to measure with Equilibrium MD or forward NEMD simulation
110 + methods. Therefore, the Reverse NEMD (RNEMD)
111 + methods\cite{MullerPlathe:1997xw,kuang:164101} would have the
112 + advantage of applying this difficult to measure flux (while measuring
113 + the resulting gradient), given that the simulation methods being able
114 + to effectively apply an unphysical flux in non-homogeneous systems.
115 + Garde and coworkers\cite{garde:nl2005,garde:PhysRevLett2009} applied
116 + this approach to various liquid interfaces and studied how thermal
117 + conductance (or resistance) is dependent on chemistry details of
118 + interfaces, e.g. hydrophobic and hydrophilic aqueous interfaces.
119 +
120 + Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS)
121 + algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
122 + retains the desirable features of RNEMD (conservation of linear
123 + momentum and total energy, compatibility with periodic boundary
124 + conditions) while establishing true thermal distributions in each of
125 + the two slabs. Furthermore, it allows effective thermal exchange
126 + between particles of different identities, and thus makes the study of
127 + interfacial conductance much simpler.
128 +
129 + The work presented here deals with the Au(111) surface covered to
130 + varying degrees by butanethiol, a capping agent with short carbon
131 + chain, and solvated with organic solvents of different molecular
132 + properties. Different models were used for both the capping agent and
133 + the solvent force field parameters. Using the NIVS algorithm, the
134 + thermal transport across these interfaces was studied and the
135 + underlying mechanism for the phenomena was investigated.
136 +
137 + \section{Methodology}
138 + \subsection{Imposd-Flux Methods in MD Simulations}
139 + Steady state MD simulations have an advantage in that not many
140 + trajectories are needed to study the relationship between thermal flux
141 + and thermal gradients. For systems with low interfacial conductance,
142 + one must have a method capable of generating or measuring relatively
143 + small fluxes, compared to those required for bulk conductivity. This
144 + requirement makes the calculation even more difficult for
145 + slowly-converging equilibrium methods.\cite{Viscardy:2007lq} Forward
146 + NEMD methods impose a gradient (and measure a flux), but at interfaces
147 + it is not clear what behavior should be imposed at the boundaries
148 + between materials.  Imposed-flux reverse non-equilibrium
149 + methods\cite{MullerPlathe:1997xw} set the flux {\it a priori} and
150 + the thermal response becomes an easy-to-measure quantity.  Although
151 + M\"{u}ller-Plathe's original momentum swapping approach can be used
152 + for exchanging energy between particles of different identity, the
153 + kinetic energy transfer efficiency is affected by the mass difference
154 + between the particles, which limits its application on heterogeneous
155 + interfacial systems.
156 +
157 + The non-isotropic velocity scaling (NIVS) \cite{kuang:164101} approach
158 + to non-equilibrium MD simulations is able to impose a wide range of
159 + kinetic energy fluxes without obvious perturbation to the velocity
160 + distributions of the simulated systems. Furthermore, this approach has
161 + the advantage in heterogeneous interfaces in that kinetic energy flux
162 + can be applied between regions of particles of arbitary identity, and
163 + the flux will not be restricted by difference in particle mass.
164 +
165 + The NIVS algorithm scales the velocity vectors in two separate regions
166 + of a simulation system with respective diagonal scaling matricies. To
167 + determine these scaling factors in the matricies, a set of equations
168 + including linear momentum conservation and kinetic energy conservation
169 + constraints and target energy flux satisfaction is solved. With the
170 + scaling operation applied to the system in a set frequency, bulk
171 + temperature gradients can be easily established, and these can be used
172 + for computing thermal conductivities. The NIVS algorithm conserves
173 + momenta and energy and does not depend on an external thermostat.
174 +
175 + \subsection{Defining Interfacial Thermal Conductivity ($G$)}
176 +
177 + For an interface with relatively low interfacial conductance, and a
178 + thermal flux between two distinct bulk regions, the regions on either
179 + side of the interface rapidly come to a state in which the two phases
180 + have relatively homogeneous (but distinct) temperatures. The
181 + interfacial thermal conductivity $G$ can therefore be approximated as:
182 + \begin{equation}
183 +  G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
184 +    \langle T_\mathrm{cold}\rangle \right)}
185 + \label{lowG}
186 + \end{equation}
187 + where ${E_{total}}$ is the total imposed non-physical kinetic energy
188 + transfer during the simulation and ${\langle T_\mathrm{hot}\rangle}$
189 + and ${\langle T_\mathrm{cold}\rangle}$ are the average observed
190 + temperature of the two separated phases.
191 +
192 + When the interfacial conductance is {\it not} small, there are two
193 + ways to define $G$. One way is to assume the temperature is discrete
194 + on the two sides of the interface. $G$ can be calculated using the
195 + applied thermal flux $J$ and the maximum temperature difference
196 + measured along the thermal gradient max($\Delta T$), which occurs at
197 + the Gibbs deviding surface (Figure \ref{demoPic}): \begin{equation}
198 +  G=\frac{J}{\Delta T} \label{discreteG} \end{equation}
199 +
200 + \begin{figure}
201 + \includegraphics[width=\linewidth]{method}
202 + \caption{Interfacial conductance can be calculated by applying an
203 +  (unphysical) kinetic energy flux between two slabs, one located
204 +  within the metal and another on the edge of the periodic box.  The
205 +  system responds by forming a thermal response or a gradient.  In
206 +  bulk liquids, this gradient typically has a single slope, but in
207 +  interfacial systems, there are distinct thermal conductivity
208 +  domains.  The interfacial conductance, $G$ is found by measuring the
209 +  temperature gap at the Gibbs dividing surface, or by using second
210 +  derivatives of the thermal profile.}
211 + \label{demoPic}
212 + \end{figure}
213 +
214 + The other approach is to assume a continuous temperature profile along
215 + the thermal gradient axis (e.g. $z$) and define $G$ at the point where
216 + the magnitude of thermal conductivity ($\lambda$) change reaches its
217 + maximum, given that $\lambda$ is well-defined throughout the space:
218 + \begin{equation}
219 + G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
220 +         = \Big|\frac{\partial}{\partial z}\left(-J_z\Big/
221 +           \left(\frac{\partial T}{\partial z}\right)\right)\Big|
222 +         = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
223 +         \Big/\left(\frac{\partial T}{\partial z}\right)^2
224 + \label{derivativeG}
225 + \end{equation}
226 +
227 + With temperature profiles obtained from simulation, one is able to
228 + approximate the first and second derivatives of $T$ with finite
229 + difference methods and calculate $G^\prime$. In what follows, both
230 + definitions have been used, and are compared in the results.
231 +
232 + To investigate the interfacial conductivity at metal / solvent
233 + interfaces, we have modeled a metal slab with its (111) surfaces
234 + perpendicular to the $z$-axis of our simulation cells. The metal slab
235 + has been prepared both with and without capping agents on the exposed
236 + surface, and has been solvated with simple organic solvents, as
237 + illustrated in Figure \ref{gradT}.
238 +
239 + With the simulation cell described above, we are able to equilibrate
240 + the system and impose an unphysical thermal flux between the liquid
241 + and the metal phase using the NIVS algorithm. By periodically applying
242 + the unphysical flux, we obtained a temperature profile and its spatial
243 + derivatives. Figure \ref{gradT} shows how an applied thermal flux can
244 + be used to obtain the 1st and 2nd derivatives of the temperature
245 + profile.
246 +
247 + \begin{figure}
248 + \includegraphics[width=\linewidth]{gradT}
249 + \caption{A sample of Au-butanethiol/hexane interfacial system and the
250 +  temperature profile after a kinetic energy flux is imposed to
251 +  it. The 1st and 2nd derivatives of the temperature profile can be
252 +  obtained with finite difference approximation (lower panel).}
253 + \label{gradT}
254 + \end{figure}
255 +
256 + \section{Computational Details}
257 + \subsection{Simulation Protocol}
258 + The NIVS algorithm has been implemented in our MD simulation code,
259 + OpenMD\cite{Meineke:2005gd,openmd}, and was used for our simulations.
260 + Metal slabs of 6 or 11 layers of Au atoms were first equilibrated
261 + under atmospheric pressure (1 atm) and 200K. After equilibration,
262 + butanethiol capping agents were placed at three-fold hollow sites on
263 + the Au(111) surfaces. These sites are either {\it fcc} or {\it
264 +  hcp} sites, although Hase {\it et al.} found that they are
265 + equivalent in a heat transfer process,\cite{hase:2010} so we did not
266 + distinguish between these sites in our study. The maximum butanethiol
267 + capacity on Au surface is $1/3$ of the total number of surface Au
268 + atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
269 + structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
270 + series of lower coverages was also prepared by eliminating
271 + butanethiols from the higher coverage surface in a regular manner. The
272 + lower coverages were prepared in order to study the relation between
273 + coverage and interfacial conductance.
274 +
275 + The capping agent molecules were allowed to migrate during the
276 + simulations. They distributed themselves uniformly and sampled a
277 + number of three-fold sites throughout out study. Therefore, the
278 + initial configuration does not noticeably affect the sampling of a
279 + variety of configurations of the same coverage, and the final
280 + conductance measurement would be an average effect of these
281 + configurations explored in the simulations.
282 +
283 + After the modified Au-butanethiol surface systems were equilibrated in
284 + the canonical (NVT) ensemble, organic solvent molecules were packed in
285 + the previously empty part of the simulation cells.\cite{packmol} Two
286 + solvents were investigated, one which has little vibrational overlap
287 + with the alkanethiol and which has a planar shape (toluene), and one
288 + which has similar vibrational frequencies to the capping agent and
289 + chain-like shape ({\it n}-hexane).
290 +
291 + The simulation cells were not particularly extensive along the
292 + $z$-axis, as a very long length scale for the thermal gradient may
293 + cause excessively hot or cold temperatures in the middle of the
294 + solvent region and lead to undesired phenomena such as solvent boiling
295 + or freezing when a thermal flux is applied. Conversely, too few
296 + solvent molecules would change the normal behavior of the liquid
297 + phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
298 + these extreme cases did not happen to our simulations. The spacing
299 + between periodic images of the gold interfaces is $35 \sim 75$\AA.
300 +
301 + The initial configurations generated are further equilibrated with the
302 + $x$ and $y$ dimensions fixed, only allowing the $z$-length scale to
303 + change. This is to ensure that the equilibration of liquid phase does
304 + not affect the metal's crystalline structure. Comparisons were made
305 + with simulations that allowed changes of $L_x$ and $L_y$ during NPT
306 + equilibration. No substantial changes in the box geometry were noticed
307 + in these simulations. After ensuring the liquid phase reaches
308 + equilibrium at atmospheric pressure (1 atm), further equilibration was
309 + carried out under canonical (NVT) and microcanonical (NVE) ensembles.
310 +
311 + After the systems reach equilibrium, NIVS was used to impose an
312 + unphysical thermal flux between the metal and the liquid phases. Most
313 + of our simulations were done under an average temperature of
314 + $\sim$200K. Therefore, thermal flux usually came from the metal to the
315 + liquid so that the liquid has a higher temperature and would not
316 + freeze due to lowered temperatures. After this induced temperature
317 + gradient had stablized, the temperature profile of the simulation cell
318 + was recorded. To do this, the simulation cell is devided evenly into
319 + $N$ slabs along the $z$-axis. The average temperatures of each slab
320 + are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
321 + the same, the derivatives of $T$ with respect to slab number $n$ can
322 + be directly used for $G^\prime$ calculations: \begin{equation}
323 +  G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
324 +         \Big/\left(\frac{\partial T}{\partial z}\right)^2
325 +         = |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
326 +         \Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
327 +         = |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big|
328 +         \Big/\left(\frac{\partial T}{\partial n}\right)^2
329 + \label{derivativeG2}
330 + \end{equation}
331 +
332 + All of the above simulation procedures use a time step of 1 fs. Each
333 + equilibration stage took a minimum of 100 ps, although in some cases,
334 + longer equilibration stages were utilized.
335 +
336 + \subsection{Force Field Parameters}
337 + Our simulations include a number of chemically distinct components.
338 + Figure \ref{demoMol} demonstrates the sites defined for both
339 + United-Atom and All-Atom models of the organic solvent and capping
340 + agents in our simulations. Force field parameters are needed for
341 + interactions both between the same type of particles and between
342 + particles of different species.
343 +
344 + \begin{figure}
345 + \includegraphics[width=\linewidth]{structures}
346 + \caption{Structures of the capping agent and solvents utilized in
347 +  these simulations. The chemically-distinct sites (a-e) are expanded
348 +  in terms of constituent atoms for both United Atom (UA) and All Atom
349 +  (AA) force fields.  Most parameters are from
350 +  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} (UA) and
351 +  \protect\cite{OPLSAA} (AA).  Cross-interactions with the Au atoms are given
352 +  in Table \ref{MnM}.}
353 + \label{demoMol}
354 + \end{figure}
355 +
356 + The Au-Au interactions in metal lattice slab is described by the
357 + quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
358 + potentials include zero-point quantum corrections and are
359 + reparametrized for accurate surface energies compared to the
360 + Sutton-Chen potentials.\cite{Chen90}
361 +
362 + For the two solvent molecules, {\it n}-hexane and toluene, two
363 + different atomistic models were utilized. Both solvents were modeled
364 + using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
365 + parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
366 + for our UA solvent molecules. In these models, sites are located at
367 + the carbon centers for alkyl groups. Bonding interactions, including
368 + bond stretches and bends and torsions, were used for intra-molecular
369 + sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones
370 + potentials are used.
371 +
372 + By eliminating explicit hydrogen atoms, the TraPPE-UA models are
373 + simple and computationally efficient, while maintaining good accuracy.
374 + However, the TraPPE-UA model for alkanes is known to predict a slighly
375 + lower boiling point than experimental values. This is one of the
376 + reasons we used a lower average temperature (200K) for our
377 + simulations. If heat is transferred to the liquid phase during the
378 + NIVS simulation, the liquid in the hot slab can actually be
379 + substantially warmer than the mean temperature in the simulation. The
380 + lower mean temperatures therefore prevent solvent boiling.
381 +
382 + For UA-toluene, the non-bonded potentials between intermolecular sites
383 + have a similar Lennard-Jones formulation. The toluene molecules were
384 + treated as a single rigid body, so there was no need for
385 + intramolecular interactions (including bonds, bends, or torsions) in
386 + this solvent model.
387 +
388 + Besides the TraPPE-UA models, AA models for both organic solvents are
389 + included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
390 + force field is used, and additional explicit hydrogen sites were
391 + included. Besides bonding and non-bonded site-site interactions,
392 + partial charges and the electrostatic interactions were added to each
393 + CT and HC site. For toluene, the United Force Field developed by
394 + Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} was adopted, and
395 + a flexible model for the toluene molecule was utilized which included
396 + bond, bend, torsion, and inversion potentials to enforce ring
397 + planarity.
398 +
399 + The butanethiol capping agent in our simulations, were also modeled
400 + with both UA and AA model. The TraPPE-UA force field includes
401 + parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
402 + UA butanethiol model in our simulations. The OPLS-AA also provides
403 + parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
404 + surfaces do not have the hydrogen atom bonded to sulfur. To derive
405 + suitable parameters for butanethiol adsorbed on Au(111) surfaces, we
406 + adopt the S parameters from Luedtke and Landman\cite{landman:1998} and
407 + modify the parameters for the CTS atom to maintain charge neutrality
408 + in the molecule.  Note that the model choice (UA or AA) for the capping
409 + agent can be different from the solvent. Regardless of model choice,
410 + the force field parameters for interactions between capping agent and
411 + solvent can be derived using Lorentz-Berthelot Mixing Rule:
412 + \begin{eqnarray}
413 +  \sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
414 +  \epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
415 + \end{eqnarray}
416 +
417 + To describe the interactions between metal (Au) and non-metal atoms,
418 + we refer to an adsorption study of alkyl thiols on gold surfaces by
419 + Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
420 + Lennard-Jones form of potential parameters for the interaction between
421 + Au and pseudo-atoms CH$_x$ and S based on a well-established and
422 + widely-used effective potential of Hautman and Klein for the Au(111)
423 + surface.\cite{hautman:4994} As our simulations require the gold slab
424 + to be flexible to accommodate thermal excitation, the pair-wise form
425 + of potentials they developed was used for our study.
426 +
427 + The potentials developed from {\it ab initio} calculations by Leng
428 + {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
429 + interactions between Au and aromatic C/H atoms in toluene. However,
430 + the Lennard-Jones parameters between Au and other types of particles,
431 + (e.g. AA alkanes) have not yet been established. For these
432 + interactions, the Lorentz-Berthelot mixing rule can be used to derive
433 + effective single-atom LJ parameters for the metal using the fit values
434 + for toluene. These are then used to construct reasonable mixing
435 + parameters for the interactions between the gold and other atoms.
436 + Table \ref{MnM} summarizes the ``metal/non-metal'' parameters used in
437 + our simulations.
438 +
439 + \begin{table*}
440 +  \begin{minipage}{\linewidth}
441 +    \begin{center}
442 +      \caption{Non-bonded interaction parameters (including cross
443 +        interactions with Au atoms) for both force fields used in this
444 +        work.}      
445 +      \begin{tabular}{lllllll}
446 +        \hline\hline
447 +        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
448 +        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
449 +        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
450 +        \hline
451 +        United Atom (UA)
452 +        &CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
453 +        &CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
454 +        &CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
455 +        &CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
456 +        \hline
457 +        All Atom (AA)
458 +        &CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
459 +        &CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
460 +        &CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
461 +        &HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
462 +        &CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
463 +        &HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
464 +        \hline
465 +        Both UA and AA
466 +        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
467 +        \hline\hline
468 +      \end{tabular}
469 +      \label{MnM}
470 +    \end{center}
471 +  \end{minipage}
472 + \end{table*}
473 +
474 + \subsection{Vibrational Power Spectrum}
475 +
476 + To investigate the mechanism of interfacial thermal conductance, the
477 + vibrational power spectrum was computed. Power spectra were taken for
478 + individual components in different simulations. To obtain these
479 + spectra, simulations were run after equilibration, in the NVE
480 + ensemble, and without a thermal gradient. Snapshots of configurations
481 + were collected at a frequency that is higher than that of the fastest
482 + vibrations occuring in the simulations. With these configurations, the
483 + velocity auto-correlation functions can be computed:
484 + \begin{equation}
485 + C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
486 + \label{vCorr}
487 + \end{equation}
488 + The power spectrum is constructed via a Fourier transform of the
489 + symmetrized velocity autocorrelation function,
490 + \begin{equation}
491 +  \hat{f}(\omega) =
492 +  \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
493 + \label{fourier}
494 + \end{equation}
495 +
496 + \section{Results and Discussions}
497 + In what follows, how the parameters and protocol of simulations would
498 + affect the measurement of $G$'s is first discussed. With a reliable
499 + protocol and set of parameters, the influence of capping agent
500 + coverage on thermal conductance is investigated. Besides, different
501 + force field models for both solvents and selected deuterated models
502 + were tested and compared. Finally, a summary of the role of capping
503 + agent in the interfacial thermal transport process is given.
504 +
505 + \subsection{How Simulation Parameters Affects $G$}
506 + We have varied our protocol or other parameters of the simulations in
507 + order to investigate how these factors would affect the measurement of
508 + $G$'s. It turned out that while some of these parameters would not
509 + affect the results substantially, some other changes to the
510 + simulations would have a significant impact on the measurement
511 + results.
512 +
513 + In some of our simulations, we allowed $L_x$ and $L_y$ to change
514 + during equilibrating the liquid phase. Due to the stiffness of the
515 + crystalline Au structure, $L_x$ and $L_y$ would not change noticeably
516 + after equilibration. Although $L_z$ could fluctuate $\sim$1\% after a
517 + system is fully equilibrated in the NPT ensemble, this fluctuation, as
518 + well as those of $L_x$ and $L_y$ (which is significantly smaller),
519 + would not be magnified on the calculated $G$'s, as shown in Table
520 + \ref{AuThiolHexaneUA}. This insensivity to $L_i$ fluctuations allows
521 + reliable measurement of $G$'s without the necessity of extremely
522 + cautious equilibration process.
523 +
524 + As stated in our computational details, the spacing filled with
525 + solvent molecules can be chosen within a range. This allows some
526 + change of solvent molecule numbers for the same Au-butanethiol
527 + surfaces. We did this study on our Au-butanethiol/hexane
528 + simulations. Nevertheless, the results obtained from systems of
529 + different $N_{hexane}$ did not indicate that the measurement of $G$ is
530 + susceptible to this parameter. For computational efficiency concern,
531 + smaller system size would be preferable, given that the liquid phase
532 + structure is not affected.
533 +
534 + Our NIVS algorithm allows change of unphysical thermal flux both in
535 + direction and in quantity. This feature extends our investigation of
536 + interfacial thermal conductance. However, the magnitude of this
537 + thermal flux is not arbitary if one aims to obtain a stable and
538 + reliable thermal gradient. A temperature profile would be
539 + substantially affected by noise when $|J_z|$ has a much too low
540 + magnitude; while an excessively large $|J_z|$ that overwhelms the
541 + conductance capacity of the interface would prevent a thermal gradient
542 + to reach a stablized steady state. NIVS has the advantage of allowing
543 + $J$ to vary in a wide range such that the optimal flux range for $G$
544 + measurement can generally be simulated by the algorithm. Within the
545 + optimal range, we were able to study how $G$ would change according to
546 + the thermal flux across the interface. For our simulations, we denote
547 + $J_z$ to be positive when the physical thermal flux is from the liquid
548 + to metal, and negative vice versa. The $G$'s measured under different
549 + $J_z$ is listed in Table \ref{AuThiolHexaneUA} and
550 + \ref{AuThiolToluene}. These results do not suggest that $G$ is
551 + dependent on $J_z$ within this flux range. The linear response of flux
552 + to thermal gradient simplifies our investigations in that we can rely
553 + on $G$ measurement with only a couple $J_z$'s and do not need to test
554 + a large series of fluxes.
555 +
556 + \begin{table*}
557 +  \begin{minipage}{\linewidth}
558 +    \begin{center}
559 +      \caption{Computed interfacial thermal conductivity ($G$ and
560 +        $G^\prime$) values for the 100\% covered Au-butanethiol/hexane
561 +        interfaces with UA model and different hexane molecule numbers
562 +        at different temperatures using a range of energy
563 +        fluxes. Error estimates indicated in parenthesis.}
564 +      
565 +      \begin{tabular}{ccccccc}
566 +        \hline\hline
567 +        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
568 +        $J_z$ & $G$ & $G^\prime$ \\
569 +        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
570 +        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
571 +        \hline
572 +        200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
573 +            & 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
574 +            &     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
575 +            &     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
576 +            &     &     &       &  1.91 & 139(10)   & 101(10)   \\
577 +            &     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
578 +            & 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
579 +            &     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
580 +            &     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
581 +            &     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
582 +        \hline
583 +        250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
584 +            &     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
585 +            & 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
586 +            &     & No  & 0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
587 +            &     &     &       &  1.44 & 76.2(5.0) & 64.8(3.8) \\
588 +            &     &     &       & -0.95 & 56.4(2.5) & 54.4(1.1) \\
589 +            &     &     &       & -1.85 & 47.8(1.1) & 53.5(1.5) \\
590 +        \hline\hline
591 +      \end{tabular}
592 +      \label{AuThiolHexaneUA}
593 +    \end{center}
594 +  \end{minipage}
595 + \end{table*}
596 +
597 + Furthermore, we also attempted to increase system average temperatures
598 + to above 200K. These simulations are first equilibrated in the NPT
599 + ensemble under normal pressure. As stated above, the TraPPE-UA model
600 + for hexane tends to predict a lower boiling point. In our simulations,
601 + hexane had diffculty to remain in liquid phase when NPT equilibration
602 + temperature is higher than 250K. Additionally, the equilibrated liquid
603 + hexane density under 250K becomes lower than experimental value. This
604 + expanded liquid phase leads to lower contact between hexane and
605 + butanethiol as well.[MAY NEED SLAB DENSITY FIGURE]
606 + And this reduced contact would
607 + probably be accountable for a lower interfacial thermal conductance,
608 + as shown in Table \ref{AuThiolHexaneUA}.
609 +
610 + A similar study for TraPPE-UA toluene agrees with the above result as
611 + well. Having a higher boiling point, toluene tends to remain liquid in
612 + our simulations even equilibrated under 300K in NPT
613 + ensembles. Furthermore, the expansion of the toluene liquid phase is
614 + not as significant as that of the hexane. This prevents severe
615 + decrease of liquid-capping agent contact and the results (Table
616 + \ref{AuThiolToluene}) show only a slightly decreased interface
617 + conductance. Therefore, solvent-capping agent contact should play an
618 + important role in the thermal transport process across the interface
619 + in that higher degree of contact could yield increased conductance.
620 +
621 + \begin{table*}
622 +  \begin{minipage}{\linewidth}
623 +    \begin{center}
624 +      \caption{Computed interfacial thermal conductivity ($G$ and
625 +        $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
626 +        interface at different temperatures using a range of energy
627 +        fluxes. Error estimates indicated in parenthesis.}
628 +      
629 +      \begin{tabular}{ccccc}
630 +        \hline\hline
631 +        $\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
632 +        (K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
633 +        \hline
634 +        200 & 0.933 &  2.15 & 204(12) & 113(12) \\
635 +            &       & -1.86 & 180(3)  & 135(21) \\
636 +            &       & -3.93 & 176(5)  & 113(12) \\
637 +        \hline
638 +        300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
639 +            &       & -4.19 & 135(9)  & 113(12) \\
640 +        \hline\hline
641 +      \end{tabular}
642 +      \label{AuThiolToluene}
643 +    \end{center}
644 +  \end{minipage}
645 + \end{table*}
646 +
647 + Besides lower interfacial thermal conductance, surfaces in relatively
648 + high temperatures are susceptible to reconstructions, when
649 + butanethiols have a full coverage on the Au(111) surface. These
650 + reconstructions include surface Au atoms migrated outward to the S
651 + atom layer, and butanethiol molecules embedded into the original
652 + surface Au layer. The driving force for this behavior is the strong
653 + Au-S interactions in our simulations. And these reconstructions lead
654 + to higher ratio of Au-S attraction and thus is energetically
655 + favorable. Furthermore, this phenomenon agrees with experimental
656 + results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. Vlugt
657 + {\it et al.} had kept their Au(111) slab rigid so that their
658 + simulations can reach 300K without surface reconstructions. Without
659 + this practice, simulating 100\% thiol covered interfaces under higher
660 + temperatures could hardly avoid surface reconstructions. However, our
661 + measurement is based on assuming homogeneity on $x$ and $y$ dimensions
662 + so that measurement of $T$ at particular $z$ would be an effective
663 + average of the particles of the same type. Since surface
664 + reconstructions could eliminate the original $x$ and $y$ dimensional
665 + homogeneity, measurement of $G$ is more difficult to conduct under
666 + higher temperatures. Therefore, most of our measurements are
667 + undertaken at $\langle T\rangle\sim$200K.
668 +
669 + However, when the surface is not completely covered by butanethiols,
670 + the simulated system is more resistent to the reconstruction
671 + above. Our Au-butanethiol/toluene system had the Au(111) surfaces 90\%
672 + covered by butanethiols, but did not see this above phenomena even at
673 + $\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
674 + capping agents could help prevent surface reconstruction in that they
675 + provide other means of capping agent relaxation. It is observed that
676 + butanethiols can migrate to their neighbor empty sites during a
677 + simulation. Therefore, we were able to obtain $G$'s for these
678 + interfaces even at a relatively high temperature without being
679 + affected by surface reconstructions.
680 +
681 + \subsection{Influence of Capping Agent Coverage on $G$}
682 + To investigate the influence of butanethiol coverage on interfacial
683 + thermal conductance, a series of different coverage Au-butanethiol
684 + surfaces is prepared and solvated with various organic
685 + molecules. These systems are then equilibrated and their interfacial
686 + thermal conductivity are measured with our NIVS algorithm. Figure
687 + \ref{coverage} demonstrates the trend of conductance change with
688 + respect to different coverages of butanethiol. To study the isotope
689 + effect in interfacial thermal conductance, deuterated UA-hexane is
690 + included as well.
691 +
692 + \begin{figure}
693 + \includegraphics[width=\linewidth]{coverage}
694 + \caption{Comparison of interfacial thermal conductivity ($G$) values
695 +  for the Au-butanethiol/solvent interface with various UA models and
696 +  different capping agent coverages at $\langle T\rangle\sim$200K
697 +  using certain energy flux respectively.}
698 + \label{coverage}
699 + \end{figure}
700 +
701 + It turned out that with partial covered butanethiol on the Au(111)
702 + surface, the derivative definition for $G^\prime$
703 + (Eq. \ref{derivativeG}) was difficult to apply, due to the difficulty
704 + in locating the maximum of change of $\lambda$. Instead, the discrete
705 + definition (Eq. \ref{discreteG}) is easier to apply, as the Gibbs
706 + deviding surface can still be well-defined. Therefore, $G$ (not
707 + $G^\prime$) was used for this section.
708 +
709 + From Figure \ref{coverage}, one can see the significance of the
710 + presence of capping agents. Even when a fraction of the Au(111)
711 + surface sites are covered with butanethiols, the conductivity would
712 + see an enhancement by at least a factor of 3. This indicates the
713 + important role cappping agent is playing for thermal transport
714 + phenomena on metal / organic solvent surfaces.
715 +
716 + Interestingly, as one could observe from our results, the maximum
717 + conductance enhancement (largest $G$) happens while the surfaces are
718 + about 75\% covered with butanethiols. This again indicates that
719 + solvent-capping agent contact has an important role of the thermal
720 + transport process. Slightly lower butanethiol coverage allows small
721 + gaps between butanethiols to form. And these gaps could be filled with
722 + solvent molecules, which acts like ``heat conductors'' on the
723 + surface. The higher degree of interaction between these solvent
724 + molecules and capping agents increases the enhancement effect and thus
725 + produces a higher $G$ than densely packed butanethiol arrays. However,
726 + once this maximum conductance enhancement is reached, $G$ decreases
727 + when butanethiol coverage continues to decrease. Each capping agent
728 + molecule reaches its maximum capacity for thermal
729 + conductance. Therefore, even higher solvent-capping agent contact
730 + would not offset this effect. Eventually, when butanethiol coverage
731 + continues to decrease, solvent-capping agent contact actually
732 + decreases with the disappearing of butanethiol molecules. In this
733 + case, $G$ decrease could not be offset but instead accelerated. [NEED
734 + SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
735 +
736 + A comparison of the results obtained from differenet organic solvents
737 + can also provide useful information of the interfacial thermal
738 + transport process. The deuterated hexane (UA) results do not appear to
739 + be much different from those of normal hexane (UA), given that
740 + butanethiol (UA) is non-deuterated for both solvents. These UA model
741 + studies, even though eliminating C-H vibration samplings, still have
742 + C-C vibrational frequencies different from each other. However, these
743 + differences in the infrared range do not seem to produce an observable
744 + difference for the results of $G$ (Figure \ref{uahxnua}).
745 +
746 + \begin{figure}
747 + \includegraphics[width=\linewidth]{uahxnua}
748 + \caption{Vibrational spectra obtained for normal (upper) and
749 +  deuterated (lower) hexane in Au-butanethiol/hexane
750 +  systems. Butanethiol spectra are shown as reference. Both hexane and
751 +  butanethiol were using United-Atom models.}
752 + \label{uahxnua}
753 + \end{figure}
754 +
755 + Furthermore, results for rigid body toluene solvent, as well as other
756 + UA-hexane solvents, are reasonable within the general experimental
757 + ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
758 + suggests that explicit hydrogen might not be a required factor for
759 + modeling thermal transport phenomena of systems such as
760 + Au-thiol/organic solvent.
761 +
762 + However, results for Au-butanethiol/toluene do not show an identical
763 + trend with those for Au-butanethiol/hexane in that $G$ remains at
764 + approximately the same magnitue when butanethiol coverage differs from
765 + 25\% to 75\%. This might be rooted in the molecule shape difference
766 + for planar toluene and chain-like {\it n}-hexane. Due to this
767 + difference, toluene molecules have more difficulty in occupying
768 + relatively small gaps among capping agents when their coverage is not
769 + too low. Therefore, the solvent-capping agent contact may keep
770 + increasing until the capping agent coverage reaches a relatively low
771 + level. This becomes an offset for decreasing butanethiol molecules on
772 + its effect to the process of interfacial thermal transport. Thus, one
773 + can see a plateau of $G$ vs. butanethiol coverage in our results.
774 +
775 + \subsection{Influence of Chosen Molecule Model on $G$}
776 + In addition to UA solvent/capping agent models, AA models are included
777 + in our simulations as well. Besides simulations of the same (UA or AA)
778 + model for solvent and capping agent, different models can be applied
779 + to different components. Furthermore, regardless of models chosen,
780 + either the solvent or the capping agent can be deuterated, similar to
781 + the previous section. Table \ref{modelTest} summarizes the results of
782 + these studies.
783 +
784 + \begin{table*}
785 +  \begin{minipage}{\linewidth}
786 +    \begin{center}
787 +      
788 +      \caption{Computed interfacial thermal conductivity ($G$ and
789 +        $G^\prime$) values for interfaces using various models for
790 +        solvent and capping agent (or without capping agent) at
791 +        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
792 +        or capping agent molecules; ``Avg.'' denotes results that are
793 +        averages of simulations under different $J_z$'s. Error
794 +        estimates indicated in parenthesis.)}
795 +      
796 +      \begin{tabular}{llccc}
797 +        \hline\hline
798 +        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
799 +        (or bare surface) & model & (GW/m$^2$) &
800 +        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
801 +        \hline
802 +        UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
803 +              & UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
804 +              & AA hexane    & Avg. & 131(6)    & 122(10)   \\
805 +              & UA toluene   & 1.96 & 187(16)   & 151(11)   \\
806 +              & AA toluene   & 1.89 & 200(36)   & 149(53)   \\
807 +        \hline
808 +        AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
809 +              & AA hexane    & Avg. & 442(14)   & 356(31)   \\
810 +              & AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
811 +              & UA toluene   & 1.98 & 125(25)   & 97(60)    \\
812 +              & AA toluene   & 3.79 & 487(56)   & 290(42)   \\
813 +        \hline
814 +        AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
815 +              & AA hexane    & 1.92 & 243(29)   & 191(11)   \\
816 +              & AA toluene   & 1.93 & 364(36)   & 322(67)   \\
817 +        \hline
818 +        bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
819 +              & UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
820 +              & AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
821 +              & UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
822 +        \hline\hline
823 +      \end{tabular}
824 +      \label{modelTest}
825 +    \end{center}
826 +  \end{minipage}
827 + \end{table*}
828 +
829 + To facilitate direct comparison, the same system with differnt models
830 + for different components uses the same length scale for their
831 + simulation cells. Without the presence of capping agent, using
832 + different models for hexane yields similar results for both $G$ and
833 + $G^\prime$, and these two definitions agree with eath other very
834 + well. This indicates very weak interaction between the metal and the
835 + solvent, and is a typical case for acoustic impedance mismatch between
836 + these two phases.
837 +
838 + As for Au(111) surfaces completely covered by butanethiols, the choice
839 + of models for capping agent and solvent could impact the measurement
840 + of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
841 + interfaces, using AA model for both butanethiol and hexane yields
842 + substantially higher conductivity values than using UA model for at
843 + least one component of the solvent and capping agent, which exceeds
844 + the general range of experimental measurement results. This is
845 + probably due to the classically treated C-H vibrations in the AA
846 + model, which should not be appreciably populated at normal
847 + temperatures. In comparison, once either the hexanes or the
848 + butanethiols are deuterated, one can see a significantly lower $G$ and
849 + $G^\prime$. In either of these cases, the C-H(D) vibrational overlap
850 + between the solvent and the capping agent is removed (Figure
851 + \ref{aahxntln}). Conclusively, the improperly treated C-H vibration in
852 + the AA model produced over-predicted results accordingly. Compared to
853 + the AA model, the UA model yields more reasonable results with higher
854 + computational efficiency.
855 +
856 + \begin{figure}
857 + \includegraphics[width=\linewidth]{aahxntln}
858 + \caption{Spectra obtained for All-Atom model Au-butanethil/solvent
859 +  systems. When butanethiol is deuterated (lower left), its
860 +  vibrational overlap with hexane would decrease significantly,
861 +  compared with normal butanethiol (upper left). However, this
862 +  dramatic change does not apply to toluene as much (right).}
863 + \label{aahxntln}
864 + \end{figure}
865 +
866 + However, for Au-butanethiol/toluene interfaces, having the AA
867 + butanethiol deuterated did not yield a significant change in the
868 + measurement results. Compared to the C-H vibrational overlap between
869 + hexane and butanethiol, both of which have alkyl chains, that overlap
870 + between toluene and butanethiol is not so significant and thus does
871 + not have as much contribution to the heat exchange
872 + process. Conversely, extra degrees of freedom such as the C-H
873 + vibrations could yield higher heat exchange rate between these two
874 + phases and result in a much higher conductivity.
875 +
876 + Although the QSC model for Au is known to predict an overly low value
877 + for bulk metal gold conductivity\cite{kuang:164101}, our computational
878 + results for $G$ and $G^\prime$ do not seem to be affected by this
879 + drawback of the model for metal. Instead, our results suggest that the
880 + modeling of interfacial thermal transport behavior relies mainly on
881 + the accuracy of the interaction descriptions between components
882 + occupying the interfaces.
883 +
884 + \subsection{Role of Capping Agent in Interfacial Thermal Conductance}
885 + The vibrational spectra for gold slabs in different environments are
886 + shown as in Figure \ref{specAu}. Regardless of the presence of
887 + solvent, the gold surfaces covered by butanethiol molecules, compared
888 + to bare gold surfaces, exhibit an additional peak observed at the
889 + frequency of $\sim$170cm$^{-1}$, which is attributed to the S-Au
890 + bonding vibration. This vibration enables efficient thermal transport
891 + from surface Au layer to the capping agents. Therefore, in our
892 + simulations, the Au/S interfaces do not appear major heat barriers
893 + compared to the butanethiol / solvent interfaces.
894 +
895 + Simultaneously, the vibrational overlap between butanethiol and
896 + organic solvents suggests higher thermal exchange efficiency between
897 + these two components. Even exessively high heat transport was observed
898 + when All-Atom models were used and C-H vibrations were treated
899 + classically. Compared to metal and organic liquid phase, the heat
900 + transfer efficiency between butanethiol and organic solvents is closer
901 + to that within bulk liquid phase.
902 +
903 + Furthermore, our observation validated previous
904 + results\cite{hase:2010} that the intramolecular heat transport of
905 + alkylthiols is highly effecient. As a combinational effects of these
906 + phenomena, butanethiol acts as a channel to expedite thermal transport
907 + process. The acoustic impedance mismatch between the metal and the
908 + liquid phase can be effectively reduced with the presence of suitable
909 + capping agents.
910 +
911 + \begin{figure}
912 + \includegraphics[width=\linewidth]{vibration}
913 + \caption{Vibrational spectra obtained for gold in different
914 +  environments.}
915 + \label{specAu}
916 + \end{figure}
917 +
918 + [MAY ADD COMPARISON OF AU SLAB WIDTHS]
919 +
920 + \section{Conclusions}
921 + The NIVS algorithm we developed has been applied to simulations of
922 + Au-butanethiol surfaces with organic solvents. This algorithm allows
923 + effective unphysical thermal flux transferred between the metal and
924 + the liquid phase. With the flux applied, we were able to measure the
925 + corresponding thermal gradient and to obtain interfacial thermal
926 + conductivities. Under steady states, single trajectory simulation
927 + would be enough for accurate measurement. This would be advantageous
928 + compared to transient state simulations, which need multiple
929 + trajectories to produce reliable average results.
930 +
931 + Our simulations have seen significant conductance enhancement with the
932 + presence of capping agent, compared to the bare gold / liquid
933 + interfaces. The acoustic impedance mismatch between the metal and the
934 + liquid phase is effectively eliminated by proper capping
935 + agent. Furthermore, the coverage precentage of the capping agent plays
936 + an important role in the interfacial thermal transport
937 + process. Moderately lower coverages allow higher contact between
938 + capping agent and solvent, and thus could further enhance the heat
939 + transfer process.
940 +
941 + Our measurement results, particularly of the UA models, agree with
942 + available experimental data. This indicates that our force field
943 + parameters have a nice description of the interactions between the
944 + particles at the interfaces. AA models tend to overestimate the
945 + interfacial thermal conductance in that the classically treated C-H
946 + vibration would be overly sampled. Compared to the AA models, the UA
947 + models have higher computational efficiency with satisfactory
948 + accuracy, and thus are preferable in interfacial thermal transport
949 + modelings. Of the two definitions for $G$, the discrete form
950 + (Eq. \ref{discreteG}) was easier to use and gives out relatively
951 + consistent results, while the derivative form (Eq. \ref{derivativeG})
952 + is not as versatile. Although $G^\prime$ gives out comparable results
953 + and follows similar trend with $G$ when measuring close to fully
954 + covered or bare surfaces, the spatial resolution of $T$ profile is
955 + limited for accurate computation of derivatives data.
956 +
957 + Vlugt {\it et al.} has investigated the surface thiol structures for
958 + nanocrystal gold and pointed out that they differs from those of the
959 + Au(111) surface\cite{vlugt:cpc2007154}. This difference might lead to
960 + change of interfacial thermal transport behavior as well. To
961 + investigate this problem, an effective means to introduce thermal flux
962 + and measure the corresponding thermal gradient is desirable for
963 + simulating structures with spherical symmetry.
964 +
965   \section{Acknowledgments}
966   Support for this project was provided by the National Science
967   Foundation under grant CHE-0848243. Computational time was provided by
968   the Center for Research Computing (CRC) at the University of Notre
969 < Dame.  \newpage
969 > Dame. \newpage
970  
971   \bibliography{interfacial}
972  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines