ViewVC Help

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

root/group/interfacial/interfacial.tex

Comparing interfacial/interfacial.tex (file contents):
Revision 3746 by skuang, Fri Jul 22 03:35:19 2011 UTC vs.
Revision 3761 by gezelter, Fri Sep 23 17:31:50 2011 UTC

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

Diff Legend

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