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 3729 by skuang, Wed Jun 29 18:14:23 2011 UTC vs.
Revision 3755 by skuang, Fri Jul 29 15:45:14 2011 UTC

# Line 23 | Line 23
23   \setlength{\belowcaptionskip}{30 pt}
24  
25   %\renewcommand\citemid{\ } % no comma in optional reference note
26 < \bibpunct{[}{]}{,}{s}{}{;}
27 < \bibliographystyle{aip}
26 > \bibpunct{[}{]}{,}{n}{}{;}
27 > \bibliographystyle{achemso}
28  
29   \begin{document}
30  
# Line 45 | Line 45 | We have developed a Non-Isotropic Velocity Scaling alg
45  
46   \begin{abstract}
47  
48 < We have developed a Non-Isotropic Velocity Scaling algorithm for
49 < setting up and maintaining stable thermal gradients in non-equilibrium
50 < molecular dynamics simulations. This approach effectively imposes
51 < unphysical thermal flux even between particles of different
52 < identities, conserves linear momentum and kinetic energy, and
53 < minimally perturbs the velocity profile of a system when compared with
54 < previous RNEMD methods. We have used this method to simulate thermal
55 < conductance at metal / organic solvent interfaces both with and
56 < without the presence of thiol-based capping agents.  We obtained
57 < values comparable with experimental values, and observed significant
58 < conductance enhancement with the presence of capping agents. Computed
59 < power spectra indicate the acoustic impedance mismatch between metal
60 < and liquid phase is greatly reduced by the capping agents and thus
61 < leads to higher interfacial thermal transfer efficiency.
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  
# Line 71 | Line 73 | leads to higher interfacial thermal transfer efficienc
73   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74  
75   \section{Introduction}
76 < [BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM]
77 < Interfacial thermal conductance is extensively studied both
78 < experimentally and computationally, and systems with interfaces
79 < present are generally heterogeneous. Although interfaces are commonly
80 < barriers to heat transfer, it has been
81 < reported\cite{doi:10.1021/la904855s} that under specific circustances,
82 < e.g. with certain capping agents present on the surface, interfacial
83 < conductance can be significantly enhanced. However, heat conductance
82 < of molecular and nano-scale interfaces will be affected by the
83 < chemical details of the surface and is challenging to
84 < experimentalist. The lower thermal flux through interfaces is even
85 < more difficult to measure with EMD and forward NEMD simulation
86 < methods. Therefore, developing good simulation methods will be
87 < desirable in order to investigate thermal transport across interfaces.
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 < Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
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 > 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
106 > this phenomena at the molecular level. Recently, Hase and coworkers
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
111 > measurements for heat conductance of interfaces between the capping
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 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
130   conditions) while establishing true thermal distributions in each of
131 < the two slabs. Furthermore, it allows more effective thermal exchange
132 < between particles of different identities, and thus enables extensive
133 < study of interfacial conductance.
131 > the two slabs. Furthermore, it allows effective thermal exchange
132 > between particles of different identities, and thus makes the study of
133 > interfacial conductance much simpler.
134  
135 + The work presented here deals with the Au(111) surface covered to
136 + varying degrees by butanethiol, a capping agent with short carbon
137 + chain, and solvated with organic solvents of different molecular
138 + properties. Different models were used for both the capping agent and
139 + the solvent force field parameters. Using the NIVS algorithm, the
140 + thermal transport across these interfaces was studied and the
141 + underlying mechanism for the phenomena was investigated.
142 +
143   \section{Methodology}
144 < \subsection{Algorithm}
145 < [BACKGROUND FOR MD METHODS]
146 < There have been many algorithms for computing thermal conductivity
147 < using molecular dynamics simulations. However, interfacial conductance
148 < is at least an order of magnitude smaller. This would make the
149 < calculation even more difficult for those slowly-converging
150 < equilibrium methods. Imposed-flux non-equilibrium
151 < methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
152 < the response of temperature or momentum gradients are easier to
153 < measure than the flux, if unknown, and thus, is a preferable way to
154 < the forward NEMD methods. Although the momentum swapping approach for
155 < flux-imposing can be used for exchanging energy between particles of
156 < different identity, the kinetic energy transfer efficiency is affected
157 < by the mass difference between the particles, which limits its
158 < application on heterogeneous interfacial systems.
144 > \subsection{Imposd-Flux Methods in MD Simulations}
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 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 in
164 < non-equilibrium MD simulations is able to impose relatively large
165 < kinetic energy flux without obvious perturbation to the velocity
166 < distribution of the simulated systems. Furthermore, this approach has
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
168   can be applied between regions of particles of arbitary identity, and
169 < the flux quantity is not restricted by particle mass difference.
169 > the flux will not be restricted by difference in particle mass.
170  
171   The NIVS algorithm scales the velocity vectors in two separate regions
172   of a simulation system with respective diagonal scaling matricies. To
173   determine these scaling factors in the matricies, a set of equations
174   including linear momentum conservation and kinetic energy conservation
175 < constraints and target momentum/energy flux satisfaction is
176 < solved. With the scaling operation applied to the system in a set
177 < frequency, corresponding momentum/temperature gradients can be built,
178 < which can be used for computing transportation properties and other
179 < applications related to momentum/temperature gradients. The NIVS
132 < algorithm conserves momenta and energy and does not depend on an
133 < external thermostat.
175 > constraints and target energy flux satisfaction is solved. With the
176 > scaling operation applied to the system in a set frequency, bulk
177 > temperature gradients can be easily established, and these can be used
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 < For interfaces with a relatively low interfacial conductance, the bulk
183 < regions on either side of an interface rapidly come to a state in
184 < which the two phases have relatively homogeneous (but distinct)
185 < temperatures. The interfacial thermal conductivity $G$ can therefore
186 < 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, two ways can be
199 < used to define $G$.
200 <
201 < One way is to assume the temperature is discretely different on two
202 < sides of the interface, $G$ can be calculated with the thermal flux
203 < applied $J$ and the maximum temperature difference measured along the
204 < thermal gradient max($\Delta T$), which occurs at the interface, as:
198 > When the interfacial conductance is {\it not} small, there are two
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  
211 + \begin{figure}
212 + \includegraphics[width=\linewidth]{method}
213 + \caption{Interfacial conductance can be calculated by applying an
214 +  (unphysical) kinetic energy flux between two slabs, one located
215 +  within the metal and another on the edge of the periodic box.  The
216 +  system responds by forming a thermal response or a gradient.  In
217 +  bulk liquids, this gradient typically has a single slope, but in
218 +  interfacial systems, there are distinct thermal conductivity
219 +  domains.  The interfacial conductance, $G$ is found by measuring the
220 +  temperature gap at the Gibbs dividing surface, or by using second
221 +  derivatives of the thermal profile.}
222 + \label{demoPic}
223 + \end{figure}
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 173 | 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 method 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 are used for calculation and comparison.
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 < [IMPOSE G DEFINITION INTO OUR SYSTEMS]
251 < To facilitate the use of the above definitions in calculating $G$ and
252 < $G^\prime$, we have a metal slab with its (111) surfaces perpendicular
253 < to the $z$-axis of our simulation cells. With or withour capping
254 < agents on the surfaces, the metal slab is solvated with organic
255 < solvents, as illustrated in Figure \ref{demoPic}.
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 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}
190 – \includegraphics[width=\linewidth]{demoPic}
191 – \caption{A sample showing how a metal slab has its (111) surface
192 –  covered by capping agent molecules and solvated by hexane.}
193 – \label{demoPic}
194 – \end{figure}
195 –
196 – With a simulation cell setup following the above manner, one is able
197 – to equilibrate the system and impose an unphysical thermal flux
198 – between the liquid and the metal phase with the NIVS algorithm. Under
199 – a stablized thermal gradient induced by periodically applying the
200 – unphysical flux, one is able to obtain a temperature profile and the
201 – physical thermal flux corresponding to it, which equals to the
202 – unphysical flux applied by NIVS. These data enables the evaluation of
203 – the interfacial thermal conductance of a surface. Figure \ref{gradT}
204 – is an example how those stablized thermal gradient can be used to
205 – obtain the 1st and 2nd derivatives of the temperature profile.
206 –
207 – \begin{figure}
259   \includegraphics[width=\linewidth]{gradT}
260 < \caption{The 1st and 2nd derivatives of temperature profile can be
261 <  obtained with finite difference approximation.}
260 > \caption{A sample of Au-butanethiol/hexane interfacial system and the
261 >  temperature profile after a kinetic energy flux is imposed to
262 >  it. The 1st and 2nd derivatives of the temperature profile can be
263 >  obtained with finite difference approximation (lower panel).}
264   \label{gradT}
265   \end{figure}
266  
267   \section{Computational Details}
268 < \subsection{System Geometry}
269 < In our simulations, Au is used to construct a metal slab with bare
270 < (111) surface perpendicular to the $z$-axis. Different slab thickness
271 < (layer numbers of Au) are simulated. This metal slab is first
272 < equilibrated under normal pressure (1 atm) and a desired
273 < temperature. After equilibration, butanethiol is used as the capping
274 < agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
275 < atoms in the butanethiol molecules would occupy the three-fold sites
276 < of the surfaces, and the maximal butanethiol capacity on Au surface is
277 < $1/3$ of the total number of surface Au atoms[CITATION]. A series of
278 < different coverage surfaces is investigated in order to study the
279 < relation between coverage and conductance.
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 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 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 < [COVERAGE DISCRIPTION] However, since the interactions between surface
287 < Au and butanethiol is non-bonded, the capping agent molecules are
288 < allowed to migrate to an empty neighbor three-fold site during a
289 < simulation. Therefore, the initial configuration would not severely
290 < affect the sampling of a variety of configurations of the same
291 < coverage, and the final conductance measurement would be an average
292 < effect of these configurations explored in the simulations. [MAY NEED FIGURES]
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 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.
293  
294 < After the modified Au-butanethiol surface systems are equilibrated
295 < under canonical ensemble, Packmol\cite{packmol} is used to pack
296 < organic solvent molecules in the previously vacuum part of the
297 < simulation cells, which guarantees that short range repulsive
298 < interactions do not disrupt the simulations. Two solvents are
299 < investigated, one which has little vibrational overlap with the
300 < alkanethiol and plane-like shape (toluene), and one which has similar
243 < vibrational frequencies and chain-like shape ({\it n}-hexane). The
244 < spacing filled by solvent molecules, i.e. the gap between periodically
245 < repeated Au-butanethiol surfaces should be carefully chosen so that it
246 < would not be too short to affect the liquid phase structure, nor too
247 < long, leading to over cooling (freezing) or heating (boiling) when a
248 < thermal flux is applied. In our simulations, this spacing is usually
249 < $35 \sim 60$\AA.
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 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 initial configurations generated by Packmol are further
303 < equilibrated with the $x$ and $y$ dimensions fixed, only allowing
304 < length scale change in $z$ dimension. This is to ensure that the
305 < equilibration of liquid phase does not affect the metal crystal
306 < structure in $x$ and $y$ dimensions. Further equilibration are run
307 < under NVT and then NVE ensembles.
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. The spacing
310 > between periodic images of the gold interfaces is $45 \sim 75$\AA.
311  
312 < After the systems reach equilibrium, NIVS is implemented to impose a
313 < periodic unphysical thermal flux between the metal and the liquid
314 < phase. Most of our simulations are under an average temperature of
315 < $\sim$200K. Therefore, this flux usually comes from the metal to the
312 > The initial configurations generated are further equilibrated with the
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 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. This induced temperature
328 < gradient is stablized and the simulation cell is devided evenly into
329 < N slabs along the $z$-axis and the temperatures of each slab are
330 < recorded. When the slab width $d$ of each slab is the same, the
331 < derivatives of $T$ with respect to slab number $n$ can be directly
332 < used for $G^\prime$ calculations:
333 < \begin{equation}
334 < G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
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: \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 276 | Line 340 | G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}
340   \label{derivativeG2}
341   \end{equation}
342  
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. Therefore, force field
349 < parameter descriptions are needed for interactions both between the
350 < same type of particles and between particles of different species.
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  
355 + \begin{figure}
356 + \includegraphics[width=\linewidth]{structures}
357 + \caption{Structures of the capping agent and solvents utilized in
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,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
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, pseudo-atoms are
378 < located at the carbon centers for alkyl groups. By eliminating
379 < explicit hydrogen atoms, these models are simple and computationally
380 < efficient, while maintains good accuracy. However, the TraPPE-UA for
381 < alkanes is known to predict a lower boiling point than experimental
299 < values. Considering that after an unphysical thermal flux is applied
300 < to a system, the temperature of ``hot'' area in the liquid phase would be
301 < significantly higher than the average, to prevent over heating and
302 < boiling of the liquid phase, the average temperature in our
303 < simulations should be much lower than the liquid boiling point. [NEED MORE DISCUSSION]
304 < For UA-toluene model, rigid body constraints are applied, so that the
305 < benzene ring and the methyl-C(aromatic) bond are kept rigid. This
306 < would save computational time.[MORE DETAILS NEEDED]
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 closer than 3 bonds. For non-bonded interactions, Lennard-Jones
381 > potentials are used.
382  
383 < Besides the TraPPE-UA models, AA models for both organic solvents are
384 < included in our studies as well. For hexane, the OPLS
385 < all-atom\cite{OPLSAA} force field is used. [MORE DETAILS]
386 < For toluene, the United Force Field developed by Rapp\'{e} {\it et
387 <  al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]
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 < The capping agent in our simulations, the butanethiol molecules can
394 < either use UA or AA model. The TraPPE-UA force fields includes
395 < parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used in
396 < our simulations corresponding to our TraPPE-UA models for solvent.
397 < and All-Atom models [NEED CITATIONS]
319 < However, the model choice (UA or AA) of capping agent can be different
320 < from the solvent. Regardless of model choice, the force field
321 < parameters for interactions between capping agent and solvent can be
322 < derived using Lorentz-Berthelot Mixing Rule.
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 < To describe the interactions between metal Au and non-metal capping
400 < agent and solvent, we refer to Vlugt\cite{vlugt:cpc2007154} and derive
401 < other interactions which are not yet finely parametrized. [can add
402 < hautman and klein's paper here and more discussion; need to put
403 < aromatic-metal interaction approximation here]\cite{doi:10.1021/jp034405s}
399 > Besides the TraPPE-UA models, AA models for both organic solvents are
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, 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 < [TABULATED FORCE FIELD PARAMETERS NEEDED]
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 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}}
424 > \end{eqnarray}
425  
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 < [SURFACE RECONSTRUCTION PREVENTS SIMULATION TEMP TO GO HIGHER]
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 >
448 > \begin{table*}
449 >  \begin{minipage}{\linewidth}
450 >    \begin{center}
451 >      \caption{Non-bonded interaction parameters (including cross
452 >        interactions with Au atoms) for both force fields used in this
453 >        work.}      
454 >      \begin{tabular}{lllllll}
455 >        \hline\hline
456 >        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
457 >        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
458 >        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
459 >        \hline
460 >        United Atom (UA)
461 >        &CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
462 >        &CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
463 >        &CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
464 >        &CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
465 >        \hline
466 >        All Atom (AA)
467 >        &CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
468 >        &CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
469 >        &CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
470 >        &HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
471 >        &CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
472 >        &HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
473 >        \hline
474 >        Both UA and AA
475 >        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
476 >        \hline\hline
477 >      \end{tabular}
478 >      \label{MnM}
479 >    \end{center}
480 >  \end{minipage}
481 > \end{table*}
482  
483  
484   \section{Results}
485 < [REARRANGEMENT NEEDED]
486 < \subsection{Toluene Solvent}
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 < The results (Table \ref{AuThiolToluene}) show a
341 < significant conductance enhancement compared to the gold/water
342 < interface without capping agent and agree with available experimental
343 < data. This indicates that the metal-metal potential, though not
344 < predicting an accurate bulk metal thermal conductivity, does not
345 < greatly interfere with the simulation of the thermal conductance
346 < behavior across a non-metal interface. The solvent model is not
347 < particularly volatile, so the simulation cell does not expand
348 < significantly under higher temperature. We did not observe a
349 < significant conductance decrease when the temperature was increased to
350 < 300K. The results show that the two definitions used for $G$ yield
351 < comparable values, though $G^\prime$ tends to be smaller.
493 > \subsection{Effects due to capping agent coverage}
494  
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}
356 –      \caption{Computed interfacial thermal conductivity ($G$ and
357 –        $G^\prime$) values for the Au/butanethiol/toluene interface at
358 –        different temperatures using a range of energy fluxes.}
528        
529 <      \begin{tabular}{cccc}
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 <        $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
540 <        (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
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 <        200 & 1.86 & 180 & 135 \\
544 <            & 2.15 & 204 & 113 \\
545 <            & 3.93 & 175 & 114 \\
546 <        300 & 1.91 & 143 & 125 \\
547 <            & 4.19 & 134 & 113 \\
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{AuThiolToluene}
565 >      \label{modelTest}
566      \end{center}
567    \end{minipage}
568   \end{table*}
569  
570 < \subsection{Hexane Solvent}
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 < Using the united-atom model, different coverages of capping agent,
380 < temperatures of simulations and numbers of solvent molecules were all
381 < investigated and Table \ref{AuThiolHexaneUA} shows the results of
382 < these computations. The number of hexane molecules in our simulations
383 < does not affect the calculations significantly. However, a very long
384 < length scale for the thermal gradient axis ($z$) may cause excessively
385 < hot or cold temperatures in the middle of the solvent region and lead
386 < to undesired phenomena such as solvent boiling or freezing, while too
387 < few solvent molecules would change the normal behavior of the liquid
388 < phase. Our $N_{hexane}$ values were chosen to ensure that these
389 < extreme cases did not happen to our simulations.
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 < Table \ref{AuThiolHexaneUA} enables direct comparison between
392 < different coverages of capping agent, when other system parameters are
393 < held constant. With high coverage of butanethiol on the gold surface,
394 < the interfacial thermal conductance is enhanced
395 < significantly. Interestingly, a slightly lower butanethiol coverage
396 < leads to a moderately higher conductivity. This is probably due to
397 < more solvent/capping agent contact when butanethiol molecules are
398 < not densely packed, which enhances the interactions between the two
399 < phases and lowers the thermal transfer barrier of this interface.
400 < % [COMPARE TO AU/WATER IN PAPER]
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 < It is also noted that the overall simulation temperature is another
403 < factor that affects the interfacial thermal conductance. One
404 < possibility of this effect may be rooted in the decrease in density of
405 < the liquid phase. We observed that when the average temperature
406 < increases from 200K to 250K, the bulk hexane density becomes lower
407 < than experimental value, as the system is equilibrated under NPT
408 < ensemble. This leads to lower contact between solvent and capping
409 < agent, and thus lower conductivity.
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 < Conductivity values are more difficult to obtain under higher
579 < temperatures. This is because the Au surface tends to undergo
580 < reconstructions in relatively high temperatures. Surface Au atoms can
581 < migrate outward to reach higher Au-S contact; and capping agent
582 < molecules can be embedded into the surface Au layer due to the same
583 < driving force. This phenomenon agrees with experimental
584 < results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. A surface
418 < fully covered in capping agent is more susceptible to reconstruction,
419 < possibly because fully coverage prevents other means of capping agent
420 < relaxation, such as migration to an empty neighbor three-fold site.
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 < %MAY ADD MORE DATA TO TABLE
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
602 > surfaces. We did this study on our Au-butanethiol/hexane
603 > simulations. Nevertheless, the results obtained from systems of
604 > different $N_{hexane}$ did not indicate that the measurement of $G$ is
605 > susceptible to this parameter. For computational efficiency concern,
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
613 > thermal flux is not arbitary if one aims to obtain a stable and
614 > reliable thermal gradient. A temperature profile would be
615 > substantially affected by noise when $|J_z|$ has a much too low
616 > magnitude; while an excessively large $|J_z|$ that overwhelms the
617 > conductance capacity of the interface would prevent a thermal gradient
618 > to reach a stablized steady state. NIVS has the advantage of allowing
619 > $J$ to vary in a wide range such that the optimal flux range for $G$
620 > measurement can generally be simulated by the algorithm. Within the
621 > optimal range, we were able to study how $G$ would change according to
622 > the thermal flux across the interface. For our simulations, we denote
623 > $J_z$ to be positive when the physical thermal flux is from the liquid
624 > to metal, and negative vice versa. The $G$'s measured under different
625 > $J_z$ is listed in Table \ref{AuThiolHexaneUA} and
626 > \ref{AuThiolToluene}. These results do not suggest that $G$ is
627 > dependent on $J_z$ within this flux range. The linear response of flux
628 > to thermal gradient simplifies our investigations in that we can rely
629 > on $G$ measurement with only a couple $J_z$'s and do not need to test
630 > a large series of fluxes.
631 >
632   \begin{table*}
633    \begin{minipage}{\linewidth}
634      \begin{center}
635        \caption{Computed interfacial thermal conductivity ($G$ and
636 <        $G^\prime$) values for the Au/butanethiol/hexane interface
637 <        with united-atom model and different capping agent coverage
638 <        and solvent molecule numbers at different temperatures using a
639 <        range of energy fluxes.}
636 >        $G^\prime$) values for the 100\% covered Au-butanethiol/hexane
637 >        interfaces with UA model and different hexane molecule numbers
638 >        at different temperatures using a range of energy
639 >        fluxes. Error estimates indicated in parenthesis.}
640        
641 <      \begin{tabular}{cccccc}
641 >      \begin{tabular}{ccccccc}
642          \hline\hline
643 <        Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\
644 <        coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) &
643 >        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
644 >        $J_z$ & $G$ & $G^\prime$ \\
645 >        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
646          \multicolumn{2}{c}{(MW/m$^2$/K)} \\
647          \hline
648 <        0.0   & 200 & 200 & 0.96 & 43.3 & 42.7 \\
649 <              &     &     & 1.91 & 45.7 & 42.9 \\
650 <              &     & 166 & 0.96 & 43.1 & 53.4 \\
651 <        88.9  & 200 & 166 & 1.94 & 172  & 108  \\
652 <        100.0 & 250 & 200 & 0.96 & 81.8 & 67.0 \\
653 <              &     & 166 & 0.98 & 79.0 & 62.9 \\
654 <              &     &     & 1.44 & 76.2 & 64.8 \\
655 <              & 200 & 200 & 1.92 & 129  & 87.3 \\
656 <              &     &     & 1.93 & 131  & 77.5 \\
657 <              &     & 166 & 0.97 & 115  & 69.3 \\
658 <              &     &     & 1.94 & 125  & 87.1 \\
648 >        200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
649 >            & 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
650 >            &     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
651 >            &     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
652 >            &     &     &       &  1.91 & 139(10)   & 101(10)   \\
653 >            &     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
654 >            & 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
655 >            &     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
656 >            &     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
657 >            &     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
658 >        \hline
659 >        250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
660 >            &     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
661 >            & 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
662 >            &     & No  & 0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
663 >            &     &     &       &  1.44 & 76.2(5.0) & 64.8(3.8) \\
664 >            &     &     &       & -0.95 & 56.4(2.5) & 54.4(1.1) \\
665 >            &     &     &       & -1.85 & 47.8(1.1) & 53.5(1.5) \\
666          \hline\hline
667        \end{tabular}
668        \label{AuThiolHexaneUA}
# Line 453 | Line 670 | For the all-atom model, the liquid hexane phase was no
670    \end{minipage}
671   \end{table*}
672  
673 < For the all-atom model, the liquid hexane phase was not stable under NPT
457 < conditions. Therefore, the simulation length scale parameters are
458 < adopted from previous equilibration results of the united-atom model
459 < at 200K. Table \ref{AuThiolHexaneAA} shows the results of these
460 < simulations. The conductivity values calculated with full capping
461 < agent coverage are substantially larger than observed in the
462 < united-atom model, and is even higher than predicted by
463 < experiments. It is possible that our parameters for metal-non-metal
464 < particle interactions lead to an overestimate of the interfacial
465 < thermal conductivity, although the active C-H vibrations in the
466 < all-atom model (which should not be appreciably populated at normal
467 < temperatures) could also account for this high conductivity. The major
468 < thermal transfer barrier of Au/butanethiol/hexane interface is between
469 < the liquid phase and the capping agent, so extra degrees of freedom
470 < such as the C-H vibrations could enhance heat exchange between these
471 < two phases and result in a much higher conductivity.
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
678 + for hexane tends to predict a lower boiling point. In our simulations,
679 + hexane had diffculty to remain in liquid phase when NPT equilibration
680 + temperature is higher than 250K. Additionally, the equilibrated liquid
681 + hexane density under 250K becomes lower than experimental value. This
682 + expanded liquid phase leads to lower contact between hexane and
683 + butanethiol as well.[MAY NEED SLAB DENSITY FIGURE]
684 + And this reduced contact would
685 + probably be accountable for a lower interfacial thermal conductance,
686 + as shown in Table \ref{AuThiolHexaneUA}.
687 +
688 + A similar study for TraPPE-UA toluene agrees with the above result as
689 + well. Having a higher boiling point, toluene tends to remain liquid in
690 + our simulations even equilibrated under 300K in NPT
691 + ensembles. Furthermore, the expansion of the toluene liquid phase is
692 + not as significant as that of the hexane. This prevents severe
693 + decrease of liquid-capping agent contact and the results (Table
694 + \ref{AuThiolToluene}) show only a slightly decreased interface
695 + conductance. Therefore, solvent-capping agent contact should play an
696 + important role in the thermal transport process across the interface
697 + in that higher degree of contact could yield increased conductance.
698 +
699   \begin{table*}
700    \begin{minipage}{\linewidth}
701      \begin{center}
476 –      
702        \caption{Computed interfacial thermal conductivity ($G$ and
703 <        $G^\prime$) values for the Au/butanethiol/hexane interface
704 <        with all-atom model and different capping agent coverage at
705 <        200K using a range of energy fluxes.}
703 >        $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
704 >        interface at different temperatures using a range of energy
705 >        fluxes. Error estimates indicated in parenthesis.}
706        
707 <      \begin{tabular}{cccc}
707 >      \begin{tabular}{ccccc}
708          \hline\hline
709 <        Thiol & $J_z$ & $G$ & $G^\prime$ \\
710 <        coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
709 >        $\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
710 >        (K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
711          \hline
712 <        0.0   & 0.95 & 28.5 & 27.2 \\
713 <              & 1.88 & 30.3 & 28.9 \\
714 <        100.0 & 2.87 & 551  & 294  \\
715 <              & 3.81 & 494  & 193  \\
712 >        200 & 0.933 &  2.15 & 204(12) & 113(12) \\
713 >            &       & -1.86 & 180(3)  & 135(21) \\
714 >            &       & -3.93 & 176(5)  & 113(12) \\
715 >        \hline
716 >        300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
717 >            &       & -4.19 & 135(9)  & 113(12) \\
718          \hline\hline
719        \end{tabular}
720 <      \label{AuThiolHexaneAA}
720 >      \label{AuThiolToluene}
721      \end{center}
722    \end{minipage}
723   \end{table*}
724  
725 < %subsubsection{Vibrational spectrum study on conductance mechanism}
726 < To investigate the mechanism of this interfacial thermal conductance,
727 < the vibrational spectra of various gold systems were obtained and are
728 < shown as in the upper panel of Fig. \ref{vibration}. To obtain these
729 < spectra, one first runs a simulation in the NVE ensemble and collects
730 < snapshots of configurations; these configurations are used to compute
731 < the velocity auto-correlation functions, which is used to construct a
732 < power spectrum via a Fourier transform. The gold surfaces covered by
733 < butanethiol molecules exhibit an additional peak observed at a
734 < frequency of $\sim$170cm$^{-1}$, which is attributed to the vibration
735 < of the S-Au bond. This vibration enables efficient thermal transport
736 < from surface Au atoms to the capping agents. Simultaneously, as shown
737 < in the lower panel of Fig. \ref{vibration}, the large overlap of the
738 < vibration spectra of butanethiol and hexane in the all-atom model,
739 < including the C-H vibration, also suggests high thermal exchange
740 < efficiency. The combination of these two effects produces the drastic
741 < interfacial thermal conductance enhancement in the all-atom model.
725 > Besides lower interfacial thermal conductance, surfaces in relatively
726 > high temperatures are susceptible to reconstructions, when
727 > butanethiols have a full coverage on the Au(111) surface. These
728 > reconstructions include surface Au atoms migrated outward to the S
729 > atom layer, and butanethiol molecules embedded into the original
730 > surface Au layer. The driving force for this behavior is the strong
731 > Au-S interactions in our simulations. And these reconstructions lead
732 > to higher ratio of Au-S attraction and thus is energetically
733 > favorable. Furthermore, this phenomenon agrees with experimental
734 > results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. Vlugt
735 > {\it et al.} had kept their Au(111) slab rigid so that their
736 > simulations can reach 300K without surface reconstructions. Without
737 > this practice, simulating 100\% thiol covered interfaces under higher
738 > temperatures could hardly avoid surface reconstructions. However, our
739 > measurement is based on assuming homogeneity on $x$ and $y$ dimensions
740 > so that measurement of $T$ at particular $z$ would be an effective
741 > average of the particles of the same type. Since surface
742 > reconstructions could eliminate the original $x$ and $y$ dimensional
743 > homogeneity, measurement of $G$ is more difficult to conduct under
744 > higher temperatures. Therefore, most of our measurements are
745 > undertaken at $\langle T\rangle\sim$200K.
746  
747 + However, when the surface is not completely covered by butanethiols,
748 + the simulated system is more resistent to the reconstruction
749 + above. Our Au-butanethiol/toluene system had the Au(111) surfaces 90\%
750 + covered by butanethiols, but did not see this above phenomena even at
751 + $\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
752 + capping agents could help prevent surface reconstruction in that they
753 + provide other means of capping agent relaxation. It is observed that
754 + butanethiols can migrate to their neighbor empty sites during a
755 + simulation. Therefore, we were able to obtain $G$'s for these
756 + interfaces even at a relatively high temperature without being
757 + affected by surface reconstructions.
758 +
759 +
760 + \section{Discussion}
761 +
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 +
784 +
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
789 + to bare gold surfaces, exhibit an additional peak observed at the
790 + frequency of $\sim$170cm$^{-1}$, which is attributed to the S-Au
791 + bonding vibration. This vibration enables efficient thermal transport
792 + from surface Au layer to the capping agents. Therefore, in our
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
800 + when All-Atom models were used and C-H vibrations were treated
801 + classically. Compared to metal and organic liquid phase, the heat
802 + transfer efficiency between butanethiol and organic solvents is closer
803 + to that within bulk liquid phase.
804 +
805 + Furthermore, our observation validated previous
806 + results\cite{hase:2010} that the intramolecular heat transport of
807 + alkylthiols is highly effecient. As a combinational effects of these
808 + phenomena, butanethiol acts as a channel to expedite thermal transport
809 + process. The acoustic impedance mismatch between the metal and the
810 + liquid phase can be effectively reduced with the presence of suitable
811 + capping agents.
812 +
813   \begin{figure}
814   \includegraphics[width=\linewidth]{vibration}
815   \caption{Vibrational spectra obtained for gold in different
816 <  environments (upper panel) and for Au/thiol/hexane simulation in
817 <  all-atom model (lower panel).}
521 < \label{vibration}
816 >  environments.}
817 > \label{specAu}
818   \end{figure}
523 – % 600dpi, letter size. too large?
819  
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
843 + effective unphysical thermal flux transferred between the metal and
844 + the liquid phase. With the flux applied, we were able to measure the
845 + corresponding thermal gradient and to obtain interfacial thermal
846 + conductivities. Under steady states, single trajectory simulation
847 + would be enough for accurate measurement. This would be advantageous
848 + compared to transient state simulations, which need multiple
849 + trajectories to produce reliable average results.
850 +
851 + Our simulations have seen significant conductance enhancement with the
852 + presence of capping agent, compared to the bare gold / liquid
853 + interfaces. The acoustic impedance mismatch between the metal and the
854 + liquid phase is effectively eliminated by proper capping
855 + agent. Furthermore, the coverage precentage of the capping agent plays
856 + an important role in the interfacial thermal transport
857 + process. Moderately lower coverages allow higher contact between
858 + capping agent and solvent, and thus could further enhance the heat
859 + transfer process.
860 +
861 + Our measurement results, particularly of the UA models, agree with
862 + available experimental data. This indicates that our force field
863 + parameters have a nice description of the interactions between the
864 + particles at the interfaces. AA models tend to overestimate the
865 + interfacial thermal conductance in that the classically treated C-H
866 + vibration would be overly sampled. Compared to the AA models, the UA
867 + models have higher computational efficiency with satisfactory
868 + accuracy, and thus are preferable in interfacial thermal transport
869 + modelings. Of the two definitions for $G$, the discrete form
870 + (Eq. \ref{discreteG}) was easier to use and gives out relatively
871 + consistent results, while the derivative form (Eq. \ref{derivativeG})
872 + is not as versatile. Although $G^\prime$ gives out comparable results
873 + and follows similar trend with $G$ when measuring close to fully
874 + covered or bare surfaces, the spatial resolution of $T$ profile is
875 + limited for accurate computation of derivatives data.
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{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