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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines