ViewVC Help

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

root/group/interfacial/interfacial.tex

Comparing interfacial/interfacial.tex (file contents):
Revision 3739 by skuang, Thu Jul 14 19:49:12 2011 UTC vs.
Revision 3767 by gezelter, Fri Sep 30 19:37:13 2011 UTC

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

Diff Legend

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