[ViewVC] Diff of: group/interfacial/interfacial.tex

Comparing interfacial/interfacial.tex (file contents):
Revision 3731 by skuang, Tue Jul 5 21:30:29 2011 UTC vs.
Revision 3749 by skuang, Mon Jul 25 19:11:33 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
#	Line 45 \| Line 45 \| We have developed a Non-Isotropic Velocity Scaling alg
45
46		\begin{abstract}
47
48	<	We have developed a Non-Isotropic Velocity Scaling algorithm for
49	<	setting up and maintaining stable thermal gradients in non-equilibrium
50	<	molecular dynamics simulations. This approach effectively imposes
51	<	unphysical thermal flux even between particles of different
52	<	identities, conserves linear momentum and kinetic energy, and
53	<	minimally perturbs the velocity profile of a system when compared with
54	<	previous RNEMD methods. We have used this method to simulate thermal
55	<	conductance at metal / organic solvent interfaces both with and
56	<	without the presence of thiol-based capping agents. We obtained
57	<	values comparable with experimental values, and observed significant
58	<	conductance enhancement with the presence of capping agents. Computed
59	<	power spectra indicate the acoustic impedance mismatch between metal
60	<	and liquid phase is greatly reduced by the capping agents and thus
61	<	leads to higher interfacial thermal transfer efficiency.
48	>	With the Non-Isotropic Velocity Scaling algorithm (NIVS) we have
49	>	developed, an unphysical thermal flux can be effectively set up even
50	>	for non-homogeneous systems like interfaces in non-equilibrium
51	>	molecular dynamics simulations. In this work, this algorithm is
52	>	applied for simulating thermal conductance at metal / organic solvent
53	>	interfaces with various coverages of butanethiol capping
54	>	agents. Different solvents and force field models were tested. Our
55	>	results suggest that the United-Atom models are able to provide an
56	>	estimate of the interfacial thermal conductivity comparable to
57	>	experiments in our simulations with satisfactory computational
58	>	efficiency. From our results, the acoustic impedance mismatch between
59	>	metal and liquid phase is effectively reduced by the capping
60	>	agents, and thus leads to interfacial thermal conductance
61	>	enhancement. Furthermore, this effect is closely related to the
62	>	capping agent coverage on the metal surfaces and the type of solvent
63	>	molecules, and is affected by the models used in the simulations.
64
65		\end{abstract}
66
#	Line 71 \| Line 73 \| leads to higher interfacial thermal transfer efficienc
73		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74
75		\section{Introduction}
74	–	[BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM]
76		Interfacial thermal conductance is extensively studied both
77	<	experimentally and computationally, and systems with interfaces
78	<	present are generally heterogeneous. Although interfaces are commonly
79	<	barriers to heat transfer, it has been
80	<	reported\cite{doi:10.1021/la904855s} that under specific circustances,
81	<	e.g. with certain capping agents present on the surface, interfacial
82	<	conductance can be significantly enhanced. However, heat conductance
83	<	of molecular and nano-scale interfaces will be affected by the
84	<	chemical details of the surface and is challenging to
85	<	experimentalist. The lower thermal flux through interfaces is even
85	<	more difficult to measure with EMD and forward NEMD simulation
86	<	methods. Therefore, developing good simulation methods will be
87	<	desirable in order to investigate thermal transport across interfaces.
77	>	experimentally and computationally\cite{cahill:793}, due to its
78	>	importance in nanoscale science and technology. Reliability of
79	>	nanoscale devices depends on their thermal transport
80	>	properties. Unlike bulk homogeneous materials, nanoscale materials
81	>	features significant presence of interfaces, and these interfaces
82	>	could dominate the heat transfer behavior of these
83	>	materials. Furthermore, these materials are generally heterogeneous,
84	>	which challenges traditional research methods for homogeneous
85	>	systems.
86
87	+	Heat conductance of molecular and nano-scale interfaces will be
88	+	affected by the chemical details of the surface. Experimentally,
89	+	various interfaces have been investigated for their thermal
90	+	conductance properties. Wang {\it et al.} studied heat transport
91	+	through long-chain hydrocarbon monolayers on gold substrate at
92	+	individual molecular level\cite{Wang10082007}; Schmidt {\it et al.}
93	+	studied the role of CTAB on thermal transport between gold nanorods
94	+	and solvent\cite{doi:10.1021/jp8051888}; Juv\'e {\it et al.} studied
95	+	the cooling dynamics, which is controlled by thermal interface
96	+	resistence of glass-embedded metal
97	+	nanoparticles\cite{PhysRevB.80.195406}. Although interfaces are
98	+	commonly barriers for heat transport, Alper {\it et al.} suggested
99	+	that specific ligands (capping agents) could completely eliminate this
100	+	barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}.
101	+
102	+	Theoretical and computational models have also been used to study the
103	+	interfacial thermal transport in order to gain an understanding of
104	+	this phenomena at the molecular level. Recently, Hase and coworkers
105	+	employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
106	+	study thermal transport from hot Au(111) substrate to a self-assembled
107	+	monolayer of alkylthiol with relatively long chain (8-20 carbon
108	+	atoms)\cite{hase:2010,hase:2011}. However, ensemble averaged
109	+	measurements for heat conductance of interfaces between the capping
110	+	monolayer on Au and a solvent phase has yet to be studied.
111	+	The comparatively low thermal flux through interfaces is
112	+	difficult to measure with Equilibrium MD or forward NEMD simulation
113	+	methods. Therefore, the Reverse NEMD (RNEMD) methods would have the
114	+	advantage of having this difficult to measure flux known when studying
115	+	the thermal transport across interfaces, given that the simulation
116	+	methods being able to effectively apply an unphysical flux in
117	+	non-homogeneous systems.
118	+
119		Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
120		algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
121		retains the desirable features of RNEMD (conservation of linear
122		momentum and total energy, compatibility with periodic boundary
123		conditions) while establishing true thermal distributions in each of
124	<	the two slabs. Furthermore, it allows more effective thermal exchange
125	<	between particles of different identities, and thus enables extensive
126	<	study of interfacial conductance.
124	>	the two slabs. Furthermore, it allows effective thermal exchange
125	>	between particles of different identities, and thus makes the study of
126	>	interfacial conductance much simpler.
127
128	+	The work presented here deals with the Au(111) surface covered to
129	+	varying degrees by butanethiol, a capping agent with short carbon
130	+	chain, and solvated with organic solvents of different molecular
131	+	properties. Different models were used for both the capping agent and
132	+	the solvent force field parameters. Using the NIVS algorithm, the
133	+	thermal transport across these interfaces was studied and the
134	+	underlying mechanism for the phenomena was investigated.
135	+
136	+	[MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]
137	+
138		\section{Methodology}
139	<	\subsection{Algorithm}
140	<	[BACKGROUND FOR MD METHODS]
141	<	There have been many algorithms for computing thermal conductivity
142	<	using molecular dynamics simulations. However, interfacial conductance
143	<	is at least an order of magnitude smaller. This would make the
144	<	calculation even more difficult for those slowly-converging
145	<	equilibrium methods. Imposed-flux non-equilibrium
139	>	\subsection{Imposd-Flux Methods in MD Simulations}
140	>	Steady state MD simulations has the advantage that not many
141	>	trajectories are needed to study the relationship between thermal flux
142	>	and thermal gradients. For systems including low conductance
143	>	interfaces one must have a method capable of generating or measuring
144	>	relatively small fluxes, compared to those required for bulk
145	>	conductivity. This requirement makes the calculation even more
146	>	difficult for those slowly-converging equilibrium
147	>	methods\cite{Viscardy:2007lq}. Forward methods may impose gradient,
148	>	but in interfacial conditions it is not clear what behavior to impose
149	>	at the interfacial boundaries. Imposed-flux reverse non-equilibrium
150		methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
151	<	the response of temperature or momentum gradients are easier to
152	<	measure than the flux, if unknown, and thus, is a preferable way to
153	<	the forward NEMD methods. Although the momentum swapping approach for
154	<	flux-imposing can be used for exchanging energy between particles of
155	<	different identity, the kinetic energy transfer efficiency is affected
156	<	by the mass difference between the particles, which limits its
113	<	application on heterogeneous interfacial systems.
151	>	the thermal response becomes easier to measure than the flux. Although
152	>	M\"{u}ller-Plathe's original momentum swapping approach can be used
153	>	for exchanging energy between particles of different identity, the
154	>	kinetic energy transfer efficiency is affected by the mass difference
155	>	between the particles, which limits its application on heterogeneous
156	>	interfacial systems.
157
158	<	The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach in
159	<	non-equilibrium MD simulations is able to impose relatively large
160	<	kinetic energy flux without obvious perturbation to the velocity
161	<	distribution of the simulated systems. Furthermore, this approach has
158	>	The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach to
159	>	non-equilibrium MD simulations is able to impose a wide range of
160	>	kinetic energy fluxes without obvious perturbation to the velocity
161	>	distributions of the simulated systems. Furthermore, this approach has
162		the advantage in heterogeneous interfaces in that kinetic energy flux
163		can be applied between regions of particles of arbitary identity, and
164	<	the flux quantity is not restricted by particle mass difference.
164	>	the flux will not be restricted by difference in particle mass.
165
166		The NIVS algorithm scales the velocity vectors in two separate regions
167		of a simulation system with respective diagonal scaling matricies. To
168		determine these scaling factors in the matricies, a set of equations
169		including linear momentum conservation and kinetic energy conservation
170	<	constraints and target momentum/energy flux satisfaction is
171	<	solved. With the scaling operation applied to the system in a set
172	<	frequency, corresponding momentum/temperature gradients can be built,
173	<	which can be used for computing transportation properties and other
174	<	applications related to momentum/temperature gradients. The NIVS
132	<	algorithm conserves momenta and energy and does not depend on an
133	<	external thermostat.
170	>	constraints and target energy flux satisfaction is solved. With the
171	>	scaling operation applied to the system in a set frequency, bulk
172	>	temperature gradients can be easily established, and these can be used
173	>	for computing thermal conductivities. The NIVS algorithm conserves
174	>	momenta and energy and does not depend on an external thermostat.
175
176		\subsection{Defining Interfacial Thermal Conductivity $G$}
177	<	For interfaces with a relatively low interfacial conductance, the bulk
178	<	regions on either side of an interface rapidly come to a state in
179	<	which the two phases have relatively homogeneous (but distinct)
180	<	temperatures. The interfacial thermal conductivity $G$ can therefore
181	<	be approximated as:
177	>	Given a system with thermal gradients and the corresponding thermal
178	>	flux, for interfaces with a relatively low interfacial conductance,
179	>	the bulk regions on either side of an interface rapidly come to a
180	>	state in which the two phases have relatively homogeneous (but
181	>	distinct) temperatures. The interfacial thermal conductivity $G$ can
182	>	therefore be approximated as:
183		\begin{equation}
184		G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
185		\langle T_\mathrm{cold}\rangle \right)}
#	Line 148 \| Line 190 \| When the interfacial conductance is {\it not} small, t
190		T_\mathrm{cold}\rangle}$ are the average observed temperature of the
191		two separated phases.
192
193	<	When the interfacial conductance is {\it not} small, two ways can be
194	<	used to define $G$.
193	>	When the interfacial conductance is {\it not} small, there are two
194	>	ways to define $G$.
195
196	<	One way is to assume the temperature is discretely different on two
197	<	sides of the interface, $G$ can be calculated with the thermal flux
198	<	applied $J$ and the maximum temperature difference measured along the
199	<	thermal gradient max($\Delta T$), which occurs at the interface, as:
196	>	One way is to assume the temperature is discrete on the two sides of
197	>	the interface. $G$ can be calculated using the applied thermal flux
198	>	$J$ and the maximum temperature difference measured along the thermal
199	>	gradient max($\Delta T$), which occurs at the Gibbs deviding surface
200	>	(Figure \ref{demoPic}):
201		\begin{equation}
202		G=\frac{J}{\Delta T}
203		\label{discreteG}
204		\end{equation}
205
206	+	\begin{figure}
207	+	\includegraphics[width=\linewidth]{method}
208	+	\caption{Interfacial conductance can be calculated by applying an
209	+	(unphysical) kinetic energy flux between two slabs, one located
210	+	within the metal and another on the edge of the periodic box. The
211	+	system responds by forming a thermal response or a gradient. In
212	+	bulk liquids, this gradient typically has a single slope, but in
213	+	interfacial systems, there are distinct thermal conductivity
214	+	domains. The interfacial conductance, $G$ is found by measuring the
215	+	temperature gap at the Gibbs dividing surface, or by using second
216	+	derivatives of the thermal profile.}
217	+	\label{demoPic}
218	+	\end{figure}
219	+
220		The other approach is to assume a continuous temperature profile along
221		the thermal gradient axis (e.g. $z$) and define $G$ at the point where
222		the magnitude of thermal conductivity $\lambda$ change reach its
#	Line 175 \| Line 232 \| difference method and thus calculate $G^\prime$.
232
233		With the temperature profile obtained from simulations, one is able to
234		approximate the first and second derivatives of $T$ with finite
235	<	difference method and thus calculate $G^\prime$.
235	>	difference methods and thus calculate $G^\prime$.
236
237	<	In what follows, both definitions are used for calculation and comparison.
237	>	In what follows, both definitions have been used for calculation and
238	>	are compared in the results.
239
240	<	[IMPOSE G DEFINITION INTO OUR SYSTEMS]
241	<	To facilitate the use of the above definitions in calculating $G$ and
242	<	$G^\prime$, we have a metal slab with its (111) surfaces perpendicular
243	<	to the $z$-axis of our simulation cells. With or withour capping
244	<	agents on the surfaces, the metal slab is solvated with organic
187	<	solvents, as illustrated in Figure \ref{demoPic}.
240	>	To compare the above definitions ($G$ and $G^\prime$), we have modeled
241	>	a metal slab with its (111) surfaces perpendicular to the $z$-axis of
242	>	our simulation cells. Both with and without capping agents on the
243	>	surfaces, the metal slab is solvated with simple organic solvents, as
244	>	illustrated in Figure \ref{gradT}.
245
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.}
250	<	\label{demoPic}
251	<	\end{figure}
246	>	With the simulation cell described above, we are able to equilibrate
247	>	the system and impose an unphysical thermal flux between the liquid
248	>	and the metal phase using the NIVS algorithm. By periodically applying
249	>	the unphysical flux, we are able to obtain a temperature profile and
250	>	its spatial derivatives. These quantities enable the evaluation of the
251	>	interfacial thermal conductance of a surface. Figure \ref{gradT} is an
252	>	example of how an applied thermal flux can be used to obtain the 1st
253	>	and 2nd derivatives of the temperature profile.
254
196	–	With a simulation cell setup following the above manner, one is able
197	–	to equilibrate the system and impose an unphysical thermal flux
198	–	between the liquid and the metal phase with the NIVS algorithm. Under
199	–	a stablized thermal gradient induced by periodically applying the
200	–	unphysical flux, one is able to obtain a temperature profile and the
201	–	physical thermal flux corresponding to it, which equals to the
202	–	unphysical flux applied by NIVS. These data enables the evaluation of
203	–	the interfacial thermal conductance of a surface. Figure \ref{gradT}
204	–	is an example how those stablized thermal gradient can be used to
205	–	obtain the 1st and 2nd derivatives of the temperature profile.
206	–
255		\begin{figure}
256		\includegraphics[width=\linewidth]{gradT}
257	<	\caption{The 1st and 2nd derivatives of temperature profile can be
258	<	obtained with finite difference approximation.}
257	>	\caption{A sample of Au-butanethiol/hexane interfacial system and the
258	>	temperature profile after a kinetic energy flux is imposed to
259	>	it. The 1st and 2nd derivatives of the temperature profile can be
260	>	obtained with finite difference approximation (lower panel).}
261		\label{gradT}
262		\end{figure}
263
214	–	[MAY INCLUDE POWER SPECTRUM PROTOCOL]
215	–
264		\section{Computational Details}
265		\subsection{Simulation Protocol}
266	<	In our simulations, Au is used to construct a metal slab with bare
267	<	(111) surface perpendicular to the $z$-axis. Different slab thickness
268	<	(layer numbers of Au) are simulated. This metal slab is first
269	<	equilibrated under normal pressure (1 atm) and a desired
270	<	temperature. After equilibration, butanethiol is used as the capping
271	<	agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
272	<	atoms in the butanethiol molecules would occupy the three-fold sites
273	<	of the surfaces, and the maximal butanethiol capacity on Au surface is
274	<	$1/3$ of the total number of surface Au atoms[CITATION]. A series of
275	<	different coverage surfaces is investigated in order to study the
276	<	relation between coverage and conductance.
266	>	The NIVS algorithm has been implemented in our MD simulation code,
267	>	OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
268	>	simulations. Different metal slab thickness (layer numbers of Au) was
269	>	simulated. Metal slabs were first equilibrated under atmospheric
270	>	pressure (1 atm) and a desired temperature (e.g. 200K). After
271	>	equilibration, butanethiol capping agents were placed at three-fold
272	>	hollow sites on the Au(111) surfaces. These sites could be either a
273	>	{\it fcc} or {\it hcp} site. However, Hase {\it et al.} found that
274	>	they are equivalent in a heat transfer process\cite{hase:2010}, so
275	>	they are not distinguished in our study. The maximum butanethiol
276	>	capacity on Au surface is $1/3$ of the total number of surface Au
277	>	atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
278	>	structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
279	>	series of different coverages was derived by evenly eliminating
280	>	butanethiols on the surfaces, and was investigated in order to study
281	>	the relation between coverage and interfacial conductance.
282
283	<	[COVERAGE DISCRIPTION] However, since the interactions between surface
284	<	Au and butanethiol is non-bonded, the capping agent molecules are
285	<	allowed to migrate to an empty neighbor three-fold site during a
286	<	simulation. Therefore, the initial configuration would not severely
287	<	affect the sampling of a variety of configurations of the same
288	<	coverage, and the final conductance measurement would be an average
289	<	effect of these configurations explored in the simulations. [MAY NEED FIGURES]
283	>	The capping agent molecules were allowed to migrate during the
284	>	simulations. They distributed themselves uniformly and sampled a
285	>	number of three-fold sites throughout out study. Therefore, the
286	>	initial configuration would not noticeably affect the sampling of a
287	>	variety of configurations of the same coverage, and the final
288	>	conductance measurement would be an average effect of these
289	>	configurations explored in the simulations. [MAY NEED SNAPSHOTS]
290
291	<	After the modified Au-butanethiol surface systems are equilibrated
292	<	under canonical ensemble, Packmol\cite{packmol} is used to pack
293	<	organic solvent molecules in the previously vacuum part of the
294	<	simulation cells, which guarantees that short range repulsive
295	<	interactions do not disrupt the simulations. Two solvents are
296	<	investigated, one which has little vibrational overlap with the
244	<	alkanethiol and plane-like shape (toluene), and one which has similar
245	<	vibrational frequencies and chain-like shape ({\it n}-hexane). [MAY
246	<	EXPLAIN WHY WE CHOOSE THEM]
291	>	After the modified Au-butanethiol surface systems were equilibrated
292	>	under canonical ensemble, organic solvent molecules were packed in the
293	>	previously empty part of the simulation cells\cite{packmol}. Two
294	>	solvents were investigated, one which has little vibrational overlap
295	>	with the alkanethiol and a planar shape (toluene), and one which has
296	>	similar vibrational frequencies and chain-like shape ({\it n}-hexane).
297
298	<	The spacing filled by solvent molecules, i.e. the gap between
298	>	The space filled by solvent molecules, i.e. the gap between
299		periodically repeated Au-butanethiol surfaces should be carefully
300		chosen. A very long length scale for the thermal gradient axis ($z$)
301		may cause excessively hot or cold temperatures in the middle of the
#	Line 254 \| Line 304 \| corresponding spacing is usually $35 \sim 60$\AA.
304		solvent molecules would change the normal behavior of the liquid
305		phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
306		these extreme cases did not happen to our simulations. And the
307	<	corresponding spacing is usually $35 \sim 60$\AA.
307	>	corresponding spacing is usually $35 \sim 75$\AA.
308
309	<	The initial configurations generated by Packmol are further
310	<	equilibrated with the $x$ and $y$ dimensions fixed, only allowing
311	<	length scale change in $z$ dimension. This is to ensure that the
312	<	equilibration of liquid phase does not affect the metal crystal
313	<	structure in $x$ and $y$ dimensions. Further equilibration are run
314	<	under NVT and then NVE ensembles.
309	>	The initial configurations generated are further equilibrated with the
310	>	$x$ and $y$ dimensions fixed, only allowing length scale change in $z$
311	>	dimension. This is to ensure that the equilibration of liquid phase
312	>	does not affect the metal crystal structure in $x$ and $y$ dimensions.
313	>	To investigate this effect, comparisons were made with simulations
314	>	that allow changes of $L_x$ and $L_y$ during NPT equilibration, and
315	>	the results are shown in later sections. After ensuring the liquid
316	>	phase reaches equilibrium at atmospheric pressure (1 atm), further
317	>	equilibration are followed under NVT and then NVE ensembles.
318
319		After the systems reach equilibrium, NIVS is implemented to impose a
320		periodic unphysical thermal flux between the metal and the liquid
321		phase. Most of our simulations are under an average temperature of
322		$\sim$200K. Therefore, this flux usually comes from the metal to the
323		liquid so that the liquid has a higher temperature and would not
324	<	freeze due to excessively low temperature. This induced temperature
325	<	gradient is stablized and the simulation cell is devided evenly into
326	<	N slabs along the $z$-axis and the temperatures of each slab are
327	<	recorded. When the slab width $d$ of each slab is the same, the
328	<	derivatives of $T$ with respect to slab number $n$ can be directly
329	<	used for $G^\prime$ calculations:
324	>	freeze due to excessively low temperature. After this induced
325	>	temperature gradient is stablized, the temperature profile of the
326	>	simulation cell is recorded. To do this, the simulation cell is
327	>	devided evenly into $N$ slabs along the $z$-axis and $N$ is maximized
328	>	for highest possible spatial resolution but not too many to have some
329	>	slabs empty most of the time. The average temperatures of each slab
330	>	are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
331	>	the same, the derivatives of $T$ with respect to slab number $n$ can
332	>	be directly used for $G^\prime$ calculations:
333		\begin{equation}
334		G^\prime = \|J_z\|\Big\|\frac{\partial^2 T}{\partial z^2}\Big\|
335		\Big/\left(\frac{\partial T}{\partial z}\right)^2
#	Line 284 \| Line 340 \| G^\prime = \|J_z\|\Big\|\frac{\partial^2 T}{\partial z^2}
340		\label{derivativeG2}
341		\end{equation}
342
343	+	All of the above simulation procedures use a time step of 1 fs. And
344	+	each equilibration / stabilization step usually takes 100 ps, or
345	+	longer, if necessary.
346	+
347		\subsection{Force Field Parameters}
348	<	Our simulations include various components. Therefore, force field
349	<	parameter descriptions are needed for interactions both between the
350	<	same type of particles and between particles of different species.
348	>	Our simulations include various components. Figure \ref{demoMol}
349	>	demonstrates the sites defined for both United-Atom and All-Atom
350	>	models of the organic solvent and capping agent molecules in our
351	>	simulations. Force field parameter descriptions are needed for
352	>	interactions both between the same type of particles and between
353	>	particles of different species.
354
355	+	\begin{figure}
356	+	\includegraphics[width=\linewidth]{structures}
357	+	\caption{Structures of the capping agent and solvents utilized in
358	+	these simulations. The chemically-distinct sites (a-e) are expanded
359	+	in terms of constituent atoms for both United Atom (UA) and All Atom
360	+	(AA) force fields. Most parameters are from
361	+	Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} (UA) and
362	+	\protect\cite{OPLSAA} (AA). Cross-interactions with the Au atoms are given
363	+	in Table \ref{MnM}.}
364	+	\label{demoMol}
365	+	\end{figure}
366	+
367		The Au-Au interactions in metal lattice slab is described by the
368	<	quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
368	>	quantum Sutton-Chen (QSC) formulation\cite{PhysRevB.59.3527}. The QSC
369		potentials include zero-point quantum corrections and are
370		reparametrized for accurate surface energies compared to the
371		Sutton-Chen potentials\cite{Chen90}.
372
298	–	Figure [REF] demonstrates how we name our pseudo-atoms of the
299	–	molecules in our simulations.
300	–	[FIGURE FOR MOLECULE NOMENCLATURE]
301	–
373		For both solvent molecules, straight chain {\it n}-hexane and aromatic
374		toluene, United-Atom (UA) and All-Atom (AA) models are used
375		respectively. The TraPPE-UA
376		parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
377	<	for our UA solvent molecules. In these models, pseudo-atoms are
378	<	located at the carbon centers for alkyl groups. By eliminating
379	<	explicit hydrogen atoms, these models are simple and computationally
380	<	efficient, while maintains good accuracy. However, the TraPPE-UA for
381	<	alkanes is known to predict a lower boiling point than experimental
311	<	values. Considering that after an unphysical thermal flux is applied
312	<	to a system, the temperature of ``hot'' area in the liquid phase would be
313	<	significantly higher than the average, to prevent over heating and
314	<	boiling of the liquid phase, the average temperature in our
315	<	simulations should be much lower than the liquid boiling point. [MORE DISCUSSION]
316	<	For UA-toluene model, rigid body constraints are applied, so that the
317	<	benzene ring and the methyl-CRar bond are kept rigid. This would save
318	<	computational time.[MORE DETAILS]
377	>	for our UA solvent molecules. In these models, sites are located at
378	>	the carbon centers for alkyl groups. Bonding interactions, including
379	>	bond stretches and bends and torsions, were used for intra-molecular
380	>	sites not separated by more than 3 bonds. Otherwise, for non-bonded
381	>	interactions, Lennard-Jones potentials are used. [CHECK CITATION]
382
383	+	By eliminating explicit hydrogen atoms, these models are simple and
384	+	computationally efficient, while maintains good accuracy. However, the
385	+	TraPPE-UA for alkanes is known to predict a lower boiling point than
386	+	experimental values. Considering that after an unphysical thermal flux
387	+	is applied to a system, the temperature of ``hot'' area in the liquid
388	+	phase would be significantly higher than the average of the system, to
389	+	prevent over heating and boiling of the liquid phase, the average
390	+	temperature in our simulations should be much lower than the liquid
391	+	boiling point.
392	+
393	+	For UA-toluene model, the non-bonded potentials between
394	+	inter-molecular sites have a similar Lennard-Jones formulation. For
395	+	intra-molecular interactions, considering the stiffness of the benzene
396	+	ring, rigid body constraints are applied for further computational
397	+	efficiency. All bonds in the benzene ring and between the ring and the
398	+	methyl group remain rigid during the progress of simulations.
399	+
400		Besides the TraPPE-UA models, AA models for both organic solvents are
401		included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
402	<	force field is used. [MORE DETAILS]
403	<	For toluene, the United Force Field developed by Rapp\'{e} {\it et
404	<	al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]
402	>	force field is used. Additional explicit hydrogen sites were
403	>	included. Besides bonding and non-bonded site-site interactions,
404	>	partial charges and the electrostatic interactions were added to each
405	>	CT and HC site. For toluene, the United Force Field developed by
406	>	Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} is
407	>	adopted. Without the rigid body constraints, bonding interactions were
408	>	included. For the aromatic ring, improper torsions (inversions) were
409	>	added as an extra potential for maintaining the planar shape.
410	>	[CHECK CITATION]
411
412		The capping agent in our simulations, the butanethiol molecules can
413		either use UA or AA model. The TraPPE-UA force fields includes
#	Line 330 \| Line 416 \| Au(111) surfaces, we adopt the S parameters from [CITA
416		parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
417		surfaces do not have the hydrogen atom bonded to sulfur. To adapt this
418		change and derive suitable parameters for butanethiol adsorbed on
419	<	Au(111) surfaces, we adopt the S parameters from [CITATION CF VLUGT]
420	<	and modify parameters for its neighbor C atom for charge balance in
421	<	the molecule. Note that the model choice (UA or AA) of capping agent
422	<	can be different from the solvent. Regardless of model choice, the
423	<	force field parameters for interactions between capping agent and
424	<	solvent can be derived using Lorentz-Berthelot Mixing Rule:
419	>	Au(111) surfaces, we adopt the S parameters from Luedtke and
420	>	Landman\cite{landman:1998}[CHECK CITATION]
421	>	and modify parameters for its neighbor C
422	>	atom for charge balance in the molecule. Note that the model choice
423	>	(UA or AA) of capping agent can be different from the
424	>	solvent. Regardless of model choice, the force field parameters for
425	>	interactions between capping agent and solvent can be derived using
426	>	Lorentz-Berthelot Mixing Rule:
427	>	\begin{eqnarray}
428	>	\sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
429	>	\epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
430	>	\end{eqnarray}
431
340	–
432		To describe the interactions between metal Au and non-metal capping
433		agent and solvent particles, we refer to an adsorption study of alkyl
434		thiols on gold surfaces by Vlugt {\it et
435		al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones
436		form of potential parameters for the interaction between Au and
437		pseudo-atoms CH$_x$ and S based on a well-established and widely-used
438	<	effective potential of Hautman and Klein[CITATION] for the Au(111)
439	<	surface. As our simulations require the gold lattice slab to be
440	<	non-rigid so that it could accommodate kinetic energy for thermal
438	>	effective potential of Hautman and Klein\cite{hautman:4994} for the
439	>	Au(111) surface. As our simulations require the gold lattice slab to
440	>	be non-rigid so that it could accommodate kinetic energy for thermal
441		transport study purpose, the pair-wise form of potentials is
442		preferred.
443
444		Besides, the potentials developed from {\it ab initio} calculations by
445		Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
446	<	interactions between Au and aromatic C/H atoms in toluene.[MORE DETAILS]
446	>	interactions between Au and aromatic C/H atoms in toluene. A set of
447	>	pseudo Lennard-Jones parameters were provided for Au in their force
448	>	fields. By using the Mixing Rule, this can be used to derive pair-wise
449	>	potentials for non-bonded interactions between Au and non-metal sites.
450
451		However, the Lennard-Jones parameters between Au and other types of
452	<	particles in our simulations are not yet well-established. For these
453	<	interactions, we attempt to derive their parameters using the Mixing
454	<	Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters
455	<	for Au is first extracted from the Au-CH$_x$ parameters by applying
456	<	the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
452	>	particles, such as All-Atom normal alkanes in our simulations are not
453	>	yet well-established. For these interactions, we attempt to derive
454	>	their parameters using the Mixing Rule. To do this, Au pseudo
455	>	Lennard-Jones parameters for ``Metal-non-Metal'' (MnM) interactions
456	>	were first extracted from the Au-CH$_x$ parameters by applying the
457	>	Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
458		parameters in our simulations.
459
460		\begin{table*}
461		\begin{minipage}{\linewidth}
462		\begin{center}
463	<	\caption{Lennard-Jones parameters for Au-non-Metal
464	<	interactions in our simulations.}
465	<
466	<	\begin{tabular}{ccc}
463	>	\caption{Non-bonded interaction parameters (including cross
464	>	interactions with Au atoms) for both force fields used in this
465	>	work.}
466	>	\begin{tabular}{lllllll}
467		\hline\hline
468	<	Non-metal & $\sigma$/\AA & $\epsilon$/kcal/mol \\
468	>	& Site & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
469	>	$\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
470	>	& & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
471		\hline
472	<	S & 2.40 & 8.465 \\
473	<	CH3 & 3.54 & 0.2146 \\
474	<	CH2 & 3.54 & 0.1749 \\
475	<	CT3 & 3.365 & 0.1373 \\
476	<	CT2 & 3.365 & 0.1373 \\
477	<	CTT & 3.365 & 0.1373 \\
478	<	HC & 2.865 & 0.09256 \\
479	<	CHar & 3.4625 & 0.1680 \\
480	<	CRar & 3.555 & 0.1604 \\
481	<	CA & 3.173 & 0.0640 \\
482	<	HA & 2.746 & 0.0414 \\
472	>	United Atom (UA)
473	>	&CH3 & 3.75 & 0.1947 & - & 3.54 & 0.2146 \\
474	>	&CH2 & 3.95 & 0.0914 & - & 3.54 & 0.1749 \\
475	>	&CHar & 3.695 & 0.1003 & - & 3.4625 & 0.1680 \\
476	>	&CRar & 3.88 & 0.04173 & - & 3.555 & 0.1604 \\
477	>	\hline
478	>	All Atom (AA)
479	>	&CT3 & 3.50 & 0.066 & -0.18 & 3.365 & 0.1373 \\
480	>	&CT2 & 3.50 & 0.066 & -0.12 & 3.365 & 0.1373 \\
481	>	&CTT & 3.50 & 0.066 & -0.065 & 3.365 & 0.1373 \\
482	>	&HC & 2.50 & 0.030 & 0.06 & 2.865 & 0.09256 \\
483	>	&CA & 3.55 & 0.070 & -0.115 & 3.173 & 0.0640 \\
484	>	&HA & 2.42 & 0.030 & 0.115 & 2.746 & 0.0414 \\
485	>	\hline
486	>	Both UA and AA
487	>	& S & 4.45 & 0.25 & - & 2.40 & 8.465 \\
488		\hline\hline
489		\end{tabular}
490		\label{MnM}
#	Line 390 \| Line 492 \| parameters in our simulations.
492		\end{minipage}
493		\end{table*}
494
495	+	\subsection{Vibrational Spectrum}
496	+	To investigate the mechanism of interfacial thermal conductance, the
497	+	vibrational spectrum is utilized as a complementary tool. Vibrational
498	+	spectra were taken for individual components in different
499	+	simulations. To obtain these spectra, simulations were run after
500	+	equilibration, in the NVE ensemble. Snapshots of configurations were
501	+	collected at a frequency that is higher than that of the fastest
502	+	vibrations occuring in the simulations. With these configurations, the
503	+	velocity auto-correlation functions can be computed:
504	+	\begin{equation}
505	+	C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
506	+	\label{vCorr}
507	+	\end{equation}
508	+	Followed by Fourier transforms, the power spectrum can be constructed:
509	+	\begin{equation}
510	+	\hat{f}(\omega) = \int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
511	+	\label{fourier}
512	+	\end{equation}
513
514		\section{Results and Discussions}
515	<	[MAY HAVE A BRIEF SUMMARY]
515	>	In what follows, how the parameters and protocol of simulations would
516	>	affect the measurement of $G$'s is first discussed. With a reliable
517	>	protocol and set of parameters, the influence of capping agent
518	>	coverage on thermal conductance is investigated. Besides, different
519	>	force field models for both solvents and selected deuterated models
520	>	were tested and compared. Finally, a summary of the role of capping
521	>	agent in the interfacial thermal transport process is given.
522	>
523		\subsection{How Simulation Parameters Affects $G$}
397	–	[MAY NOT PUT AT FIRST]
524		We have varied our protocol or other parameters of the simulations in
525		order to investigate how these factors would affect the measurement of
526		$G$'s. It turned out that while some of these parameters would not
#	Line 403 \| Line 529 \| during equilibrating the liquid phase. Due to the stif
529		results.
530
531		In some of our simulations, we allowed $L_x$ and $L_y$ to change
532	<	during equilibrating the liquid phase. Due to the stiffness of the Au
533	<	slab, $L_x$ and $L_y$ would not change noticeably after
534	<	equilibration. Although $L_z$ could fluctuate $\sim$1\% after a system
535	<	is fully equilibrated in the NPT ensemble, this fluctuation, as well
536	<	as those comparably smaller to $L_x$ and $L_y$, would not be magnified
537	<	on the calculated $G$'s, as shown in Table \ref{AuThiolHexaneUA}. This
538	<	insensivity to $L_i$ fluctuations allows reliable measurement of $G$'s
539	<	without the necessity of extremely cautious equilibration process.
532	>	during equilibrating the liquid phase. Due to the stiffness of the
533	>	crystalline Au structure, $L_x$ and $L_y$ would not change noticeably
534	>	after equilibration. Although $L_z$ could fluctuate $\sim$1\% after a
535	>	system is fully equilibrated in the NPT ensemble, this fluctuation, as
536	>	well as those of $L_x$ and $L_y$ (which is significantly smaller),
537	>	would not be magnified on the calculated $G$'s, as shown in Table
538	>	\ref{AuThiolHexaneUA}. This insensivity to $L_i$ fluctuations allows
539	>	reliable measurement of $G$'s without the necessity of extremely
540	>	cautious equilibration process.
541
542		As stated in our computational details, the spacing filled with
543		solvent molecules can be chosen within a range. This allows some
#	Line 437 \| Line 564 \| $J_z$ is listed in Table \ref{AuThiolHexaneUA} and [RE
564		the thermal flux across the interface. For our simulations, we denote
565		$J_z$ to be positive when the physical thermal flux is from the liquid
566		to metal, and negative vice versa. The $G$'s measured under different
567	<	$J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These
568	<	results do not suggest that $G$ is dependent on $J_z$ within this flux
569	<	range. The linear response of flux to thermal gradient simplifies our
570	<	investigations in that we can rely on $G$ measurement with only a
571	<	couple $J_z$'s and do not need to test a large series of fluxes.
567	>	$J_z$ is listed in Table \ref{AuThiolHexaneUA} and
568	>	\ref{AuThiolToluene}. These results do not suggest that $G$ is
569	>	dependent on $J_z$ within this flux range. The linear response of flux
570	>	to thermal gradient simplifies our investigations in that we can rely
571	>	on $G$ measurement with only a couple $J_z$'s and do not need to test
572	>	a large series of fluxes.
573
446	–	%ADD MORE TO TABLE
574		\begin{table*}
575		\begin{minipage}{\linewidth}
576		\begin{center}
577		\caption{Computed interfacial thermal conductivity ($G$ and
578		$G^\prime$) values for the 100\% covered Au-butanethiol/hexane
579		interfaces with UA model and different hexane molecule numbers
580	<	at different temperatures using a range of energy fluxes.}
580	>	at different temperatures using a range of energy
581	>	fluxes. Error estimates indicated in parenthesis.}
582
583	<	\begin{tabular}{cccccccc}
583	>	\begin{tabular}{ccccccc}
584		\hline\hline
585	<	$\langle T\rangle$ & & $L_x$ & $L_y$ & $L_z$ & $J_z$ &
586	<	$G$ & $G^\prime$ \\
587	<	(K) & $N_{hexane}$ & \multicolumn{3}{c}\AA & (GW/m$^2$) &
585	>	$\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
586	>	$J_z$ & $G$ & $G^\prime$ \\
587	>	(K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
588		\multicolumn{2}{c}{(MW/m$^2$/K)} \\
589		\hline
590	<	200 & 266 & 29.86 & 25.80 & 113.1 & -0.96 &
591	<	102() & 80.0() \\
592	<	& 200 & 29.84 & 25.81 & 93.9 & 1.92 &
593	<	129() & 87.3() \\
594	<	& & 29.84 & 25.81 & 95.3 & 1.93 &
595	<	131() & 77.5() \\
596	<	& 166 & 29.84 & 25.81 & 85.7 & 0.97 &
597	<	115() & 69.3() \\
598	<	& & & & & 1.94 &
599	<	125() & 87.1() \\
600	<	250 & 200 & 29.84 & 25.87 & 106.8 & 0.96 &
601	<	81.8() & 67.0() \\
602	<	& 166 & 29.87 & 25.84 & 94.8 & 0.98 &
603	<	79.0() & 62.9() \\
604	<	& & 29.84 & 25.85 & 95.0 & 1.44 &
605	<	76.2() & 64.8() \\
590	>	200 & 266 & No & 0.672 & -0.96 & 102(3) & 80.0(0.8) \\
591	>	& 200 & Yes & 0.694 & 1.92 & 129(11) & 87.3(0.3) \\
592	>	& & Yes & 0.672 & 1.93 & 131(16) & 78(13) \\
593	>	& & No & 0.688 & 0.96 & 125(16) & 90.2(15) \\
594	>	& & & & 1.91 & 139(10) & 101(10) \\
595	>	& & & & 2.83 & 141(6) & 89.9(9.8) \\
596	>	& 166 & Yes & 0.679 & 0.97 & 115(19) & 69(18) \\
597	>	& & & & 1.94 & 125(9) & 87.1(0.2) \\
598	>	& & No & 0.681 & 0.97 & 141(30) & 78(22) \\
599	>	& & & & 1.92 & 138(4) & 98.9(9.5) \\
600	>	\hline
601	>	250 & 200 & No & 0.560 & 0.96 & 75(10) & 61.8(7.3) \\
602	>	& & & & -0.95 & 49.4(0.3) & 45.7(2.1) \\
603	>	& 166 & Yes & 0.570 & 0.98 & 79.0(3.5) & 62.9(3.0) \\
604	>	& & No & 0.569 & 0.97 & 80.3(0.6) & 67(11) \\
605	>	& & & & 1.44 & 76.2(5.0) & 64.8(3.8) \\
606	>	& & & & -0.95 & 56.4(2.5) & 54.4(1.1) \\
607	>	& & & & -1.85 & 47.8(1.1) & 53.5(1.5) \\
608		\hline\hline
609		\end{tabular}
610		\label{AuThiolHexaneUA}
#	Line 490 \| Line 620 \| butanethiol as well.[MAY NEED FIGURE] And this reduced
620		temperature is higher than 250K. Additionally, the equilibrated liquid
621		hexane density under 250K becomes lower than experimental value. This
622		expanded liquid phase leads to lower contact between hexane and
623	<	butanethiol as well.[MAY NEED FIGURE] And this reduced contact would
623	>	butanethiol as well.[MAY NEED SLAB DENSITY FIGURE]
624	>	And this reduced contact would
625		probably be accountable for a lower interfacial thermal conductance,
626		as shown in Table \ref{AuThiolHexaneUA}.
627
#	Line 505 \| Line 636 \| in that higher degree of contact could yield increased
636		important role in the thermal transport process across the interface
637		in that higher degree of contact could yield increased conductance.
638
508	–	[ADD Lxyz AND ERROR ESTIMATE TO TABLE]
639		\begin{table*}
640		\begin{minipage}{\linewidth}
641		\begin{center}
642		\caption{Computed interfacial thermal conductivity ($G$ and
643		$G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
644		interface at different temperatures using a range of energy
645	<	fluxes.}
645	>	fluxes. Error estimates indicated in parenthesis.}
646
647	<	\begin{tabular}{cccc}
647	>	\begin{tabular}{ccccc}
648		\hline\hline
649	<	$\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
650	<	(K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
649	>	$\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
650	>	(K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
651		\hline
652	<	200 & -1.86 & 180() & 135() \\
653	<	& 2.15 & 204() & 113() \\
654	<	& -3.93 & 175() & 114() \\
655	<	300 & -1.91 & 143() & 125() \\
656	<	& -4.19 & 134() & 113() \\
652	>	200 & 0.933 & 2.15 & 204(12) & 113(12) \\
653	>	& & -1.86 & 180(3) & 135(21) \\
654	>	& & -3.93 & 176(5) & 113(12) \\
655	>	\hline
656	>	300 & 0.855 & -1.91 & 143(5) & 125(2) \\
657	>	& & -4.19 & 135(9) & 113(12) \\
658		\hline\hline
659		\end{tabular}
660		\label{AuThiolToluene}
#	Line 551 \| Line 682 \| undertaken at $<T>\sim$200K.
682		reconstructions could eliminate the original $x$ and $y$ dimensional
683		homogeneity, measurement of $G$ is more difficult to conduct under
684		higher temperatures. Therefore, most of our measurements are
685	<	undertaken at $<T>\sim$200K.
685	>	undertaken at $\langle T\rangle\sim$200K.
686
687		However, when the surface is not completely covered by butanethiols,
688		the simulated system is more resistent to the reconstruction
689	<	above. Our Au-butanethiol/toluene system did not see this phenomena
690	<	even at $<T>\sim$300K. The Au(111) surfaces have a 90\% coverage of
691	<	butanethiols and have empty three-fold sites. These empty sites could
692	<	help prevent surface reconstruction in that they provide other means
693	<	of capping agent relaxation. It is observed that butanethiols can
694	<	migrate to their neighbor empty sites during a simulation. Therefore,
695	<	we were able to obtain $G$'s for these interfaces even at a relatively
696	<	high temperature without being affected by surface reconstructions.
689	>	above. Our Au-butanethiol/toluene system had the Au(111) surfaces 90\%
690	>	covered by butanethiols, but did not see this above phenomena even at
691	>	$\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
692	>	capping agents could help prevent surface reconstruction in that they
693	>	provide other means of capping agent relaxation. It is observed that
694	>	butanethiols can migrate to their neighbor empty sites during a
695	>	simulation. Therefore, we were able to obtain $G$'s for these
696	>	interfaces even at a relatively high temperature without being
697	>	affected by surface reconstructions.
698
699		\subsection{Influence of Capping Agent Coverage on $G$}
700		To investigate the influence of butanethiol coverage on interfacial
701		thermal conductance, a series of different coverage Au-butanethiol
702		surfaces is prepared and solvated with various organic
703		molecules. These systems are then equilibrated and their interfacial
704	<	thermal conductivity are measured with our NIVS algorithm. Table
705	<	\ref{tlnUhxnUhxnD} lists these results for direct comparison between
706	<	different coverages of butanethiol. To study the isotope effect in
707	<	interfacial thermal conductance, deuterated UA-hexane is included as
708	<	well.
704	>	thermal conductivity are measured with our NIVS algorithm. Figure
705	>	\ref{coverage} demonstrates the trend of conductance change with
706	>	respect to different coverages of butanethiol. To study the isotope
707	>	effect in interfacial thermal conductance, deuterated UA-hexane is
708	>	included as well.
709
710	+	\begin{figure}
711	+	\includegraphics[width=\linewidth]{coverage}
712	+	\caption{Comparison of interfacial thermal conductivity ($G$) values
713	+	for the Au-butanethiol/solvent interface with various UA models and
714	+	different capping agent coverages at $\langle T\rangle\sim$200K
715	+	using certain energy flux respectively.}
716	+	\label{coverage}
717	+	\end{figure}
718	+
719		It turned out that with partial covered butanethiol on the Au(111)
720	<	surface, the derivative definition for $G$ (Eq. \ref{derivativeG}) has
721	<	difficulty to apply, due to the difficulty in locating the maximum of
722	<	change of $\lambda$. Instead, the discrete definition
723	<	(Eq. \ref{discreteG}) is easier to apply, as max($\Delta T$) can still
724	<	be well-defined. Therefore, $G$'s (not $G^\prime$) are used for this
725	<	section.
720	>	surface, the derivative definition for $G^\prime$
721	>	(Eq. \ref{derivativeG}) was difficult to apply, due to the difficulty
722	>	in locating the maximum of change of $\lambda$. Instead, the discrete
723	>	definition (Eq. \ref{discreteG}) is easier to apply, as the Gibbs
724	>	deviding surface can still be well-defined. Therefore, $G$ (not
725	>	$G^\prime$) was used for this section.
726
727	<	From Table \ref{tlnUhxnUhxnD}, one can see the significance of the
727	>	From Figure \ref{coverage}, one can see the significance of the
728		presence of capping agents. Even when a fraction of the Au(111)
729		surface sites are covered with butanethiols, the conductivity would
730		see an enhancement by at least a factor of 3. This indicates the
731		important role cappping agent is playing for thermal transport
732	<	phenomena on metal/organic solvent surfaces.
732	>	phenomena on metal / organic solvent surfaces.
733
734		Interestingly, as one could observe from our results, the maximum
735		conductance enhancement (largest $G$) happens while the surfaces are
#	Line 607 \| Line 748 \| case, $G$ decrease could not be offset but instead acc
748		would not offset this effect. Eventually, when butanethiol coverage
749		continues to decrease, solvent-capping agent contact actually
750		decreases with the disappearing of butanethiol molecules. In this
751	<	case, $G$ decrease could not be offset but instead accelerated.
751	>	case, $G$ decrease could not be offset but instead accelerated. [NEED
752	>	SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
753
754		A comparison of the results obtained from differenet organic solvents
755		can also provide useful information of the interfacial thermal
#	Line 616 \| Line 758 \| differences in the IR range do not seem to produce an
758		butanethiol (UA) is non-deuterated for both solvents. These UA model
759		studies, even though eliminating C-H vibration samplings, still have
760		C-C vibrational frequencies different from each other. However, these
761	<	differences in the IR range do not seem to produce an observable
762	<	difference for the results of $G$. [MAY NEED FIGURE]
761	>	differences in the infrared range do not seem to produce an observable
762	>	difference for the results of $G$ (Figure \ref{uahxnua}).
763
764	+	\begin{figure}
765	+	\includegraphics[width=\linewidth]{uahxnua}
766	+	\caption{Vibrational spectra obtained for normal (upper) and
767	+	deuterated (lower) hexane in Au-butanethiol/hexane
768	+	systems. Butanethiol spectra are shown as reference. Both hexane and
769	+	butanethiol were using United-Atom models.}
770	+	\label{uahxnua}
771	+	\end{figure}
772	+
773		Furthermore, results for rigid body toluene solvent, as well as other
774		UA-hexane solvents, are reasonable within the general experimental
775	<	ranges[CITATIONS]. This suggests that explicit hydrogen might not be a
776	<	required factor for modeling thermal transport phenomena of systems
777	<	such as Au-thiol/organic solvent.
775	>	ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
776	>	suggests that explicit hydrogen might not be a required factor for
777	>	modeling thermal transport phenomena of systems such as
778	>	Au-thiol/organic solvent.
779
780		However, results for Au-butanethiol/toluene do not show an identical
781	<	trend with those for Au-butanethiol/hexane in that $G$'s remain at
781	>	trend with those for Au-butanethiol/hexane in that $G$ remains at
782		approximately the same magnitue when butanethiol coverage differs from
783		25\% to 75\%. This might be rooted in the molecule shape difference
784	<	for plane-like toluene and chain-like {\it n}-hexane. Due to this
784	>	for planar toluene and chain-like {\it n}-hexane. Due to this
785		difference, toluene molecules have more difficulty in occupying
786		relatively small gaps among capping agents when their coverage is not
787		too low. Therefore, the solvent-capping agent contact may keep
#	Line 638 \| Line 790 \| can see a plateau of $G$ vs. butanethiol coverage in o
790		its effect to the process of interfacial thermal transport. Thus, one
791		can see a plateau of $G$ vs. butanethiol coverage in our results.
792
641	–	[NEED ERROR ESTIMATE, MAY ALSO PUT J HERE]
642	–	\begin{table*}
643	–	\begin{minipage}{\linewidth}
644	–	\begin{center}
645	–	\caption{Computed interfacial thermal conductivity ($G$ in
646	–	MW/m$^2$/K) values for the Au-butanethiol/solvent interface
647	–	with various UA models and different capping agent coverages
648	–	at $<T>\sim$200K using certain energy flux respectively.}
649	–
650	–	\begin{tabular}{cccc}
651	–	\hline\hline
652	–	Thiol & & & \\
653	–	coverage (\%) & hexane & hexane-D & toluene \\
654	–	\hline
655	–	0.0 & 46.5 & 43.9 & 70.1 \\
656	–	25.0 & 151 & 153 & 249 \\
657	–	50.0 & 172 & 182 & 214 \\
658	–	75.0 & 242 & 229 & 244 \\
659	–	88.9 & 178 & - & - \\
660	–	100.0 & 137 & 153 & 187 \\
661	–	\hline\hline
662	–	\end{tabular}
663	–	\label{tlnUhxnUhxnD}
664	–	\end{center}
665	–	\end{minipage}
666	–	\end{table*}
667	–
793		\subsection{Influence of Chosen Molecule Model on $G$}
794	<	[MAY COMBINE W MECHANISM STUDY]
794	>	In addition to UA solvent/capping agent models, AA models are included
795	>	in our simulations as well. Besides simulations of the same (UA or AA)
796	>	model for solvent and capping agent, different models can be applied
797	>	to different components. Furthermore, regardless of models chosen,
798	>	either the solvent or the capping agent can be deuterated, similar to
799	>	the previous section. Table \ref{modelTest} summarizes the results of
800	>	these studies.
801
671	–	For the all-atom model, the liquid hexane phase was not stable under NPT
672	–	conditions. Therefore, the simulation length scale parameters are
673	–	adopted from previous equilibration results of the united-atom model
674	–	at 200K. Table \ref{AuThiolHexaneAA} shows the results of these
675	–	simulations. The conductivity values calculated with full capping
676	–	agent coverage are substantially larger than observed in the
677	–	united-atom model, and is even higher than predicted by
678	–	experiments. It is possible that our parameters for metal-non-metal
679	–	particle interactions lead to an overestimate of the interfacial
680	–	thermal conductivity, although the active C-H vibrations in the
681	–	all-atom model (which should not be appreciably populated at normal
682	–	temperatures) could also account for this high conductivity. The major
683	–	thermal transfer barrier of Au/butanethiol/hexane interface is between
684	–	the liquid phase and the capping agent, so extra degrees of freedom
685	–	such as the C-H vibrations could enhance heat exchange between these
686	–	two phases and result in a much higher conductivity.
687	–
802		\begin{table*}
803		\begin{minipage}{\linewidth}
804		\begin{center}
805
806		\caption{Computed interfacial thermal conductivity ($G$ and
807	<	$G^\prime$) values for the Au/butanethiol/hexane interface
808	<	with all-atom model and different capping agent coverage at
809	<	200K using a range of energy fluxes.}
807	>	$G^\prime$) values for interfaces using various models for
808	>	solvent and capping agent (or without capping agent) at
809	>	$\langle T\rangle\sim$200K. (D stands for deuterated solvent
810	>	or capping agent molecules; ``Avg.'' denotes results that are
811	>	averages of simulations under different $J_z$'s. Error
812	>	estimates indicated in parenthesis.)}
813
814	<	\begin{tabular}{cccc}
814	>	\begin{tabular}{llccc}
815		\hline\hline
816	<	Thiol & $J_z$ & $G$ & $G^\prime$ \\
817	<	coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
816	>	Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
817	>	(or bare surface) & model & (GW/m$^2$) &
818	>	\multicolumn{2}{c}{(MW/m$^2$/K)} \\
819		\hline
820	<	0.0 & 0.95 & 28.5 & 27.2 \\
821	<	& 1.88 & 30.3 & 28.9 \\
822	<	100.0 & 2.87 & 551 & 294 \\
823	<	& 3.81 & 494 & 193 \\
820	>	UA & UA hexane & Avg. & 131(9) & 87(10) \\
821	>	& UA hexane(D) & 1.95 & 153(5) & 136(13) \\
822	>	& AA hexane & Avg. & 131(6) & 122(10) \\
823	>	& UA toluene & 1.96 & 187(16) & 151(11) \\
824	>	& AA toluene & 1.89 & 200(36) & 149(53) \\
825	>	\hline
826	>	AA & UA hexane & 1.94 & 116(9) & 129(8) \\
827	>	& AA hexane & Avg. & 442(14) & 356(31) \\
828	>	& AA hexane(D) & 1.93 & 222(12) & 234(54) \\
829	>	& UA toluene & 1.98 & 125(25) & 97(60) \\
830	>	& AA toluene & 3.79 & 487(56) & 290(42) \\
831	>	\hline
832	>	AA(D) & UA hexane & 1.94 & 158(25) & 172(4) \\
833	>	& AA hexane & 1.92 & 243(29) & 191(11) \\
834	>	& AA toluene & 1.93 & 364(36) & 322(67) \\
835	>	\hline
836	>	bare & UA hexane & Avg. & 46.5(3.2) & 49.4(4.5) \\
837	>	& UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
838	>	& AA hexane & 0.96 & 31.0(1.4) & 29.4(1.3) \\
839	>	& UA toluene & 1.99 & 70.1(1.3) & 65.8(0.5) \\
840		\hline\hline
841		\end{tabular}
842	<	\label{AuThiolHexaneAA}
842	>	\label{modelTest}
843		\end{center}
844		\end{minipage}
845		\end{table*}
846
847	+	To facilitate direct comparison, the same system with differnt models
848	+	for different components uses the same length scale for their
849	+	simulation cells. Without the presence of capping agent, using
850	+	different models for hexane yields similar results for both $G$ and
851	+	$G^\prime$, and these two definitions agree with eath other very
852	+	well. This indicates very weak interaction between the metal and the
853	+	solvent, and is a typical case for acoustic impedance mismatch between
854	+	these two phases.
855
856	<	significant conductance enhancement compared to the gold/water
857	<	interface without capping agent and agree with available experimental
858	<	data. This indicates that the metal-metal potential, though not
859	<	predicting an accurate bulk metal thermal conductivity, does not
860	<	greatly interfere with the simulation of the thermal conductance
861	<	behavior across a non-metal interface.
856	>	As for Au(111) surfaces completely covered by butanethiols, the choice
857	>	of models for capping agent and solvent could impact the measurement
858	>	of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
859	>	interfaces, using AA model for both butanethiol and hexane yields
860	>	substantially higher conductivity values than using UA model for at
861	>	least one component of the solvent and capping agent, which exceeds
862	>	the general range of experimental measurement results. This is
863	>	probably due to the classically treated C-H vibrations in the AA
864	>	model, which should not be appreciably populated at normal
865	>	temperatures. In comparison, once either the hexanes or the
866	>	butanethiols are deuterated, one can see a significantly lower $G$ and
867	>	$G^\prime$. In either of these cases, the C-H(D) vibrational overlap
868	>	between the solvent and the capping agent is removed (Figure
869	>	\ref{aahxntln}). Conclusively, the improperly treated C-H vibration in
870	>	the AA model produced over-predicted results accordingly. Compared to
871	>	the AA model, the UA model yields more reasonable results with higher
872	>	computational efficiency.
873
874	<	% The results show that the two definitions used for $G$ yield
875	<	% comparable values, though $G^\prime$ tends to be smaller.
874	>	\begin{figure}
875	>	\includegraphics[width=\linewidth]{aahxntln}
876	>	\caption{Spectra obtained for All-Atom model Au-butanethil/solvent
877	>	systems. When butanethiol is deuterated (lower left), its
878	>	vibrational overlap with hexane would decrease significantly,
879	>	compared with normal butanethiol (upper left). However, this
880	>	dramatic change does not apply to toluene as much (right).}
881	>	\label{aahxntln}
882	>	\end{figure}
883
884	<	\subsection{Mechanism of Interfacial Thermal Conductance Enhancement
885	<	by Capping Agent}
886	<	[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL]
884	>	However, for Au-butanethiol/toluene interfaces, having the AA
885	>	butanethiol deuterated did not yield a significant change in the
886	>	measurement results. Compared to the C-H vibrational overlap between
887	>	hexane and butanethiol, both of which have alkyl chains, that overlap
888	>	between toluene and butanethiol is not so significant and thus does
889	>	not have as much contribution to the heat exchange
890	>	process. Conversely, extra degrees of freedom such as the C-H
891	>	vibrations could yield higher heat exchange rate between these two
892	>	phases and result in a much higher conductivity.
893
894	+	Although the QSC model for Au is known to predict an overly low value
895	+	for bulk metal gold conductivity\cite{kuang:164101}, our computational
896	+	results for $G$ and $G^\prime$ do not seem to be affected by this
897	+	drawback of the model for metal. Instead, our results suggest that the
898	+	modeling of interfacial thermal transport behavior relies mainly on
899	+	the accuracy of the interaction descriptions between components
900	+	occupying the interfaces.
901
902	<	%subsubsection{Vibrational spectrum study on conductance mechanism}
903	<	To investigate the mechanism of this interfacial thermal conductance,
904	<	the vibrational spectra of various gold systems were obtained and are
905	<	shown as in the upper panel of Fig. \ref{vibration}. To obtain these
906	<	spectra, one first runs a simulation in the NVE ensemble and collects
907	<	snapshots of configurations; these configurations are used to compute
908	<	the velocity auto-correlation functions, which is used to construct a
909	<	power spectrum via a Fourier transform. The gold surfaces covered by
910	<	butanethiol molecules exhibit an additional peak observed at a
911	<	frequency of $\sim$170cm$^{-1}$, which is attributed to the vibration
739	<	of the S-Au bond. This vibration enables efficient thermal transport
740	<	from surface Au atoms to the capping agents. Simultaneously, as shown
741	<	in the lower panel of Fig. \ref{vibration}, the large overlap of the
742	<	vibration spectra of butanethiol and hexane in the all-atom model,
743	<	including the C-H vibration, also suggests high thermal exchange
744	<	efficiency. The combination of these two effects produces the drastic
745	<	interfacial thermal conductance enhancement in the all-atom model.
902	>	\subsection{Role of Capping Agent in Interfacial Thermal Conductance}
903	>	The vibrational spectra for gold slabs in different environments are
904	>	shown as in Figure \ref{specAu}. Regardless of the presence of
905	>	solvent, the gold surfaces covered by butanethiol molecules, compared
906	>	to bare gold surfaces, exhibit an additional peak observed at the
907	>	frequency of $\sim$170cm$^{-1}$, which is attributed to the S-Au
908	>	bonding vibration. This vibration enables efficient thermal transport
909	>	from surface Au layer to the capping agents. Therefore, in our
910	>	simulations, the Au/S interfaces do not appear major heat barriers
911	>	compared to the butanethiol / solvent interfaces.
912
913	+	Simultaneously, the vibrational overlap between butanethiol and
914	+	organic solvents suggests higher thermal exchange efficiency between
915	+	these two components. Even exessively high heat transport was observed
916	+	when All-Atom models were used and C-H vibrations were treated
917	+	classically. Compared to metal and organic liquid phase, the heat
918	+	transfer efficiency between butanethiol and organic solvents is closer
919	+	to that within bulk liquid phase.
920	+
921	+	Furthermore, our observation validated previous
922	+	results\cite{hase:2010} that the intramolecular heat transport of
923	+	alkylthiols is highly effecient. As a combinational effects of these
924	+	phenomena, butanethiol acts as a channel to expedite thermal transport
925	+	process. The acoustic impedance mismatch between the metal and the
926	+	liquid phase can be effectively reduced with the presence of suitable
927	+	capping agents.
928	+
929		\begin{figure}
930		\includegraphics[width=\linewidth]{vibration}
931		\caption{Vibrational spectra obtained for gold in different
932	<	environments (upper panel) and for Au/thiol/hexane simulation in
933	<	all-atom model (lower panel).}
752	<	\label{vibration}
932	>	environments.}
933	>	\label{specAu}
934		\end{figure}
754	–	% MAY NEED TO CONVERT TO JPEG
935
936	+	[MAY ADD COMPARISON OF AU SLAB WIDTHS]
937	+
938		\section{Conclusions}
939	+	The NIVS algorithm we developed has been applied to simulations of
940	+	Au-butanethiol surfaces with organic solvents. This algorithm allows
941	+	effective unphysical thermal flux transferred between the metal and
942	+	the liquid phase. With the flux applied, we were able to measure the
943	+	corresponding thermal gradient and to obtain interfacial thermal
944	+	conductivities. Under steady states, single trajectory simulation
945	+	would be enough for accurate measurement. This would be advantageous
946	+	compared to transient state simulations, which need multiple
947	+	trajectories to produce reliable average results.
948
949	+	Our simulations have seen significant conductance enhancement with the
950	+	presence of capping agent, compared to the bare gold / liquid
951	+	interfaces. The acoustic impedance mismatch between the metal and the
952	+	liquid phase is effectively eliminated by proper capping
953	+	agent. Furthermore, the coverage precentage of the capping agent plays
954	+	an important role in the interfacial thermal transport
955	+	process. Moderately lower coverages allow higher contact between
956	+	capping agent and solvent, and thus could further enhance the heat
957	+	transfer process.
958
959	<	[NECESSITY TO STUDY THERMAL CONDUCTANCE IN NANOCRYSTAL STRUCTURE]\cite{vlugt:cpc2007154}
959	>	Our measurement results, particularly of the UA models, agree with
960	>	available experimental data. This indicates that our force field
961	>	parameters have a nice description of the interactions between the
962	>	particles at the interfaces. AA models tend to overestimate the
963	>	interfacial thermal conductance in that the classically treated C-H
964	>	vibration would be overly sampled. Compared to the AA models, the UA
965	>	models have higher computational efficiency with satisfactory
966	>	accuracy, and thus are preferable in interfacial thermal transport
967	>	modelings. 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 is
973	>	limited for accurate computation of derivatives data.
974
975	+	Vlugt {\it et al.} has investigated the surface thiol structures for
976	+	nanocrystal gold and pointed out that they differs from those of the
977	+	Au(111) surface\cite{vlugt:cpc2007154}. This difference might lead to
978	+	change of interfacial thermal transport behavior as well. To
979	+	investigate this problem, an effective means to introduce thermal flux
980	+	and measure the corresponding thermal gradient is desirable for
981	+	simulating structures with spherical symmetry.
982	+
983		\section{Acknowledgments}
984		Support for this project was provided by the National Science
985		Foundation under grant CHE-0848243. Computational time was provided by

Diff Legend

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

Comparing interfacial/interfacial.tex (file contents): Revision 3731 by skuang, Tue Jul 5 21:30:29 2011 UTC vs. Revision 3749 by skuang, Mon Jul 25 19:11:33 2011 UTC

Diff Legend

Comparing interfacial/interfacial.tex (file contents):
Revision 3731 by skuang, Tue Jul 5 21:30:29 2011 UTC vs.
Revision 3749 by skuang, Mon Jul 25 19:11:33 2011 UTC