ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/interfacial/interfacial.tex

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

Diff Legend

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