ViewVC Help

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

root/group/interfacial/interfacial.tex

Comparing interfacial/interfacial.tex (file contents):
Revision 3741 by gezelter, Thu Jul 14 21:55:24 2011 UTC vs.
Revision 3760 by skuang, Mon Aug 1 17:57:34 2011 UTC

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

Diff Legend

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