ViewVC Help

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

root/group/interfacial/interfacial.tex

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

Diff Legend

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