ViewVC Help

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

root/group/interfacial/interfacial.tex

Comparing interfacial/interfacial.tex (file contents):
Revision 3743 by skuang, Fri Jul 15 19:36:02 2011 UTC vs.
Revision 3764 by skuang, Wed Sep 28 22:26:23 2011 UTC

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

Diff Legend

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