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 3734 by skuang, Mon Jul 11 18:19:26 2011 UTC vs.
Revision 3749 by skuang, Mon Jul 25 19:11:33 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 74 | Line 74 | experimentally and computationally, due to its importa
74  
75   \section{Introduction}
76   Interfacial thermal conductance is extensively studied both
77 < experimentally and computationally, due to its importance in nanoscale
78 < science and technology. Reliability of nanoscale devices depends on
79 < their thermal transport properties. Unlike bulk homogeneous materials,
80 < nanoscale materials features significant presence of interfaces, and
81 < these interfaces could dominate the heat transfer behavior of these
77 > experimentally and computationally\cite{cahill:793}, due to its
78 > importance in nanoscale science and technology. Reliability of
79 > nanoscale devices depends on their thermal transport
80 > properties. Unlike bulk homogeneous materials, nanoscale materials
81 > features significant presence of interfaces, and these interfaces
82 > could dominate the heat transfer behavior of these
83   materials. Furthermore, these materials are generally heterogeneous,
84 < which challenges traditional research methods for homogeneous systems.
84 > which challenges traditional research methods for homogeneous
85 > systems.
86  
87   Heat conductance of molecular and nano-scale interfaces will be
88   affected by the chemical details of the surface. Experimentally,
# Line 97 | Line 99 | Theoretical and computational studies were also engage
99   that specific ligands (capping agents) could completely eliminate this
100   barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}.
101  
102 < Theoretical and computational studies were also engaged in the
103 < interfacial thermal transport research in order to gain an
104 < understanding of this phenomena at the molecular level. Hase and
105 < coworkers employed Non-Equilibrium Molecular Dynamics (NEMD)
106 < simulations to study thermal transport from hot Au(111) substrate to a
107 < self-assembled monolayer of alkylthiolate with relatively long chain
108 < (8-20 carbon atoms)[CITE TWO PAPERS]. However, emsemble average measurements for heat
109 < conductance of interfaces between the capping monolayer on Au and a
110 < solvent phase has yet to be studied. The relatively low thermal flux
111 < through interfaces is difficult to measure with Equilibrium MD or
112 < forward NEMD simulation methods. Therefore, the Reverse NEMD (RNEMD)
113 < methods would have the advantage of having this difficult to measure
114 < flux known when studying the thermal transport
115 < across interfaces, given that the simulation
102 > Theoretical and computational models have also been used to study the
103 > interfacial thermal transport in order to gain an understanding of
104 > this phenomena at the molecular level. Recently, Hase and coworkers
105 > employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
106 > study thermal transport from hot Au(111) substrate to a self-assembled
107 > monolayer of alkylthiol with relatively long chain (8-20 carbon
108 > atoms)\cite{hase:2010,hase:2011}. However, ensemble averaged
109 > measurements for heat conductance of interfaces between the capping
110 > monolayer on Au and a solvent phase has yet to be studied.
111 > The comparatively low thermal flux through interfaces is
112 > difficult to measure with Equilibrium MD or forward NEMD simulation
113 > methods. Therefore, the Reverse NEMD (RNEMD) methods would have the
114 > advantage of having this difficult to measure flux known when studying
115 > the thermal transport across interfaces, given that the simulation
116   methods being able to effectively apply an unphysical flux in
117   non-homogeneous systems.
118  
# Line 119 | Line 121 | the two slabs. Furthermore, it allows more effective t
121   retains the desirable features of RNEMD (conservation of linear
122   momentum and total energy, compatibility with periodic boundary
123   conditions) while establishing true thermal distributions in each of
124 < the two slabs. Furthermore, it allows more effective thermal exchange
125 < between particles of different identities, and thus enables extensive
126 < study of interfacial conductance under steady states.
124 > the two slabs. Furthermore, it allows effective thermal exchange
125 > between particles of different identities, and thus makes the study of
126 > interfacial conductance much simpler.
127  
128 < Our work presented here investigated the Au(111) surface with various
129 < coverage of butanethiol, a capping agent with shorter carbon chain,
130 < solvated with organic solvents of different molecular shapes. And
131 < different models were used for both the capping agent and the solvent
132 < force field parameters. With the NIVS algorithm applied, the thermal
133 < transport through these interfacial systems was studied and the
134 < underlying mechanism for this phenomena was investigated.
128 > The work presented here deals with the Au(111) surface covered to
129 > varying degrees by butanethiol, a capping agent with short carbon
130 > chain, and solvated with organic solvents of different molecular
131 > properties. Different models were used for both the capping agent and
132 > the solvent force field parameters. Using the NIVS algorithm, the
133 > thermal transport across these interfaces was studied and the
134 > underlying mechanism for the phenomena was investigated.
135  
136 < [WHY STUDY AU-THIOL SURFACE; MAY CITE SHAOYI JIANG]
136 > [MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]
137  
138   \section{Methodology}
139 < \subsection{Algorithm}
140 < [BACKGROUND FOR MD METHODS]
141 < There have been many algorithms for computing thermal conductivity
142 < using molecular dynamics simulations. However, interfacial conductance
143 < is at least an order of magnitude smaller. This would make the
144 < calculation even more difficult for those slowly-converging
145 < equilibrium methods. Imposed-flux non-equilibrium
139 > \subsection{Imposd-Flux Methods in MD Simulations}
140 > Steady state MD simulations has the advantage that not many
141 > trajectories are needed to study the relationship between thermal flux
142 > and thermal gradients. For systems including low conductance
143 > interfaces one must have a method capable of generating or measuring
144 > relatively small fluxes, compared to those required for bulk
145 > conductivity. This requirement makes the calculation even more
146 > difficult for those slowly-converging equilibrium
147 > methods\cite{Viscardy:2007lq}. Forward methods may impose gradient,
148 > but in interfacial conditions it is not clear what behavior to impose
149 > at the interfacial boundaries. Imposed-flux reverse non-equilibrium
150   methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
151 < the response of temperature or momentum gradients are easier to
152 < measure than the flux, if unknown, and thus, is a preferable way to
153 < the forward NEMD methods. Although the momentum swapping approach for
154 < flux-imposing can be used for exchanging energy between particles of
155 < different identity, the kinetic energy transfer efficiency is affected
156 < by the mass difference between the particles, which limits its
151 < application on heterogeneous interfacial systems.
151 > the thermal response becomes easier to measure than the flux. Although
152 > M\"{u}ller-Plathe's original momentum swapping approach can be used
153 > for exchanging energy between particles of different identity, the
154 > kinetic energy transfer efficiency is affected by the mass difference
155 > between the particles, which limits its application on heterogeneous
156 > interfacial systems.
157  
158 < The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach in
159 < non-equilibrium MD simulations is able to impose relatively large
160 < kinetic energy flux without obvious perturbation to the velocity
161 < distribution of the simulated systems. Furthermore, this approach has
158 > The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach to
159 > non-equilibrium MD simulations is able to impose a wide range of
160 > kinetic energy fluxes without obvious perturbation to the velocity
161 > distributions of the simulated systems. Furthermore, this approach has
162   the advantage in heterogeneous interfaces in that kinetic energy flux
163   can be applied between regions of particles of arbitary identity, and
164 < the flux quantity is not restricted by particle mass difference.
164 > the flux will not be restricted by difference in particle mass.
165  
166   The NIVS algorithm scales the velocity vectors in two separate regions
167   of a simulation system with respective diagonal scaling matricies. To
168   determine these scaling factors in the matricies, a set of equations
169   including linear momentum conservation and kinetic energy conservation
170 < constraints and target momentum/energy flux satisfaction is
171 < solved. With the scaling operation applied to the system in a set
172 < frequency, corresponding momentum/temperature gradients can be built,
173 < which can be used for computing transportation properties and other
174 < applications related to momentum/temperature gradients. The NIVS
170 < algorithm conserves momenta and energy and does not depend on an
171 < external thermostat.
170 > constraints and target energy flux satisfaction is solved. With the
171 > scaling operation applied to the system in a set frequency, bulk
172 > temperature gradients can be easily established, and these can be used
173 > for computing thermal conductivities. The NIVS algorithm conserves
174 > momenta and energy and does not depend on an external thermostat.
175  
176   \subsection{Defining Interfacial Thermal Conductivity $G$}
177 < For interfaces with a relatively low interfacial conductance, the bulk
178 < regions on either side of an interface rapidly come to a state in
179 < which the two phases have relatively homogeneous (but distinct)
180 < temperatures. The interfacial thermal conductivity $G$ can therefore
181 < be approximated as:
177 > Given a system with thermal gradients and the corresponding thermal
178 > flux, for interfaces with a relatively low interfacial conductance,
179 > the bulk regions on either side of an interface rapidly come to a
180 > state in which the two phases have relatively homogeneous (but
181 > distinct) temperatures. The interfacial thermal conductivity $G$ can
182 > therefore be approximated as:
183   \begin{equation}
184   G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
185      \langle T_\mathrm{cold}\rangle \right)}
# Line 186 | Line 190 | When the interfacial conductance is {\it not} small, t
190    T_\mathrm{cold}\rangle}$ are the average observed temperature of the
191   two separated phases.
192  
193 < When the interfacial conductance is {\it not} small, two ways can be
194 < used to define $G$.
193 > When the interfacial conductance is {\it not} small, there are two
194 > ways to define $G$.
195  
196 < One way is to assume the temperature is discretely different on two
197 < sides of the interface, $G$ can be calculated with the thermal flux
198 < applied $J$ and the maximum temperature difference measured along the
199 < thermal gradient max($\Delta T$), which occurs at the interface, as:
196 > One way is to assume the temperature is discrete on the two sides of
197 > the interface. $G$ can be calculated using the applied thermal flux
198 > $J$ and the maximum temperature difference measured along the thermal
199 > gradient max($\Delta T$), which occurs at the Gibbs deviding surface
200 > (Figure \ref{demoPic}):
201   \begin{equation}
202   G=\frac{J}{\Delta T}
203   \label{discreteG}
204   \end{equation}
205  
206 + \begin{figure}
207 + \includegraphics[width=\linewidth]{method}
208 + \caption{Interfacial conductance can be calculated by applying an
209 +  (unphysical) kinetic energy flux between two slabs, one located
210 +  within the metal and another on the edge of the periodic box.  The
211 +  system responds by forming a thermal response or a gradient.  In
212 +  bulk liquids, this gradient typically has a single slope, but in
213 +  interfacial systems, there are distinct thermal conductivity
214 +  domains.  The interfacial conductance, $G$ is found by measuring the
215 +  temperature gap at the Gibbs dividing surface, or by using second
216 +  derivatives of the thermal profile.}
217 + \label{demoPic}
218 + \end{figure}
219 +
220   The other approach is to assume a continuous temperature profile along
221   the thermal gradient axis (e.g. $z$) and define $G$ at the point where
222   the magnitude of thermal conductivity $\lambda$ change reach its
# Line 213 | Line 232 | difference method and thus calculate $G^\prime$.
232  
233   With the temperature profile obtained from simulations, one is able to
234   approximate the first and second derivatives of $T$ with finite
235 < difference method and thus calculate $G^\prime$.
235 > difference methods and thus calculate $G^\prime$.
236  
237 < In what follows, both definitions are used for calculation and comparison.
237 > In what follows, both definitions have been used for calculation and
238 > are compared in the results.
239  
240 < [IMPOSE G DEFINITION INTO OUR SYSTEMS]
241 < To facilitate the use of the above definitions in calculating $G$ and
242 < $G^\prime$, we have a metal slab with its (111) surfaces perpendicular
243 < to the $z$-axis of our simulation cells. With or withour capping
244 < agents on the surfaces, the metal slab is solvated with organic
225 < solvents, as illustrated in Figure \ref{demoPic}.
240 > To compare the above definitions ($G$ and $G^\prime$), we have modeled
241 > a metal slab with its (111) surfaces perpendicular to the $z$-axis of
242 > our simulation cells. Both with and without capping agents on the
243 > surfaces, the metal slab is solvated with simple organic solvents, as
244 > illustrated in Figure \ref{gradT}.
245  
246 < \begin{figure}
247 < \includegraphics[width=\linewidth]{demoPic}
248 < \caption{A sample showing how a metal slab has its (111) surface
249 <  covered by capping agent molecules and solvated by hexane.}
250 < \label{demoPic}
251 < \end{figure}
246 > With the simulation cell described above, we are able to equilibrate
247 > the system and impose an unphysical thermal flux between the liquid
248 > and the metal phase using the NIVS algorithm. By periodically applying
249 > the unphysical flux, we are able to obtain a temperature profile and
250 > its spatial derivatives. These quantities enable the evaluation of the
251 > interfacial thermal conductance of a surface. Figure \ref{gradT} is an
252 > example of how an applied thermal flux can be used to obtain the 1st
253 > and 2nd derivatives of the temperature profile.
254  
234 With a simulation cell setup following the above manner, one is able
235 to equilibrate the system and impose an unphysical thermal flux
236 between the liquid and the metal phase with the NIVS algorithm. Under
237 a stablized thermal gradient induced by periodically applying the
238 unphysical flux, one is able to obtain a temperature profile and the
239 physical thermal flux corresponding to it, which equals to the
240 unphysical flux applied by NIVS. These data enables the evaluation of
241 the interfacial thermal conductance of a surface. Figure \ref{gradT}
242 is an example how those stablized thermal gradient can be used to
243 obtain the 1st and 2nd derivatives of the temperature profile.
244
255   \begin{figure}
256   \includegraphics[width=\linewidth]{gradT}
257 < \caption{The 1st and 2nd derivatives of temperature profile can be
258 <  obtained with finite difference approximation.}
257 > \caption{A sample of Au-butanethiol/hexane interfacial system and the
258 >  temperature profile after a kinetic energy flux is imposed to
259 >  it. The 1st and 2nd derivatives of the temperature profile can be
260 >  obtained with finite difference approximation (lower panel).}
261   \label{gradT}
262   \end{figure}
263  
252 [MAY INCLUDE POWER SPECTRUM PROTOCOL]
253
264   \section{Computational Details}
265   \subsection{Simulation Protocol}
266 < In our simulations, Au is used to construct a metal slab with bare
267 < (111) surface perpendicular to the $z$-axis. Different slab thickness
268 < (layer numbers of Au) are simulated. This metal slab is first
269 < equilibrated under normal pressure (1 atm) and a desired
270 < temperature. After equilibration, butanethiol is used as the capping
271 < agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
272 < atoms in the butanethiol molecules would occupy the three-fold sites
273 < of the surfaces, and the maximal butanethiol capacity on Au surface is
274 < $1/3$ of the total number of surface Au atoms[CITATION]. A series of
275 < different coverage surfaces is investigated in order to study the
276 < relation between coverage and conductance.
277 <
278 < [COVERAGE DISCRIPTION] However, since the interactions between surface
279 < Au and butanethiol is non-bonded, the capping agent molecules are
280 < allowed to migrate to an empty neighbor three-fold site during a
281 < simulation. Therefore, the initial configuration would not severely
272 < affect the sampling of a variety of configurations of the same
273 < coverage, and the final conductance measurement would be an average
274 < effect of these configurations explored in the simulations. [MAY NEED FIGURES]
266 > The NIVS algorithm has been implemented in our MD simulation code,
267 > OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
268 > simulations. Different metal slab thickness (layer numbers of Au) was
269 > simulated. Metal slabs were first equilibrated under atmospheric
270 > pressure (1 atm) and a desired temperature (e.g. 200K). After
271 > equilibration, butanethiol capping agents were placed at three-fold
272 > hollow sites on the Au(111) surfaces. These sites could be either a
273 > {\it fcc} or {\it hcp} site. However, Hase {\it et al.} found that
274 > they are equivalent in a heat transfer process\cite{hase:2010}, so
275 > they are not distinguished in our study. The maximum butanethiol
276 > capacity on Au surface is $1/3$ of the total number of surface Au
277 > atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
278 > structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
279 > series of different coverages was derived by evenly eliminating
280 > butanethiols on the surfaces, and was investigated in order to study
281 > the relation between coverage and interfacial conductance.
282  
283 < After the modified Au-butanethiol surface systems are equilibrated
284 < under canonical ensemble, Packmol\cite{packmol} is used to pack
285 < organic solvent molecules in the previously vacuum part of the
286 < simulation cells, which guarantees that short range repulsive
287 < interactions do not disrupt the simulations. Two solvents are
288 < investigated, one which has little vibrational overlap with the
289 < alkanethiol and plane-like shape (toluene), and one which has similar
283 < vibrational frequencies and chain-like shape ({\it n}-hexane). [MAY
284 < EXPLAIN WHY WE CHOOSE THEM]
283 > The capping agent molecules were allowed to migrate during the
284 > simulations. They distributed themselves uniformly and sampled a
285 > number of three-fold sites throughout out study. Therefore, the
286 > initial configuration would not noticeably affect the sampling of a
287 > variety of configurations of the same coverage, and the final
288 > conductance measurement would be an average effect of these
289 > configurations explored in the simulations. [MAY NEED SNAPSHOTS]
290  
291 < The spacing filled by solvent molecules, i.e. the gap between
291 > After the modified Au-butanethiol surface systems were equilibrated
292 > under canonical ensemble, organic solvent molecules were packed in the
293 > previously empty part of the simulation cells\cite{packmol}. Two
294 > solvents were investigated, one which has little vibrational overlap
295 > with the alkanethiol and a planar shape (toluene), and one which has
296 > similar vibrational frequencies and chain-like shape ({\it n}-hexane).
297 >
298 > The space filled by solvent molecules, i.e. the gap between
299   periodically repeated Au-butanethiol surfaces should be carefully
300   chosen. A very long length scale for the thermal gradient axis ($z$)
301   may cause excessively hot or cold temperatures in the middle of the
# Line 292 | Line 304 | corresponding spacing is usually $35 \sim 60$\AA.
304   solvent molecules would change the normal behavior of the liquid
305   phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
306   these extreme cases did not happen to our simulations. And the
307 < corresponding spacing is usually $35 \sim 60$\AA.
307 > corresponding spacing is usually $35 \sim 75$\AA.
308  
309 < The initial configurations generated by Packmol are further
310 < equilibrated with the $x$ and $y$ dimensions fixed, only allowing
311 < length scale change in $z$ dimension. This is to ensure that the
312 < equilibration of liquid phase does not affect the metal crystal
313 < structure in $x$ and $y$ dimensions. Further equilibration are run
314 < under NVT and then NVE ensembles.
309 > The initial configurations generated are further equilibrated with the
310 > $x$ and $y$ dimensions fixed, only allowing length scale change in $z$
311 > dimension. This is to ensure that the equilibration of liquid phase
312 > does not affect the metal crystal structure in $x$ and $y$ dimensions.
313 > To investigate this effect, comparisons were made with simulations
314 > that allow changes of $L_x$ and $L_y$ during NPT equilibration, and
315 > the results are shown in later sections. After ensuring the liquid
316 > phase reaches equilibrium at atmospheric pressure (1 atm), further
317 > equilibration are followed under NVT and then NVE ensembles.
318  
319   After the systems reach equilibrium, NIVS is implemented to impose a
320   periodic unphysical thermal flux between the metal and the liquid
321   phase. Most of our simulations are under an average temperature of
322   $\sim$200K. Therefore, this flux usually comes from the metal to the
323   liquid so that the liquid has a higher temperature and would not
324 < freeze due to excessively low temperature. This induced temperature
325 < gradient is stablized and the simulation cell is devided evenly into
326 < N slabs along the $z$-axis and the temperatures of each slab are
327 < recorded. When the slab width $d$ of each slab is the same, the
328 < derivatives of $T$ with respect to slab number $n$ can be directly
329 < used for $G^\prime$ calculations:
324 > freeze due to excessively low temperature. After this induced
325 > temperature gradient is stablized, the temperature profile of the
326 > simulation cell is recorded. To do this, the simulation cell is
327 > devided evenly into $N$ slabs along the $z$-axis and $N$ is maximized
328 > for highest possible spatial resolution but not too many to have some
329 > slabs empty most of the time. The average temperatures of each slab
330 > are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
331 > the same, the derivatives of $T$ with respect to slab number $n$ can
332 > be directly used for $G^\prime$ calculations:
333   \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
# Line 322 | 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. And
344 + each equilibration / stabilization step usually takes 100 ps, or
345 + longer, if necessary.
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 various components. Figure \ref{demoMol}
349 > demonstrates the sites defined for both United-Atom and All-Atom
350 > models of the organic solvent and capping agent molecules in our
351 > simulations. Force field parameter descriptions are needed for
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} (UA) and
362 +  \protect\cite{OPLSAA} (AA).  Cross-interactions with the Au atoms are given
363 +  in Table \ref{MnM}.}
364 + \label{demoMol}
365 + \end{figure}
366 +
367   The Au-Au interactions in metal lattice slab is described by the
368 < quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
368 > quantum Sutton-Chen (QSC) formulation\cite{PhysRevB.59.3527}. The QSC
369   potentials include zero-point quantum corrections and are
370   reparametrized for accurate surface energies compared to the
371   Sutton-Chen potentials\cite{Chen90}.
372  
336 Figure [REF] demonstrates how we name our pseudo-atoms of the
337 molecules in our simulations.
338 [FIGURE FOR MOLECULE NOMENCLATURE]
339
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
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
349 < values. Considering that after an unphysical thermal flux is applied
350 < to a system, the temperature of ``hot'' area in the liquid phase would be
351 < significantly higher than the average, to prevent over heating and
352 < boiling of the liquid phase, the average temperature in our
353 < simulations should be much lower than the liquid boiling point. [MORE DISCUSSION]
354 < For UA-toluene model, rigid body constraints are applied, so that the
355 < benzene ring and the methyl-CRar bond are kept rigid. This would save
356 < computational time.[MORE DETAILS]
377 > for our UA solvent molecules. In these models, sites are located at
378 > the carbon centers for alkyl groups. Bonding interactions, including
379 > bond stretches and bends and torsions, were used for intra-molecular
380 > sites not separated by more than 3 bonds. Otherwise, for non-bonded
381 > interactions, Lennard-Jones potentials are used. [CHECK CITATION]
382  
383 + By eliminating explicit hydrogen atoms, these models are simple and
384 + computationally efficient, while maintains good accuracy. However, the
385 + TraPPE-UA for alkanes is known to predict a lower boiling point than
386 + experimental values. Considering that after an unphysical thermal flux
387 + is applied to a system, the temperature of ``hot'' area in the liquid
388 + phase would be significantly higher than the average of the system, to
389 + prevent over heating and boiling of the liquid phase, the average
390 + temperature in our simulations should be much lower than the liquid
391 + boiling point.
392 +
393 + For UA-toluene model, the non-bonded potentials between
394 + inter-molecular sites have a similar Lennard-Jones formulation. For
395 + intra-molecular interactions, considering the stiffness of the benzene
396 + ring, rigid body constraints are applied for further computational
397 + efficiency. All bonds in the benzene ring and between the ring and the
398 + methyl group remain rigid during the progress of simulations.
399 +
400   Besides the TraPPE-UA models, AA models for both organic solvents are
401   included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
402 < force field is used. [MORE DETAILS]
403 < For toluene, the United Force Field developed by Rapp\'{e} {\it et
404 <  al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]
402 > force field is used. Additional explicit hydrogen sites were
403 > included. Besides bonding and non-bonded site-site interactions,
404 > partial charges and the electrostatic interactions were added to each
405 > CT and HC site. For toluene, the United Force Field developed by
406 > Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} is
407 > adopted. Without the rigid body constraints, bonding interactions were
408 > included. For the aromatic ring, improper torsions (inversions) were
409 > added as an extra potential for maintaining the planar shape.
410 > [CHECK CITATION]
411  
412   The capping agent in our simulations, the butanethiol molecules can
413   either use UA or AA model. The TraPPE-UA force fields includes
# Line 368 | Line 416 | Au(111) surfaces, we adopt the S parameters from [CITA
416   parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
417   surfaces do not have the hydrogen atom bonded to sulfur. To adapt this
418   change and derive suitable parameters for butanethiol adsorbed on
419 < Au(111) surfaces, we adopt the S parameters from [CITATION CF VLUGT]
420 < and modify parameters for its neighbor C atom for charge balance in
421 < the molecule. Note that the model choice (UA or AA) of capping agent
422 < can be different from the solvent. Regardless of model choice, the
423 < force field parameters for interactions between capping agent and
424 < solvent can be derived using Lorentz-Berthelot Mixing Rule:
419 > Au(111) surfaces, we adopt the S parameters from Luedtke and
420 > Landman\cite{landman:1998}[CHECK CITATION]
421 > and modify parameters for its neighbor C
422 > atom for charge balance in the molecule. Note that the model choice
423 > (UA or AA) of capping agent can be different from the
424 > solvent. Regardless of model choice, the force field parameters for
425 > interactions between capping agent and solvent can be derived using
426 > Lorentz-Berthelot Mixing Rule:
427 > \begin{eqnarray}
428 > \sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
429 > \epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
430 > \end{eqnarray}
431  
378
432   To describe the interactions between metal Au and non-metal capping
433   agent and solvent particles, we refer to an adsorption study of alkyl
434   thiols on gold surfaces by Vlugt {\it et
435    al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones
436   form of potential parameters for the interaction between Au and
437   pseudo-atoms CH$_x$ and S based on a well-established and widely-used
438 < effective potential of Hautman and Klein[CITATION] for the Au(111)
439 < surface. As our simulations require the gold lattice slab to be
440 < non-rigid so that it could accommodate kinetic energy for thermal
438 > effective potential of Hautman and Klein\cite{hautman:4994} for the
439 > Au(111) surface. As our simulations require the gold lattice slab to
440 > be non-rigid so that it could accommodate kinetic energy for thermal
441   transport study purpose, the pair-wise form of potentials is
442   preferred.
443  
444   Besides, the potentials developed from {\it ab initio} calculations by
445   Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
446 < interactions between Au and aromatic C/H atoms in toluene.[MORE DETAILS]
446 > interactions between Au and aromatic C/H atoms in toluene. A set of
447 > pseudo Lennard-Jones parameters were provided for Au in their force
448 > fields. By using the Mixing Rule, this can be used to derive pair-wise
449 > potentials for non-bonded interactions between Au and non-metal sites.
450  
451   However, the Lennard-Jones parameters between Au and other types of
452 < particles in our simulations are not yet well-established. For these
453 < interactions, we attempt to derive their parameters using the Mixing
454 < Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters
455 < for Au is first extracted from the Au-CH$_x$ parameters by applying
456 < the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
452 > particles, such as All-Atom normal alkanes in our simulations are not
453 > yet well-established. For these interactions, we attempt to derive
454 > their parameters using the Mixing Rule. To do this, Au pseudo
455 > Lennard-Jones parameters for ``Metal-non-Metal'' (MnM) interactions
456 > were first extracted from the Au-CH$_x$ parameters by applying the
457 > Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
458   parameters in our simulations.
459  
460   \begin{table*}
461    \begin{minipage}{\linewidth}
462      \begin{center}
463 <      \caption{Lennard-Jones parameters for Au-non-Metal
464 <        interactions in our simulations.}
465 <      
466 <      \begin{tabular}{ccc}
463 >      \caption{Non-bonded interaction parameters (including cross
464 >        interactions with Au atoms) for both force fields used in this
465 >        work.}      
466 >      \begin{tabular}{lllllll}
467          \hline\hline
468 <        Non-metal atom   & $\sigma$ & $\epsilon$ \\
469 <        (or pseudo-atom) & \AA      & kcal/mol  \\
468 >        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
469 >        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
470 >        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
471          \hline
472 <        S    & 2.40   & 8.465   \\
473 <        CH3  & 3.54   & 0.2146  \\
474 <        CH2  & 3.54   & 0.1749  \\
475 <        CT3  & 3.365  & 0.1373  \\
476 <        CT2  & 3.365  & 0.1373  \\
477 <        CTT  & 3.365  & 0.1373  \\
478 <        HC   & 2.865  & 0.09256 \\
479 <        CHar & 3.4625 & 0.1680  \\
480 <        CRar & 3.555  & 0.1604  \\
481 <        CA   & 3.173  & 0.0640  \\
482 <        HA   & 2.746  & 0.0414  \\
472 >        United Atom (UA)
473 >        &CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
474 >        &CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
475 >        &CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
476 >        &CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
477 >        \hline
478 >        All Atom (AA)
479 >        &CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
480 >        &CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
481 >        &CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
482 >        &HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
483 >        &CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
484 >        &HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
485 >        \hline
486 >        Both UA and AA
487 >        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
488          \hline\hline
489        \end{tabular}
490        \label{MnM}
# Line 429 | Line 492 | parameters in our simulations.
492    \end{minipage}
493   \end{table*}
494  
495 + \subsection{Vibrational Spectrum}
496 + To investigate the mechanism of interfacial thermal conductance, the
497 + vibrational spectrum is utilized as a complementary tool. Vibrational
498 + spectra were taken for individual components in different
499 + simulations. To obtain these spectra, simulations were run after
500 + equilibration, in the NVE ensemble. Snapshots of configurations were
501 + collected at a frequency that is higher than that of the fastest
502 + vibrations occuring in the simulations. With these configurations, the
503 + velocity auto-correlation functions can be computed:
504 + \begin{equation}
505 + C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
506 + \label{vCorr}
507 + \end{equation}
508 + Followed by Fourier transforms, the power spectrum can be constructed:
509 + \begin{equation}
510 + \hat{f}(\omega) = \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
511 + \label{fourier}
512 + \end{equation}
513  
514   \section{Results and Discussions}
515 < [MAY HAVE A BRIEF SUMMARY]
515 > In what follows, how the parameters and protocol of simulations would
516 > affect the measurement of $G$'s is first discussed. With a reliable
517 > protocol and set of parameters, the influence of capping agent
518 > coverage on thermal conductance is investigated. Besides, different
519 > force field models for both solvents and selected deuterated models
520 > were tested and compared. Finally, a summary of the role of capping
521 > agent in the interfacial thermal transport process is given.
522 >
523   \subsection{How Simulation Parameters Affects $G$}
436 [MAY NOT PUT AT FIRST]
524   We have varied our protocol or other parameters of the simulations in
525   order to investigate how these factors would affect the measurement of
526   $G$'s. It turned out that while some of these parameters would not
# Line 442 | Line 529 | during equilibrating the liquid phase. Due to the stif
529   results.
530  
531   In some of our simulations, we allowed $L_x$ and $L_y$ to change
532 < during equilibrating the liquid phase. Due to the stiffness of the Au
533 < slab, $L_x$ and $L_y$ would not change noticeably after
534 < equilibration. Although $L_z$ could fluctuate $\sim$1\% after a system
535 < is fully equilibrated in the NPT ensemble, this fluctuation, as well
536 < as those comparably smaller to $L_x$ and $L_y$, would not be magnified
537 < on the calculated $G$'s, as shown in Table \ref{AuThiolHexaneUA}. This
538 < insensivity to $L_i$ fluctuations allows reliable measurement of $G$'s
539 < without the necessity of extremely cautious equilibration process.
532 > during equilibrating the liquid phase. Due to the stiffness of the
533 > crystalline Au structure, $L_x$ and $L_y$ would not change noticeably
534 > after equilibration. Although $L_z$ could fluctuate $\sim$1\% after a
535 > system is fully equilibrated in the NPT ensemble, this fluctuation, as
536 > well as those of $L_x$ and $L_y$ (which is significantly smaller),
537 > would not be magnified on the calculated $G$'s, as shown in Table
538 > \ref{AuThiolHexaneUA}. This insensivity to $L_i$ fluctuations allows
539 > reliable measurement of $G$'s without the necessity of extremely
540 > cautious equilibration process.
541  
542   As stated in our computational details, the spacing filled with
543   solvent molecules can be chosen within a range. This allows some
# Line 476 | Line 564 | $J_z$ is listed in Table \ref{AuThiolHexaneUA} and [RE
564   the thermal flux across the interface. For our simulations, we denote
565   $J_z$ to be positive when the physical thermal flux is from the liquid
566   to metal, and negative vice versa. The $G$'s measured under different
567 < $J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These
568 < results do not suggest that $G$ is dependent on $J_z$ within this flux
569 < range. The linear response of flux to thermal gradient simplifies our
570 < investigations in that we can rely on $G$ measurement with only a
571 < couple $J_z$'s and do not need to test a large series of fluxes.
567 > $J_z$ is listed in Table \ref{AuThiolHexaneUA} and
568 > \ref{AuThiolToluene}. These results do not suggest that $G$ is
569 > dependent on $J_z$ within this flux range. The linear response of flux
570 > to thermal gradient simplifies our investigations in that we can rely
571 > on $G$ measurement with only a couple $J_z$'s and do not need to test
572 > a large series of fluxes.
573  
485 %ADD MORE TO TABLE
574   \begin{table*}
575    \begin{minipage}{\linewidth}
576      \begin{center}
577        \caption{Computed interfacial thermal conductivity ($G$ and
578          $G^\prime$) values for the 100\% covered Au-butanethiol/hexane
579          interfaces with UA model and different hexane molecule numbers
580 <        at different temperatures using a range of energy fluxes.}
580 >        at different temperatures using a range of energy
581 >        fluxes. Error estimates indicated in parenthesis.}
582        
583 <      \begin{tabular}{cccccccc}
583 >      \begin{tabular}{ccccccc}
584          \hline\hline
585 <        $\langle T\rangle$ & & $L_x$ & $L_y$ & $L_z$ & $J_z$ &
586 <        $G$ & $G^\prime$ \\
587 <        (K) & $N_{hexane}$ & \multicolumn{3}{c}{(\AA)} & (GW/m$^2$) &
585 >        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
586 >        $J_z$ & $G$ & $G^\prime$ \\
587 >        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
588          \multicolumn{2}{c}{(MW/m$^2$/K)} \\
589          \hline
590 <        200 & 266 & 29.86 & 25.80 & 113.1 & -0.96 &
591 <        102()  & 80.0() \\
592 <            & 200 & 29.84 & 25.81 &  93.9 &  1.92 &
593 <        129()  & 87.3() \\
594 <            &     & 29.84 & 25.81 &  95.3 &  1.93 &
595 <        131()  & 77.5() \\
596 <            & 166 & 29.84 & 25.81 &  85.7 &  0.97 &
597 <        115()  & 69.3() \\
598 <            &     &       &       &       &  1.94 &
599 <        125()  & 87.1() \\
600 <        250 & 200 & 29.84 & 25.87 & 106.8 &  0.96 &
601 <        81.8() & 67.0() \\
602 <            & 166 & 29.87 & 25.84 &  94.8 &  0.98 &
603 <        79.0() & 62.9() \\
604 <            &     & 29.84 & 25.85 &  95.0 &  1.44 &
605 <        76.2() & 64.8() \\
590 >        200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
591 >            & 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
592 >            &     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
593 >            &     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
594 >            &     &     &       &  1.91 & 139(10)   & 101(10)   \\
595 >            &     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
596 >            & 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
597 >            &     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
598 >            &     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
599 >            &     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
600 >        \hline
601 >        250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
602 >            &     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
603 >            & 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
604 >            &     & No  & 0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
605 >            &     &     &       &  1.44 & 76.2(5.0) & 64.8(3.8) \\
606 >            &     &     &       & -0.95 & 56.4(2.5) & 54.4(1.1) \\
607 >            &     &     &       & -1.85 & 47.8(1.1) & 53.5(1.5) \\
608          \hline\hline
609        \end{tabular}
610        \label{AuThiolHexaneUA}
# Line 529 | Line 620 | butanethiol as well.[MAY NEED FIGURE] And this reduced
620   temperature is higher than 250K. Additionally, the equilibrated liquid
621   hexane density under 250K becomes lower than experimental value. This
622   expanded liquid phase leads to lower contact between hexane and
623 < butanethiol as well.[MAY NEED FIGURE] And this reduced contact would
623 > butanethiol as well.[MAY NEED SLAB DENSITY FIGURE]
624 > And this reduced contact would
625   probably be accountable for a lower interfacial thermal conductance,
626   as shown in Table \ref{AuThiolHexaneUA}.
627  
# Line 544 | Line 636 | in that higher degree of contact could yield increased
636   important role in the thermal transport process across the interface
637   in that higher degree of contact could yield increased conductance.
638  
547 [ADD Lxyz AND ERROR ESTIMATE TO TABLE]
639   \begin{table*}
640    \begin{minipage}{\linewidth}
641      \begin{center}
642        \caption{Computed interfacial thermal conductivity ($G$ and
643          $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
644          interface at different temperatures using a range of energy
645 <        fluxes.}
645 >        fluxes. Error estimates indicated in parenthesis.}
646        
647 <      \begin{tabular}{cccc}
647 >      \begin{tabular}{ccccc}
648          \hline\hline
649 <        $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
650 <        (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
649 >        $\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
650 >        (K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
651          \hline
652 <        200 & -1.86 & 180() & 135() \\
653 <            &  2.15 & 204() & 113() \\
654 <            & -3.93 & 175() & 114() \\
655 <        300 & -1.91 & 143() & 125() \\
656 <            & -4.19 & 134() & 113() \\
652 >        200 & 0.933 &  2.15 & 204(12) & 113(12) \\
653 >            &       & -1.86 & 180(3)  & 135(21) \\
654 >            &       & -3.93 & 176(5)  & 113(12) \\
655 >        \hline
656 >        300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
657 >            &       & -4.19 & 135(9)  & 113(12) \\
658          \hline\hline
659        \end{tabular}
660        \label{AuThiolToluene}
# Line 594 | Line 686 | above. Our Au-butanethiol/toluene system did not see t
686  
687   However, when the surface is not completely covered by butanethiols,
688   the simulated system is more resistent to the reconstruction
689 < above. Our Au-butanethiol/toluene system did not see this phenomena
690 < even at $\langle T\rangle\sim$300K. The Au(111) surfaces have a 90\% coverage of
691 < butanethiols and have empty three-fold sites. These empty sites could
692 < help prevent surface reconstruction in that they provide other means
693 < of capping agent relaxation. It is observed that butanethiols can
694 < migrate to their neighbor empty sites during a simulation. Therefore,
695 < we were able to obtain $G$'s for these interfaces even at a relatively
696 < high temperature without being affected by surface reconstructions.
689 > above. Our Au-butanethiol/toluene system had the Au(111) surfaces 90\%
690 > covered by butanethiols, but did not see this above phenomena even at
691 > $\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
692 > capping agents could help prevent surface reconstruction in that they
693 > provide other means of capping agent relaxation. It is observed that
694 > butanethiols can migrate to their neighbor empty sites during a
695 > simulation. Therefore, we were able to obtain $G$'s for these
696 > interfaces even at a relatively high temperature without being
697 > affected by surface reconstructions.
698  
699   \subsection{Influence of Capping Agent Coverage on $G$}
700   To investigate the influence of butanethiol coverage on interfacial
701   thermal conductance, a series of different coverage Au-butanethiol
702   surfaces is prepared and solvated with various organic
703   molecules. These systems are then equilibrated and their interfacial
704 < thermal conductivity are measured with our NIVS algorithm. Table
705 < \ref{tlnUhxnUhxnD} lists these results for direct comparison between
706 < different coverages of butanethiol. To study the isotope effect in
707 < interfacial thermal conductance, deuterated UA-hexane is included as
708 < well.
704 > thermal conductivity are measured with our NIVS algorithm. Figure
705 > \ref{coverage} demonstrates the trend of conductance change with
706 > respect to different coverages of butanethiol. To study the isotope
707 > effect in interfacial thermal conductance, deuterated UA-hexane is
708 > included as well.
709  
710 + \begin{figure}
711 + \includegraphics[width=\linewidth]{coverage}
712 + \caption{Comparison of interfacial thermal conductivity ($G$) values
713 +  for the Au-butanethiol/solvent interface with various UA models and
714 +  different capping agent coverages at $\langle T\rangle\sim$200K
715 +  using certain energy flux respectively.}
716 + \label{coverage}
717 + \end{figure}
718 +
719   It turned out that with partial covered butanethiol on the Au(111)
720 < surface, the derivative definition for $G$ (Eq. \ref{derivativeG}) has
721 < difficulty to apply, due to the difficulty in locating the maximum of
722 < change of $\lambda$. Instead, the discrete definition
723 < (Eq. \ref{discreteG}) is easier to apply, as max($\Delta T$) can still
724 < be well-defined. Therefore, $G$'s (not $G^\prime$) are used for this
725 < section.
720 > surface, the derivative definition for $G^\prime$
721 > (Eq. \ref{derivativeG}) was difficult to apply, due to the difficulty
722 > in locating the maximum of change of $\lambda$. Instead, the discrete
723 > definition (Eq. \ref{discreteG}) is easier to apply, as the Gibbs
724 > deviding surface can still be well-defined. Therefore, $G$ (not
725 > $G^\prime$) was used for this section.
726  
727 < From Table \ref{tlnUhxnUhxnD}, one can see the significance of the
727 > From Figure \ref{coverage}, one can see the significance of the
728   presence of capping agents. Even when a fraction of the Au(111)
729   surface sites are covered with butanethiols, the conductivity would
730   see an enhancement by at least a factor of 3. This indicates the
731   important role cappping agent is playing for thermal transport
732 < phenomena on metal/organic solvent surfaces.
732 > phenomena on metal / organic solvent surfaces.
733  
734   Interestingly, as one could observe from our results, the maximum
735   conductance enhancement (largest $G$) happens while the surfaces are
# Line 646 | Line 748 | case, $G$ decrease could not be offset but instead acc
748   would not offset this effect. Eventually, when butanethiol coverage
749   continues to decrease, solvent-capping agent contact actually
750   decreases with the disappearing of butanethiol molecules. In this
751 < case, $G$ decrease could not be offset but instead accelerated.
751 > case, $G$ decrease could not be offset but instead accelerated. [NEED
752 > SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
753  
754   A comparison of the results obtained from differenet organic solvents
755   can also provide useful information of the interfacial thermal
# Line 656 | Line 759 | difference for the results of $G$. [MAY NEED FIGURE]
759   studies, even though eliminating C-H vibration samplings, still have
760   C-C vibrational frequencies different from each other. However, these
761   differences in the infrared range do not seem to produce an observable
762 < difference for the results of $G$. [MAY NEED FIGURE]
762 > difference for the results of $G$ (Figure \ref{uahxnua}).
763  
764 + \begin{figure}
765 + \includegraphics[width=\linewidth]{uahxnua}
766 + \caption{Vibrational spectra obtained for normal (upper) and
767 +  deuterated (lower) hexane in Au-butanethiol/hexane
768 +  systems. Butanethiol spectra are shown as reference. Both hexane and
769 +  butanethiol were using United-Atom models.}
770 + \label{uahxnua}
771 + \end{figure}
772 +
773   Furthermore, results for rigid body toluene solvent, as well as other
774   UA-hexane solvents, are reasonable within the general experimental
775 < ranges[CITATIONS]. This suggests that explicit hydrogen might not be a
776 < required factor for modeling thermal transport phenomena of systems
777 < such as Au-thiol/organic solvent.
775 > ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
776 > suggests that explicit hydrogen might not be a required factor for
777 > modeling thermal transport phenomena of systems such as
778 > Au-thiol/organic solvent.
779  
780   However, results for Au-butanethiol/toluene do not show an identical
781 < trend with those for Au-butanethiol/hexane in that $G$'s remain at
781 > trend with those for Au-butanethiol/hexane in that $G$ remains at
782   approximately the same magnitue when butanethiol coverage differs from
783   25\% to 75\%. This might be rooted in the molecule shape difference
784 < for plane-like toluene and chain-like {\it n}-hexane. Due to this
784 > for planar toluene and chain-like {\it n}-hexane. Due to this
785   difference, toluene molecules have more difficulty in occupying
786   relatively small gaps among capping agents when their coverage is not
787   too low. Therefore, the solvent-capping agent contact may keep
# Line 676 | Line 789 | can see a plateau of $G$ vs. butanethiol coverage in o
789   level. This becomes an offset for decreasing butanethiol molecules on
790   its effect to the process of interfacial thermal transport. Thus, one
791   can see a plateau of $G$ vs. butanethiol coverage in our results.
679
680 [NEED ERROR ESTIMATE, MAY ALSO PUT J HERE]
681 \begin{table*}
682  \begin{minipage}{\linewidth}
683    \begin{center}
684      \caption{Computed interfacial thermal conductivity ($G$) values
685        for the Au-butanethiol/solvent interface with various UA
686        models and different capping agent coverages at $\langle
687        T\rangle\sim$200K using certain energy flux respectively.}
688      
689      \begin{tabular}{cccc}
690        \hline\hline
691        Thiol & \multicolumn{3}{c}{$G$(MW/m$^2$/K)} \\
692        coverage (\%) & hexane & hexane(D) & toluene \\
693        \hline
694        0.0   & 46.5() & 43.9() & 70.1() \\
695        25.0  & 151()  & 153()  & 249()  \\
696        50.0  & 172()  & 182()  & 214()  \\
697        75.0  & 242()  & 229()  & 244()  \\
698        88.9  & 178()  & -      & -      \\
699        100.0 & 137()  & 153()  & 187()  \\
700        \hline\hline
701      \end{tabular}
702      \label{tlnUhxnUhxnD}
703    \end{center}
704  \end{minipage}
705 \end{table*}
792  
793   \subsection{Influence of Chosen Molecule Model on $G$}
708 [MAY COMBINE W MECHANISM STUDY]
709
794   In addition to UA solvent/capping agent models, AA models are included
795   in our simulations as well. Besides simulations of the same (UA or AA)
796   model for solvent and capping agent, different models can be applied
# Line 715 | Line 799 | these studies.
799   the previous section. Table \ref{modelTest} summarizes the results of
800   these studies.
801  
718 [MORE DATA; ERROR ESTIMATE]
802   \begin{table*}
803    \begin{minipage}{\linewidth}
804      \begin{center}
# Line 723 | Line 806 | these studies.
806        \caption{Computed interfacial thermal conductivity ($G$ and
807          $G^\prime$) values for interfaces using various models for
808          solvent and capping agent (or without capping agent) at
809 <        $\langle T\rangle\sim$200K.}
809 >        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
810 >        or capping agent molecules; ``Avg.'' denotes results that are
811 >        averages of simulations under different $J_z$'s. Error
812 >        estimates indicated in parenthesis.)}
813        
814 <      \begin{tabular}{ccccc}
814 >      \begin{tabular}{llccc}
815          \hline\hline
816          Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
817          (or bare surface) & model & (GW/m$^2$) &
818          \multicolumn{2}{c}{(MW/m$^2$/K)} \\
819          \hline
820 <        UA    & AA hexane  & 1.94 & 135()  & 129()  \\
821 <              &            & 2.86 & 126()  & 115()  \\
822 <              & AA toluene & 1.89 & 200()  & 149()  \\
823 <        AA    & UA hexane  & 1.94 & 116()  & 129()  \\
824 <              & AA hexane  & 3.76 & 451()  & 378()  \\
825 <              &            & 4.71 & 432()  & 334()  \\
826 <              & AA toluene & 3.79 & 487()  & 290()  \\
827 <        AA(D) & UA hexane  & 1.94 & 158()  & 172()  \\
828 <        bare  & AA hexane  & 0.96 & 31.0() & 29.4() \\
820 >        UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
821 >              & UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
822 >              & AA hexane    & Avg. & 131(6)    & 122(10)   \\
823 >              & UA toluene   & 1.96 & 187(16)   & 151(11)   \\
824 >              & AA toluene   & 1.89 & 200(36)   & 149(53)   \\
825 >        \hline
826 >        AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
827 >              & AA hexane    & Avg. & 442(14)   & 356(31)   \\
828 >              & AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
829 >              & UA toluene   & 1.98 & 125(25)   & 97(60)    \\
830 >              & AA toluene   & 3.79 & 487(56)   & 290(42)   \\
831 >        \hline
832 >        AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
833 >              & AA hexane    & 1.92 & 243(29)   & 191(11)   \\
834 >              & AA toluene   & 1.93 & 364(36)   & 322(67)   \\
835 >        \hline
836 >        bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
837 >              & UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
838 >              & AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
839 >              & UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
840          \hline\hline
841        \end{tabular}
842        \label{modelTest}
# Line 762 | Line 859 | the upper bond of experimental value range. This is pr
859   interfaces, using AA model for both butanethiol and hexane yields
860   substantially higher conductivity values than using UA model for at
861   least one component of the solvent and capping agent, which exceeds
862 < the upper bond of experimental value range. This is probably due to
863 < the classically treated C-H vibrations in the AA model, which should
864 < not be appreciably populated at normal temperatures. In comparison,
865 < once either the hexanes or the butanethiols are deuterated, one can
866 < see a significantly lower $G$ and $G^\prime$. In either of these
867 < cases, the C-H(D) vibrational overlap between the solvent and the
868 < capping agent is removed. [MAY NEED FIGURE] Conclusively, the
869 < improperly treated C-H vibration in the AA model produced
870 < over-predicted results accordingly. Compared to the AA model, the UA
871 < model yields more reasonable results with higher computational
872 < efficiency.
862 > the general range of experimental measurement results. This is
863 > probably due to the classically treated C-H vibrations in the AA
864 > model, which should not be appreciably populated at normal
865 > temperatures. In comparison, once either the hexanes or the
866 > butanethiols are deuterated, one can see a significantly lower $G$ and
867 > $G^\prime$. In either of these cases, the C-H(D) vibrational overlap
868 > between the solvent and the capping agent is removed (Figure
869 > \ref{aahxntln}). Conclusively, the improperly treated C-H vibration in
870 > the AA model produced over-predicted results accordingly. Compared to
871 > the AA model, the UA model yields more reasonable results with higher
872 > computational efficiency.
873  
874 + \begin{figure}
875 + \includegraphics[width=\linewidth]{aahxntln}
876 + \caption{Spectra obtained for All-Atom model Au-butanethil/solvent
877 +  systems. When butanethiol is deuterated (lower left), its
878 +  vibrational overlap with hexane would decrease significantly,
879 +  compared with normal butanethiol (upper left). However, this
880 +  dramatic change does not apply to toluene as much (right).}
881 + \label{aahxntln}
882 + \end{figure}
883 +
884   However, for Au-butanethiol/toluene interfaces, having the AA
885   butanethiol deuterated did not yield a significant change in the
886 < measurement results.
887 < . , so extra degrees of freedom
888 < such as the C-H vibrations could enhance heat exchange between these
889 < two phases and result in a much higher conductivity.
886 > measurement results. Compared to the C-H vibrational overlap between
887 > hexane and butanethiol, both of which have alkyl chains, that overlap
888 > between toluene and butanethiol is not so significant and thus does
889 > not have as much contribution to the heat exchange
890 > process. Conversely, extra degrees of freedom such as the C-H
891 > vibrations could yield higher heat exchange rate between these two
892 > phases and result in a much higher conductivity.
893  
784
894   Although the QSC model for Au is known to predict an overly low value
895 < for bulk metal gold conductivity[CITE NIVSRNEMD], our computational
895 > for bulk metal gold conductivity\cite{kuang:164101}, our computational
896   results for $G$ and $G^\prime$ do not seem to be affected by this
897 < drawback of the model for metal. Instead, the modeling of interfacial
898 < thermal transport behavior relies mainly on an accurate description of
899 < the interactions between components occupying the interfaces.
897 > drawback of the model for metal. Instead, our results suggest that the
898 > modeling of interfacial thermal transport behavior relies mainly on
899 > the accuracy of the interaction descriptions between components
900 > occupying the interfaces.
901  
902 < \subsection{Mechanism of Interfacial Thermal Conductance Enhancement
903 <  by Capping Agent}
904 < %OR\subsection{Vibrational spectrum study on conductance mechanism}
902 > \subsection{Role of Capping Agent in Interfacial Thermal Conductance}
903 > The vibrational spectra for gold slabs in different environments are
904 > shown as in Figure \ref{specAu}. Regardless of the presence of
905 > solvent, the gold surfaces covered by butanethiol molecules, compared
906 > to bare gold surfaces, exhibit an additional peak observed at the
907 > frequency of $\sim$170cm$^{-1}$, which is attributed to the S-Au
908 > bonding vibration. This vibration enables efficient thermal transport
909 > from surface Au layer to the capping agents. Therefore, in our
910 > simulations, the Au/S interfaces do not appear major heat barriers
911 > compared to the butanethiol / solvent interfaces.
912  
913 < [MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL, EQN'S]
913 > Simultaneously, the vibrational overlap between butanethiol and
914 > organic solvents suggests higher thermal exchange efficiency between
915 > these two components. Even exessively high heat transport was observed
916 > when All-Atom models were used and C-H vibrations were treated
917 > classically. Compared to metal and organic liquid phase, the heat
918 > transfer efficiency between butanethiol and organic solvents is closer
919 > to that within bulk liquid phase.
920  
921 < To investigate the mechanism of this interfacial thermal conductance,
922 < the vibrational spectra of various gold systems were obtained and are
923 < shown as in the upper panel of Fig. \ref{vibration}. To obtain these
924 < spectra, one first runs a simulation in the NVE ensemble and collects
925 < snapshots of configurations; these configurations are used to compute
926 < the velocity auto-correlation functions, which is used to construct a
927 < power spectrum via a Fourier transform.
921 > Furthermore, our observation validated previous
922 > results\cite{hase:2010} that the intramolecular heat transport of
923 > alkylthiols is highly effecient. As a combinational effects of these
924 > phenomena, butanethiol acts as a channel to expedite thermal transport
925 > process. The acoustic impedance mismatch between the metal and the
926 > liquid phase can be effectively reduced with the presence of suitable
927 > capping agents.
928  
806 The gold surfaces covered by
807 butanethiol molecules, compared to bare gold surfaces, exhibit an
808 additional peak observed at a frequency of $\sim$170cm$^{-1}$, which
809 is attributed to the vibration of the S-Au bond. This vibration
810 enables efficient thermal transport from surface Au atoms to the
811 capping agents. Simultaneously, as shown in the lower panel of
812 Fig. \ref{vibration}, the large overlap of the vibration spectra of
813 butanethiol and hexane in the all-atom model, including the C-H
814 vibration, also suggests high thermal exchange efficiency. The
815 combination of these two effects produces the drastic interfacial
816 thermal conductance enhancement in the all-atom model.
817
818 [MAY NEED TO CONVERT TO JPEG]
929   \begin{figure}
930   \includegraphics[width=\linewidth]{vibration}
931   \caption{Vibrational spectra obtained for gold in different
932 <  environments (upper panel) and for Au/thiol/hexane simulation in
933 <  all-atom model (lower panel).}
824 < \label{vibration}
932 >  environments.}
933 > \label{specAu}
934   \end{figure}
935  
936 < [COMPARISON OF TWO G'S; AU SLAB WIDTHS; ETC]
828 < % The results show that the two definitions used for $G$ yield
829 < % comparable values, though $G^\prime$ tends to be smaller.
936 > [MAY ADD COMPARISON OF AU SLAB WIDTHS]
937  
938   \section{Conclusions}
939   The NIVS algorithm we developed has been applied to simulations of
# Line 834 | Line 941 | conductivities. Our simulations have seen significant
941   effective unphysical thermal flux transferred between the metal and
942   the liquid phase. With the flux applied, we were able to measure the
943   corresponding thermal gradient and to obtain interfacial thermal
944 < conductivities. Our simulations have seen significant conductance
945 < enhancement with the presence of capping agent, compared to the bare
946 < gold/liquid interfaces. The acoustic impedance mismatch between the
947 < metal and the liquid phase is effectively eliminated by proper capping
944 > conductivities. Under steady states, single trajectory simulation
945 > would be enough for accurate measurement. This would be advantageous
946 > compared to transient state simulations, which need multiple
947 > trajectories to produce reliable average results.
948 >
949 > Our simulations have seen significant conductance enhancement with the
950 > presence of capping agent, compared to the bare gold / liquid
951 > interfaces. The acoustic impedance mismatch between the metal and the
952 > liquid phase is effectively eliminated by proper capping
953   agent. Furthermore, the coverage precentage of the capping agent plays
954 < an important role in the interfacial thermal transport process.
954 > an important role in the interfacial thermal transport
955 > process. Moderately lower coverages allow higher contact between
956 > capping agent and solvent, and thus could further enhance the heat
957 > transfer process.
958  
959   Our measurement results, particularly of the UA models, agree with
960   available experimental data. This indicates that our force field
# Line 849 | Line 964 | modelings.
964   vibration would be overly sampled. Compared to the AA models, the UA
965   models have higher computational efficiency with satisfactory
966   accuracy, and thus are preferable in interfacial thermal transport
967 < modelings.
967 > modelings. Of the two definitions for $G$, the discrete form
968 > (Eq. \ref{discreteG}) was easier to use and gives out relatively
969 > consistent results, while the derivative form (Eq. \ref{derivativeG})
970 > is not as versatile. Although $G^\prime$ gives out comparable results
971 > and follows similar trend with $G$ when measuring close to fully
972 > covered or bare surfaces, the spatial resolution of $T$ profile is
973 > limited for accurate computation of derivatives data.
974  
975   Vlugt {\it et al.} has investigated the surface thiol structures for
976   nanocrystal gold and pointed out that they differs from those of the
# Line 859 | Line 980 | simulating structures with spherical symmetry.
980   and measure the corresponding thermal gradient is desirable for
981   simulating structures with spherical symmetry.
982  
862
983   \section{Acknowledgments}
984   Support for this project was provided by the National Science
985   Foundation under grant CHE-0848243. Computational time was provided by

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines