ViewVC Help

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

root/group/interfacial/interfacial.tex

Comparing interfacial/interfacial.tex (file contents):
Revision 3730 by skuang, Tue Jul 5 17:39:21 2011 UTC vs.
Revision 3754 by gezelter, Thu Jul 28 20:25:08 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}
76	<	[BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM]
77	<	Interfacial thermal conductance is extensively studied both
78	<	experimentally and computationally, and systems with interfaces
79	<	present are generally heterogeneous. Although interfaces are commonly
80	<	barriers to heat transfer, it has been
81	<	reported\cite{doi:10.1021/la904855s} that under specific circustances,
82	<	e.g. with certain capping agents present on the surface, interfacial
83	<	conductance can be significantly enhanced. However, heat conductance
82	<	of molecular and nano-scale interfaces will be affected by the
83	<	chemical details of the surface and is challenging to
84	<	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.
76	>	Due to the importance of heat flow in nanotechnology, interfacial
77	>	thermal conductance has been studied extensively both experimentally
78	>	and computationally.\cite{cahill:793} Unlike bulk materials, nanoscale
79	>	materials have a significant fraction of their atoms at interfaces,
80	>	and the chemical details of these interfaces govern the heat transfer
81	>	behavior. Furthermore, the interfaces are
82	>	heterogeneous (e.g. solid - liquid), which provides a challenge to
83	>	traditional methods developed for homogeneous systems.
84
85	<	Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
85	>	Experimentally, various interfaces have been investigated for their
86	>	thermal conductance. Wang {\it et al.} studied heat transport through
87	>	long-chain hydrocarbon monolayers on gold substrate at individual
88	>	molecular level,\cite{Wang10082007} Schmidt {\it et al.} studied the
89	>	role of CTAB on thermal transport between gold nanorods and
90	>	solvent,\cite{doi:10.1021/jp8051888} and Juv\'e {\it et al.} studied
91	>	the cooling dynamics, which is controlled by thermal interface
92	>	resistence of glass-embedded metal
93	>	nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
94	>	normally considered barriers for heat transport, Alper {\it et al.}
95	>	suggested that specific ligands (capping agents) could completely
96	>	eliminate this barrier
97	>	($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}
98	>
99	>	Theoretical and computational models have also been used to study the
100	>	interfacial thermal transport in order to gain an understanding of
101	>	this phenomena at the molecular level. Recently, Hase and coworkers
102	>	employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
103	>	study thermal transport from hot Au(111) substrate to a self-assembled
104	>	monolayer of alkylthiol with relatively long chain (8-20 carbon
105	>	atoms).\cite{hase:2010,hase:2011} However, ensemble averaged
106	>	measurements for heat conductance of interfaces between the capping
107	>	monolayer on Au and a solvent phase have yet to be studied with their
108	>	approach. The comparatively low thermal flux through interfaces is
109	>	difficult to measure with Equilibrium MD or forward NEMD simulation
110	>	methods. Therefore, the Reverse NEMD (RNEMD)
111	>	methods\cite{MullerPlathe:1997xw,kuang:164101} would have the
112	>	advantage of applying this difficult to measure flux (while measuring
113	>	the resulting gradient), given that the simulation methods being able
114	>	to effectively apply an unphysical flux in non-homogeneous systems.
115	>	Garde and coworkers\cite{garde:nl2005,garde:PhysRevLett2009} applied
116	>	this approach to various liquid interfaces and studied how thermal
117	>	conductance (or resistance) is dependent on chemistry details of
118	>	interfaces, e.g. hydrophobic and hydrophilic aqueous interfaces.
119	>
120	>	Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS)
121		algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
122		retains the desirable features of RNEMD (conservation of linear
123		momentum and total energy, compatibility with periodic boundary
124		conditions) while establishing true thermal distributions in each of
125	<	the two slabs. Furthermore, it allows more effective thermal exchange
126	<	between particles of different identities, and thus enables extensive
127	<	study of interfacial conductance.
125	>	the two slabs. Furthermore, it allows effective thermal exchange
126	>	between particles of different identities, and thus makes the study of
127	>	interfacial conductance much simpler.
128
129	+	The work presented here deals with the Au(111) surface covered to
130	+	varying degrees by butanethiol, a capping agent with short carbon
131	+	chain, and solvated with organic solvents of different molecular
132	+	properties. Different models were used for both the capping agent and
133	+	the solvent force field parameters. Using the NIVS algorithm, the
134	+	thermal transport across these interfaces was studied and the
135	+	underlying mechanism for the phenomena was investigated.
136	+
137		\section{Methodology}
138	<	\subsection{Algorithm}
139	<	[BACKGROUND FOR MD METHODS]
140	<	There have been many algorithms for computing thermal conductivity
141	<	using molecular dynamics simulations. However, interfacial conductance
142	<	is at least an order of magnitude smaller. This would make the
143	<	calculation even more difficult for those slowly-converging
144	<	equilibrium methods. Imposed-flux non-equilibrium
145	<	methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
146	<	the response of temperature or momentum gradients are easier to
147	<	measure than the flux, if unknown, and thus, is a preferable way to
148	<	the forward NEMD methods. Although the momentum swapping approach for
149	<	flux-imposing can be used for exchanging energy between particles of
150	<	different identity, the kinetic energy transfer efficiency is affected
151	<	by the mass difference between the particles, which limits its
152	<	application on heterogeneous interfacial systems.
138	>	\subsection{Imposd-Flux Methods in MD Simulations}
139	>	Steady state MD simulations have an advantage in that not many
140	>	trajectories are needed to study the relationship between thermal flux
141	>	and thermal gradients. For systems with low interfacial conductance,
142	>	one must have a method capable of generating or measuring relatively
143	>	small fluxes, compared to those required for bulk conductivity. This
144	>	requirement makes the calculation even more difficult for
145	>	slowly-converging equilibrium methods.\cite{Viscardy:2007lq} Forward
146	>	NEMD methods impose a gradient (and measure a flux), but at interfaces
147	>	it is not clear what behavior should be imposed at the boundaries
148	>	between materials.  Imposed-flux reverse non-equilibrium
149	>	methods\cite{MullerPlathe:1997xw} set the flux {\it a priori} and
150	>	the thermal response becomes an easy-to-measure quantity.  Although
151	>	M\"{u}ller-Plathe's original momentum swapping approach can be used
152	>	for exchanging energy between particles of different identity, the
153	>	kinetic energy transfer efficiency is affected by the mass difference
154	>	between the particles, which limits its application on heterogeneous
155	>	interfacial systems.
156
157	<	The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach in
158	<	non-equilibrium MD simulations is able to impose relatively large
159	<	kinetic energy flux without obvious perturbation to the velocity
160	<	distribution of the simulated systems. Furthermore, this approach has
157	>	The non-isotropic velocity scaling (NIVS) \cite{kuang:164101} approach
158	>	to non-equilibrium MD simulations is able to impose a wide range of
159	>	kinetic energy fluxes without obvious perturbation to the velocity
160	>	distributions of the simulated systems. Furthermore, this approach has
161		the advantage in heterogeneous interfaces in that kinetic energy flux
162		can be applied between regions of particles of arbitary identity, and
163	<	the flux quantity is not restricted by particle mass difference.
163	>	the flux will not be restricted by difference in particle mass.
164
165		The NIVS algorithm scales the velocity vectors in two separate regions
166		of a simulation system with respective diagonal scaling matricies. To
167		determine these scaling factors in the matricies, a set of equations
168		including linear momentum conservation and kinetic energy conservation
169	<	constraints and target momentum/energy flux satisfaction is
170	<	solved. With the scaling operation applied to the system in a set
171	<	frequency, corresponding momentum/temperature gradients can be built,
172	<	which can be used for computing transportation properties and other
173	<	applications related to momentum/temperature gradients. The NIVS
132	<	algorithm conserves momenta and energy and does not depend on an
133	<	external thermostat.
169	>	constraints and target energy flux satisfaction is solved. With the
170	>	scaling operation applied to the system in a set frequency, bulk
171	>	temperature gradients can be easily established, and these can be used
172	>	for computing thermal conductivities. The NIVS algorithm conserves
173	>	momenta and energy and does not depend on an external thermostat.
174
175	<	\subsection{Defining Interfacial Thermal Conductivity $G$}
176	<	For interfaces with a relatively low interfacial conductance, the bulk
177	<	regions on either side of an interface rapidly come to a state in
178	<	which the two phases have relatively homogeneous (but distinct)
179	<	temperatures. The interfacial thermal conductivity $G$ can therefore
180	<	be approximated as:
175	>	\subsection{Defining Interfacial Thermal Conductivity ($G$)}
176	>
177	>	For an interface with relatively low interfacial conductance, and a
178	>	thermal flux between two distinct bulk regions, the regions on either
179	>	side of the interface rapidly come to a state in which the two phases
180	>	have relatively homogeneous (but distinct) temperatures. The
181	>	interfacial thermal conductivity $G$ can therefore be approximated as:
182		\begin{equation}
183	<	G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
183	>	G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
184		\langle T_\mathrm{cold}\rangle \right)}
185		\label{lowG}
186		\end{equation}
187	<	where ${E_{total}}$ is the imposed non-physical kinetic energy
188	<	transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle
189	<	T_\mathrm{cold}\rangle}$ are the average observed temperature of the
190	<	two separated phases.
187	>	where ${E_{total}}$ is the total imposed non-physical kinetic energy
188	>	transfer during the simulation and ${\langle T_\mathrm{hot}\rangle}$
189	>	and ${\langle T_\mathrm{cold}\rangle}$ are the average observed
190	>	temperature of the two separated phases.
191
192	<	When the interfacial conductance is {\it not} small, two ways can be
193	<	used to define $G$.
194	<
195	<	One way is to assume the temperature is discretely different on two
196	<	sides of the interface, $G$ can be calculated with the thermal flux
197	<	applied $J$ and the maximum temperature difference measured along the
157	<	thermal gradient max($\Delta T$), which occurs at the interface, as:
192	>	When the interfacial conductance is {\it not} small, there are two
193	>	ways to define $G$. One common way is to assume the temperature is
194	>	discrete on the two sides of the interface. $G$ can be calculated
195	>	using the applied thermal flux $J$ and the maximum temperature
196	>	difference measured along the thermal gradient max($\Delta T$), which
197	>	occurs at the Gibbs deviding surface (Figure \ref{demoPic}):
198		\begin{equation}
199	<	G=\frac{J}{\Delta T}
199	>	G=\frac{J}{\Delta T}
200		\label{discreteG}
201		\end{equation}
202
203	+	\begin{figure}
204	+	\includegraphics[width=\linewidth]{method}
205	+	\caption{Interfacial conductance can be calculated by applying an
206	+	(unphysical) kinetic energy flux between two slabs, one located
207	+	within the metal and another on the edge of the periodic box.  The
208	+	system responds by forming a thermal response or a gradient.  In
209	+	bulk liquids, this gradient typically has a single slope, but in
210	+	interfacial systems, there are distinct thermal conductivity
211	+	domains.  The interfacial conductance, $G$ is found by measuring the
212	+	temperature gap at the Gibbs dividing surface, or by using second
213	+	derivatives of the thermal profile.}
214	+	\label{demoPic}
215	+	\end{figure}
216	+
217		The other approach is to assume a continuous temperature profile along
218		the thermal gradient axis (e.g. $z$) and define $G$ at the point where
219	<	the magnitude of thermal conductivity $\lambda$ change reach its
219	>	the magnitude of thermal conductivity ($\lambda$) change reaches its
220		maximum, given that $\lambda$ is well-defined throughout the space:
221		\begin{equation}
222		G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
#	Line 173 | Line 227 | With the temperature profile obtained from simulations
227		\label{derivativeG}
228		\end{equation}
229
230	<	With the temperature profile obtained from simulations, one is able to
230	>	With temperature profiles obtained from simulation, one is able to
231		approximate the first and second derivatives of $T$ with finite
232	<	difference method and thus calculate $G^\prime$.
232	>	difference methods and calculate $G^\prime$. In what follows, both
233	>	definitions have been used, and are compared in the results.
234
235	<	In what follows, both definitions are used for calculation and comparison.
235	>	To investigate the interfacial conductivity at metal / solvent
236	>	interfaces, we have modeled a metal slab with its (111) surfaces
237	>	perpendicular to the $z$-axis of our simulation cells. The metal slab
238	>	has been prepared both with and without capping agents on the exposed
239	>	surface, and has been solvated with simple organic solvents, as
240	>	illustrated in Figure \ref{gradT}.
241
242	<	[IMPOSE G DEFINITION INTO OUR SYSTEMS]
243	<	To facilitate the use of the above definitions in calculating $G$ and
244	<	$G^\prime$, we have a metal slab with its (111) surfaces perpendicular
245	<	to the $z$-axis of our simulation cells. With or withour capping
246	<	agents on the surfaces, the metal slab is solvated with organic
247	<	solvents, as illustrated in Figure \ref{demoPic}.
242	>	With the simulation cell described above, we are able to equilibrate
243	>	the system and impose an unphysical thermal flux between the liquid
244	>	and the metal phase using the NIVS algorithm. By periodically applying
245	>	the unphysical flux, we obtained a temperature profile and its spatial
246	>	derivatives. Figure \ref{gradT} shows how an applied thermal flux can
247	>	be used to obtain the 1st and 2nd derivatives of the temperature
248	>	profile.
249
250		\begin{figure}
190	–	\includegraphics[width=\linewidth]{demoPic}
191	–	\caption{A sample showing how a metal slab has its (111) surface
192	–	covered by capping agent molecules and solvated by hexane.}
193	–	\label{demoPic}
194	–	\end{figure}
195	–
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	–
207	–	\begin{figure}
251		\includegraphics[width=\linewidth]{gradT}
252	<	\caption{The 1st and 2nd derivatives of temperature profile can be
253	<	obtained with finite difference approximation.}
252	>	\caption{A sample of Au-butanethiol/hexane interfacial system and the
253	>	temperature profile after a kinetic energy flux is imposed to
254	>	it. The 1st and 2nd derivatives of the temperature profile can be
255	>	obtained with finite difference approximation (lower panel).}
256		\label{gradT}
257		\end{figure}
258
214	–	[MAY INCLUDE POWER SPECTRUM PROTOCOL]
215	–
259		\section{Computational Details}
260		\subsection{Simulation Protocol}
261	<	In our simulations, Au is used to construct a metal slab with bare
262	<	(111) surface perpendicular to the $z$-axis. Different slab thickness
263	<	(layer numbers of Au) are simulated. This metal slab is first
264	<	equilibrated under normal pressure (1 atm) and a desired
265	<	temperature. After equilibration, butanethiol is used as the capping
266	<	agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
267	<	atoms in the butanethiol molecules would occupy the three-fold sites
268	<	of the surfaces, and the maximal butanethiol capacity on Au surface is
269	<	$1/3$ of the total number of surface Au atoms[CITATION]. A series of
270	<	different coverage surfaces is investigated in order to study the
271	<	relation between coverage and conductance.
261	>	The NIVS algorithm has been implemented in our MD simulation code,
262	>	OpenMD\cite{Meineke:2005gd,openmd}, and was used for our simulations.
263	>	Metal slabs of 6 or 11 layers of Au atoms were first equilibrated
264	>	under atmospheric pressure (1 atm) and 200K. After equilibration,
265	>	butanethiol capping agents were placed at three-fold hollow sites on
266	>	the Au(111) surfaces. These sites are either {\it fcc} or {\it
267	>	hcp} sites, although Hase {\it et al.} found that they are
268	>	equivalent in a heat transfer process,\cite{hase:2010} so we did not
269	>	distinguish between these sites in our study. The maximum butanethiol
270	>	capacity on Au surface is $1/3$ of the total number of surface Au
271	>	atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
272	>	structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
273	>	series of lower coverages was also prepared by eliminating
274	>	butanethiols from the higher coverage surface in a regular manner. The
275	>	lower coverages were prepared in order to study the relation between
276	>	coverage and interfacial conductance.
277
278	<	[COVERAGE DISCRIPTION] However, since the interactions between surface
279	<	Au and butanethiol is non-bonded, the capping agent molecules are
280	<	allowed to migrate to an empty neighbor three-fold site during a
281	<	simulation. Therefore, the initial configuration would not severely
282	<	affect the sampling of a variety of configurations of the same
283	<	coverage, and the final conductance measurement would be an average
284	<	effect of these configurations explored in the simulations. [MAY NEED FIGURES]
278	>	The capping agent molecules were allowed to migrate during the
279	>	simulations. They distributed themselves uniformly and sampled a
280	>	number of three-fold sites throughout out study. Therefore, the
281	>	initial configuration does not noticeably affect the sampling of a
282	>	variety of configurations of the same coverage, and the final
283	>	conductance measurement would be an average effect of these
284	>	configurations explored in the simulations.
285
286	<	After the modified Au-butanethiol surface systems are equilibrated
287	<	under canonical ensemble, Packmol\cite{packmol} is used to pack
288	<	organic solvent molecules in the previously vacuum part of the
289	<	simulation cells, which guarantees that short range repulsive
290	<	interactions do not disrupt the simulations. Two solvents are
291	<	investigated, one which has little vibrational overlap with the
292	<	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]
286	>	After the modified Au-butanethiol surface systems were equilibrated in
287	>	the canonical (NVT) ensemble, organic solvent molecules were packed in
288	>	the previously empty part of the simulation cells.\cite{packmol} Two
289	>	solvents were investigated, one which has little vibrational overlap
290	>	with the alkanethiol and which has a planar shape (toluene), and one
291	>	which has similar vibrational frequencies to the capping agent and
292	>	chain-like shape ({\it n}-hexane).
293
294	<	The spacing filled by solvent molecules, i.e. the gap between
295	<	periodically repeated Au-butanethiol surfaces should be carefully
296	<	chosen. A very long length scale for the thermal gradient axis ($z$)
251	<	may cause excessively hot or cold temperatures in the middle of the
294	>	The simulation cells were not particularly extensive along the
295	>	$z$-axis, as a very long length scale for the thermal gradient may
296	>	cause excessively hot or cold temperatures in the middle of the
297		solvent region and lead to undesired phenomena such as solvent boiling
298		or freezing when a thermal flux is applied. Conversely, too few
299		solvent molecules would change the normal behavior of the liquid
300		phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
301	<	these extreme cases did not happen to our simulations. And the
302	<	corresponding spacing is usually $35 \sim 60$\AA.
301	>	these extreme cases did not happen to our simulations. The spacing
302	>	between periodic images of the gold interfaces is $45 \sim 75$\AA.
303
304	<	The initial configurations generated by Packmol are further
305	<	equilibrated with the $x$ and $y$ dimensions fixed, only allowing
306	<	length scale change in $z$ dimension. This is to ensure that the
307	<	equilibration of liquid phase does not affect the metal crystal
308	<	structure in $x$ and $y$ dimensions. Further equilibration are run
309	<	under NVT and then NVE ensembles.
304	>	The initial configurations generated are further equilibrated with the
305	>	$x$ and $y$ dimensions fixed, only allowing the $z$-length scale to
306	>	change. This is to ensure that the equilibration of liquid phase does
307	>	not affect the metal's crystalline structure. Comparisons were made
308	>	with simulations that allowed changes of $L_x$ and $L_y$ during NPT
309	>	equilibration. No substantial changes in the box geometry were noticed
310	>	in these simulations. After ensuring the liquid phase reaches
311	>	equilibrium at atmospheric pressure (1 atm), further equilibration was
312	>	carried out under canonical (NVT) and microcanonical (NVE) ensembles.
313
314	<	After the systems reach equilibrium, NIVS is implemented to impose a
315	<	periodic unphysical thermal flux between the metal and the liquid
316	<	phase. Most of our simulations are under an average temperature of
317	<	$\sim$200K. Therefore, this flux usually comes from the metal to the
314	>	After the systems reach equilibrium, NIVS was used to impose an
315	>	unphysical thermal flux between the metal and the liquid phases. Most
316	>	of our simulations were done under an average temperature of
317	>	$\sim$200K. Therefore, thermal flux usually came from the metal to the
318		liquid so that the liquid has a higher temperature and would not
319	<	freeze due to excessively low temperature. This induced temperature
320	<	gradient is stablized and the simulation cell is devided evenly into
321	<	N slabs along the $z$-axis and the temperatures of each slab are
322	<	recorded. When the slab width $d$ of each slab is the same, the
323	<	derivatives of $T$ with respect to slab number $n$ can be directly
324	<	used for $G^\prime$ calculations:
325	<	\begin{equation}
326	<	G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
319	>	freeze due to lowered temperatures. After this induced temperature
320	>	gradient had stablized, the temperature profile of the simulation cell
321	>	was recorded. To do this, the simulation cell is devided evenly into
322	>	$N$ slabs along the $z$-axis. The average temperatures of each slab
323	>	are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
324	>	the same, the derivatives of $T$ with respect to slab number $n$ can
325	>	be directly used for $G^\prime$ calculations: \begin{equation}
326	>	G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
327		\Big/\left(\frac{\partial T}{\partial z}\right)^2
328		= |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
329		\Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
#	Line 284 | Line 332 | G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}
332		\label{derivativeG2}
333		\end{equation}
334
335	+	All of the above simulation procedures use a time step of 1 fs. Each
336	+	equilibration stage took a minimum of 100 ps, although in some cases,
337	+	longer equilibration stages were utilized.
338	+
339		\subsection{Force Field Parameters}
340	<	Our simulations include various components. Therefore, force field
341	<	parameter descriptions are needed for interactions both between the
342	<	same type of particles and between particles of different species.
340	>	Our simulations include a number of chemically distinct components.
341	>	Figure \ref{demoMol} demonstrates the sites defined for both
342	>	United-Atom and All-Atom models of the organic solvent and capping
343	>	agents in our simulations. Force field parameters are needed for
344	>	interactions both between the same type of particles and between
345	>	particles of different species.
346
347	+	\begin{figure}
348	+	\includegraphics[width=\linewidth]{structures}
349	+	\caption{Structures of the capping agent and solvents utilized in
350	+	these simulations. The chemically-distinct sites (a-e) are expanded
351	+	in terms of constituent atoms for both United Atom (UA) and All Atom
352	+	(AA) force fields.  Most parameters are from
353	+	Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes,TraPPE-UA.thiols} (UA) and \protect\cite{OPLSAA} (AA). Cross-interactions with the Au atoms are given in Table \ref{MnM}.}
354	+	\label{demoMol}
355	+	\end{figure}
356	+
357		The Au-Au interactions in metal lattice slab is described by the
358		quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
359		potentials include zero-point quantum corrections and are
360		reparametrized for accurate surface energies compared to the
361	<	Sutton-Chen potentials\cite{Chen90}.
361	>	Sutton-Chen potentials.\cite{Chen90}
362
363	<	Figure [REF] demonstrates how we name our pseudo-atoms of the
364	<	molecules in our simulations.
365	<	[FIGURE FOR MOLECULE NOMENCLATURE]
301	<
302	<	For both solvent molecules, straight chain {\it n}-hexane and aromatic
303	<	toluene, United-Atom (UA) and All-Atom (AA) models are used
304	<	respectively. The TraPPE-UA
363	>	For the two solvent molecules, {\it n}-hexane and toluene, two
364	>	different atomistic models were utilized. Both solvents were modeled
365	>	using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
366		parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
367	<	for our UA solvent molecules. In these models, pseudo-atoms are
368	<	located at the carbon centers for alkyl groups. By eliminating
369	<	explicit hydrogen atoms, these models are simple and computationally
370	<	efficient, while maintains good accuracy. However, the TraPPE-UA for
371	<	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]
367	>	for our UA solvent molecules. In these models, sites are located at
368	>	the carbon centers for alkyl groups. Bonding interactions, including
369	>	bond stretches and bends and torsions, were used for intra-molecular
370	>	sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones
371	>	potentials are used.
372
373	+	By eliminating explicit hydrogen atoms, the TraPPE-UA models are
374	+	simple and computationally efficient, while maintaining good accuracy.
375	+	However, the TraPPE-UA model for alkanes is known to predict a slighly
376	+	lower boiling point than experimental values. This is one of the
377	+	reasons we used a lower average temperature (200K) for our
378	+	simulations. If heat is transferred to the liquid phase during the
379	+	NIVS simulation, the liquid in the hot slab can actually be
380	+	substantially warmer than the mean temperature in the simulation. The
381	+	lower mean temperatures therefore prevent solvent boiling.
382	+
383	+	For UA-toluene, the non-bonded potentials between intermolecular sites
384	+	have a similar Lennard-Jones formulation. The toluene molecules were
385	+	treated as a single rigid body, so there was no need for
386	+	intramolecular interactions (including bonds, bends, or torsions) in
387	+	this solvent model.
388	+
389		Besides the TraPPE-UA models, AA models for both organic solvents are
390	<	included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
391	<	force field is used. [MORE DETAILS]
392	<	For toluene, the United Force Field developed by Rapp\'{e} {\it et
393	<	al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]
390	>	included in our studies as well. The OPLS-AA\cite{OPLSAA} force fields
391	>	were used. For hexane, additional explicit hydrogen sites were
392	>	included. Besides bonding and non-bonded site-site interactions,
393	>	partial charges and the electrostatic interactions were added to each
394	>	CT and HC site. For toluene, a flexible model for the toluene molecule
395	>	was utilized which included bond, bend, torsion, and inversion
396	>	potentials to enforce ring planarity.
397
398	<	The capping agent in our simulations, the butanethiol molecules can
399	<	either use UA or AA model. The TraPPE-UA force fields includes
398	>	The butanethiol capping agent in our simulations, were also modeled
399	>	with both UA and AA model. The TraPPE-UA force field includes
400		parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
401		UA butanethiol model in our simulations. The OPLS-AA also provides
402		parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
403	<	surfaces do not have the hydrogen atom bonded to sulfur. To adapt this
404	<	change and derive suitable parameters for butanethiol adsorbed on
405	<	Au(111) surfaces, we adopt the S parameters from [CITATION CF VLUGT]
406	<	and modify parameters for its neighbor C atom for charge balance in
407	<	the molecule. Note that the model choice (UA or AA) of capping agent
408	<	can be different from the solvent. Regardless of model choice, the
409	<	force field parameters for interactions between capping agent and
410	<	solvent can be derived using Lorentz-Berthelot Mixing Rule:
403	>	surfaces do not have the hydrogen atom bonded to sulfur. To derive
404	>	suitable parameters for butanethiol adsorbed on Au(111) surfaces, we
405	>	adopt the S parameters from Luedtke and Landman\cite{landman:1998} and
406	>	modify the parameters for the CTS atom to maintain charge neutrality
407	>	in the molecule.  Note that the model choice (UA or AA) for the capping
408	>	agent can be different from the solvent. Regardless of model choice,
409	>	the force field parameters for interactions between capping agent and
410	>	solvent can be derived using Lorentz-Berthelot Mixing Rule:
411	>	\begin{eqnarray}
412	>	\sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
413	>	\epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
414	>	\end{eqnarray}
415
416	+	To describe the interactions between metal (Au) and non-metal atoms,
417	+	we refer to an adsorption study of alkyl thiols on gold surfaces by
418	+	Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
419	+	Lennard-Jones form of potential parameters for the interaction between
420	+	Au and pseudo-atoms CH$_x$ and S based on a well-established and
421	+	widely-used effective potential of Hautman and Klein for the Au(111)
422	+	surface.\cite{hautman:4994} As our simulations require the gold slab
423	+	to be flexible to accommodate thermal excitation, the pair-wise form
424	+	of potentials they developed was used for our study.
425
426	<	To describe the interactions between metal Au and non-metal capping
427	<	agent and solvent particles, we refer to an adsorption study of alkyl
428	<	thiols on gold surfaces by Vlugt {\it et
429	<	al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones
430	<	form of potential parameters for the interaction between Au and
431	<	pseudo-atoms CH$_x$ and S based on a well-established and widely-used
432	<	effective potential of Hautman and Klein[CITATION] for the Au(111)
433	<	surface. As our simulations require the gold lattice slab to be
434	<	non-rigid so that it could accommodate kinetic energy for thermal
435	<	transport study purpose, the pair-wise form of potentials is
436	<	preferred.
426	>	The potentials developed from {\it ab initio} calculations by Leng
427	>	{\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
428	>	interactions between Au and aromatic C/H atoms in toluene. However,
429	>	the Lennard-Jones parameters between Au and other types of particles,
430	>	(e.g. AA alkanes) have not yet been established. For these
431	>	interactions, the Lorentz-Berthelot mixing rule can be used to derive
432	>	effective single-atom LJ parameters for the metal using the fit values
433	>	for toluene. These are then used to construct reasonable mixing
434	>	parameters for the interactions between the gold and other atoms.
435	>	Table \ref{MnM} summarizes the ``metal/non-metal'' parameters used in
436	>	our simulations.
437
353	–	Besides, the potentials developed from {\it ab initio} calculations by
354	–	Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
355	–	interactions between Au and aromatic C/H atoms in toluene.[MORE DETAILS]
356	–
357	–	However, the Lennard-Jones parameters between Au and other types of
358	–	particles in our simulations are not yet well-established. For these
359	–	interactions, we attempt to derive their parameters using the Mixing
360	–	Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters
361	–	for Au is first extracted from the Au-CH$_x$ parameters by applying
362	–	the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
363	–	parameters in our simulations.
364	–
438		\begin{table*}
439		\begin{minipage}{\linewidth}
440		\begin{center}
441	<	\caption{Lennard-Jones parameters for Au-non-Metal
442	<	interactions in our simulations.}
443	<
444	<	\begin{tabular}{ccc}
441	>	\caption{Non-bonded interaction parameters (including cross
442	>	interactions with Au atoms) for both force fields used in this
443	>	work.}
444	>	\begin{tabular}{lllllll}
445		\hline\hline
446	<	Non-metal & $\sigma$/\AA & $\epsilon$/kcal/mol \\
446	>	& Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
447	>	$\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
448	>	& & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
449		\hline
450	<	S    & 2.40   & 8.465   \\
451	<	CH3  & 3.54   & 0.2146  \\
452	<	CH2  & 3.54   & 0.1749  \\
453	<	CT3  & 3.365  & 0.1373  \\
454	<	CT2  & 3.365  & 0.1373  \\
455	<	CTT  & 3.365  & 0.1373  \\
456	<	HC   & 2.865  & 0.09256 \\
457	<	CHar & 3.4625 & 0.1680  \\
458	<	CRar & 3.555  & 0.1604  \\
459	<	CA   & 3.173  & 0.0640  \\
460	<	HA   & 2.746  & 0.0414  \\
450	>	United Atom (UA)
451	>	&CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
452	>	&CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
453	>	&CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
454	>	&CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
455	>	\hline
456	>	All Atom (AA)
457	>	&CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
458	>	&CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
459	>	&CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
460	>	&HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
461	>	&CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
462	>	&HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
463	>	\hline
464	>	Both UA and AA
465	>	& S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
466		\hline\hline
467		\end{tabular}
468		\label{MnM}
#	Line 391 | Line 471 | parameters in our simulations.
471		\end{table*}
472
473
474	<	\section{Results and Discussions}
475	<	[MAY HAVE A BRIEF SUMMARY]
476	<	\subsection{How Simulation Parameters Affects $G$}
477	<	[MAY NOT PUT AT FIRST]
474	>	\section{Results}
475	>	There are many factors contributing to the measured interfacial
476	>	conductance; some of these factors are physically motivated
477	>	(e.g. coverage of the surface by the capping agent coverage and
478	>	solvent identity), while some are governed by parameters of the
479	>	methodology (e.g. applied flux and the formulas used to obtain the
480	>	conductance). In this section we discuss the major physical and
481	>	calculational effects on the computed conductivity.
482	>
483	>	\subsection{Effects due to capping agent coverage}
484	>
485	>	A series of different initial conditions with a range of surface
486	>	coverages was prepared and solvated with various with both of the
487	>	solvent molecules. These systems were then equilibrated and their
488	>	interfacial thermal conductivity was measured with our NIVS
489	>	algorithm. Figure \ref{coverage} demonstrates the trend of conductance
490	>	with respect to surface coverage.
491	>
492	>	\begin{figure}
493	>	\includegraphics[width=\linewidth]{coverage}
494	>	\caption{Comparison of interfacial thermal conductivity ($G$) values
495	>	for the Au-butanethiol/solvent interface with various UA models and
496	>	different capping agent coverages at $\langle T\rangle\sim$200K.}
497	>	\label{coverage}
498	>	\end{figure}
499	>
500	>
501	>	In partially covered butanethiol on the Au(111) surface, the
502	>	derivative definition for $G^\prime$ (Eq. \ref{derivativeG}) becomes
503	>	difficult to apply, as the location of maximum change of $\lambda$
504	>	becomes washed out.  The discrete definition (Eq. \ref{discreteG}) is
505	>	easier to apply, as the Gibbs dividing surface is still
506	>	well-defined. Therefore, $G$ (not $G^\prime$) was used in this
507	>	section.
508	>
509	>	From Figure \ref{coverage}, one can see the significance of the
510	>	presence of capping agents. When even a small fraction of the Au(111)
511	>	surface sites are covered with butanethiols, the conductivity exhibits
512	>	an enhancement by at least a factor of 3. This indicates the important
513	>	role cappping agents are playing for thermal transport at metal /
514	>	organic solvent surfaces.
515	>
516	>	We note a non-monotonic behavior in the interfacial conductance as a
517	>	function of surface coverage. The maximum conductance (largest $G$)
518	>	happens when the surfaces are about 75\% covered with butanethiol
519	>	caps.  The reason for this behavior is not entirely clear.  One
520	>	explanation is that incomplete butanethiol coverage allows small gaps
521	>	between butanethiols to form. These gaps can be filled by transient
522	>	solvent molecules.  These solvent molecules couple very strongly with
523	>	the hot capping agent molecules near the surface, and can then carry
524	>	(diffusively) the excess thermal energy away from the surface.
525	>
526	>	There appears to be a competition between the conduction of the
527	>	thermal energy away from the surface by the capping agents (enhanced
528	>	by greater coverage) and the coupling of the capping agents with the
529	>	solvent (enhanced by physical contact at lower coverages).  This
530	>	competition would lead to the non-monotonic coverage behavior observed
531	>	here.
532	>
533	>	A comparison of the results obtained from the two different organic
534	>	solvents can also provide useful information of the interfacial
535	>	thermal transport process. The deuterated hexane (UA) results do not
536	>	appear to be substantially different from those of normal hexane (UA),
537	>	given that butanethiol (UA) is non-deuterated for both solvents. The
538	>	UA models, even though they have eliminated C-H vibrational overlap,
539	>	still have significant overlap in the infrared spectra.  Because
540	>	differences in the infrared range do not seem to produce an observable
541	>	difference for the results of $G$ (Figure \ref{uahxnua}).
542	>
543	>	\begin{figure}
544	>	\includegraphics[width=\linewidth]{uahxnua}
545	>	\caption{Vibrational spectra obtained for normal (upper) and
546	>	deuterated (lower) hexane in Au-butanethiol/hexane
547	>	systems. Butanethiol spectra are shown as reference. Both hexane and
548	>	butanethiol were using United-Atom models.}
549	>	\label{uahxnua}
550	>	\end{figure}
551	>
552	>	Furthermore, results for rigid body toluene solvent, as well as other
553	>	UA-hexane solvents, are reasonable within the general experimental
554	>	ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
555	>	suggests that explicit hydrogen might not be a required factor for
556	>	modeling thermal transport phenomena of systems such as
557	>	Au-thiol/organic solvent.
558	>
559	>	However, results for Au-butanethiol/toluene do not show an identical
560	>	trend with those for Au-butanethiol/hexane in that $G$ remains at
561	>	approximately the same magnitue when butanethiol coverage differs from
562	>	25\% to 75\%. This might be rooted in the molecule shape difference
563	>	for planar toluene and chain-like {\it n}-hexane. Due to this
564	>	difference, toluene molecules have more difficulty in occupying
565	>	relatively small gaps among capping agents when their coverage is not
566	>	too low. Therefore, the solvent-capping agent contact may keep
567	>	increasing until the capping agent coverage reaches a relatively low
568	>	level. This becomes an offset for decreasing butanethiol molecules on
569	>	its effect to the process of interfacial thermal transport. Thus, one
570	>	can see a plateau of $G$ vs. butanethiol coverage in our results.
571	>
572	>	\subsection{Effects due to Solvent \& Solvent Models}
573	>	In addition to UA solvent/capping agent models, AA models are included
574	>	in our simulations as well. Besides simulations of the same (UA or AA)
575	>	model for solvent and capping agent, different models can be applied
576	>	to different components. Furthermore, regardless of models chosen,
577	>	either the solvent or the capping agent can be deuterated, similar to
578	>	the previous section. Table \ref{modelTest} summarizes the results of
579	>	these studies.
580	>
581	>	\begin{table*}
582	>	\begin{minipage}{\linewidth}
583	>	\begin{center}
584	>
585	>	\caption{Computed interfacial thermal conductivity ($G$ and
586	>	$G^\prime$) values for interfaces using various models for
587	>	solvent and capping agent (or without capping agent) at
588	>	$\langle T\rangle\sim$200K. (D stands for deuterated solvent
589	>	or capping agent molecules; ``Avg.'' denotes results that are
590	>	averages of simulations under different $J_z$'s. Error
591	>	estimates indicated in parenthesis.)}
592	>
593	>	\begin{tabular}{llccc}
594	>	\hline\hline
595	>	Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
596	>	(or bare surface) & model & (GW/m$^2$) &
597	>	\multicolumn{2}{c}{(MW/m$^2$/K)} \\
598	>	\hline
599	>	UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
600	>	& UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
601	>	& AA hexane    & Avg. & 131(6)    & 122(10)   \\
602	>	& UA toluene   & 1.96 & 187(16)   & 151(11)   \\
603	>	& AA toluene   & 1.89 & 200(36)   & 149(53)   \\
604	>	\hline
605	>	AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
606	>	& AA hexane    & Avg. & 442(14)   & 356(31)   \\
607	>	& AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
608	>	& UA toluene   & 1.98 & 125(25)   & 97(60)    \\
609	>	& AA toluene   & 3.79 & 487(56)   & 290(42)   \\
610	>	\hline
611	>	AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
612	>	& AA hexane    & 1.92 & 243(29)   & 191(11)   \\
613	>	& AA toluene   & 1.93 & 364(36)   & 322(67)   \\
614	>	\hline
615	>	bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
616	>	& UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
617	>	& AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
618	>	& UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
619	>	\hline\hline
620	>	\end{tabular}
621	>	\label{modelTest}
622	>	\end{center}
623	>	\end{minipage}
624	>	\end{table*}
625	>
626	>	To facilitate direct comparison, the same system with differnt models
627	>	for different components uses the same length scale for their
628	>	simulation cells. Without the presence of capping agent, using
629	>	different models for hexane yields similar results for both $G$ and
630	>	$G^\prime$, and these two definitions agree with eath other very
631	>	well. This indicates very weak interaction between the metal and the
632	>	solvent, and is a typical case for acoustic impedance mismatch between
633	>	these two phases.
634	>
635	>	As for Au(111) surfaces completely covered by butanethiols, the choice
636	>	of models for capping agent and solvent could impact the measurement
637	>	of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
638	>	interfaces, using AA model for both butanethiol and hexane yields
639	>	substantially higher conductivity values than using UA model for at
640	>	least one component of the solvent and capping agent, which exceeds
641	>	the general range of experimental measurement results. This is
642	>	probably due to the classically treated C-H vibrations in the AA
643	>	model, which should not be appreciably populated at normal
644	>	temperatures. In comparison, once either the hexanes or the
645	>	butanethiols are deuterated, one can see a significantly lower $G$ and
646	>	$G^\prime$. In either of these cases, the C-H(D) vibrational overlap
647	>	between the solvent and the capping agent is removed (Figure
648	>	\ref{aahxntln}). Conclusively, the improperly treated C-H vibration in
649	>	the AA model produced over-predicted results accordingly. Compared to
650	>	the AA model, the UA model yields more reasonable results with higher
651	>	computational efficiency.
652	>
653	>	\begin{figure}
654	>	\includegraphics[width=\linewidth]{aahxntln}
655	>	\caption{Spectra obtained for All-Atom model Au-butanethil/solvent
656	>	systems. When butanethiol is deuterated (lower left), its
657	>	vibrational overlap with hexane would decrease significantly,
658	>	compared with normal butanethiol (upper left). However, this
659	>	dramatic change does not apply to toluene as much (right).}
660	>	\label{aahxntln}
661	>	\end{figure}
662	>
663	>	However, for Au-butanethiol/toluene interfaces, having the AA
664	>	butanethiol deuterated did not yield a significant change in the
665	>	measurement results. Compared to the C-H vibrational overlap between
666	>	hexane and butanethiol, both of which have alkyl chains, that overlap
667	>	between toluene and butanethiol is not so significant and thus does
668	>	not have as much contribution to the heat exchange
669	>	process. Conversely, extra degrees of freedom such as the C-H
670	>	vibrations could yield higher heat exchange rate between these two
671	>	phases and result in a much higher conductivity.
672	>
673	>	Although the QSC model for Au is known to predict an overly low value
674	>	for bulk metal gold conductivity\cite{kuang:164101}, our computational
675	>	results for $G$ and $G^\prime$ do not seem to be affected by this
676	>	drawback of the model for metal. Instead, our results suggest that the
677	>	modeling of interfacial thermal transport behavior relies mainly on
678	>	the accuracy of the interaction descriptions between components
679	>	occupying the interfaces.
680	>
681	>	\subsection{Effects due to methodology and simulation parameters}
682	>
683		We have varied our protocol or other parameters of the simulations in
684		order to investigate how these factors would affect the measurement of
685		$G$'s. It turned out that while some of these parameters would not
#	Line 403 | Line 688 | during equilibrating the liquid phase. Due to the stif
688		results.
689
690		In some of our simulations, we allowed $L_x$ and $L_y$ to change
691	<	during equilibrating the liquid phase. Due to the stiffness of the Au
692	<	slab, $L_x$ and $L_y$ would not change noticeably after
693	<	equilibration. Although $L_z$ could fluctuate $\sim$1\% after a system
694	<	is fully equilibrated in the NPT ensemble, this fluctuation, as well
695	<	as those comparably smaller to $L_x$ and $L_y$, would not be magnified
696	<	on the calculated $G$'s, as shown in Table \ref{AuThiolHexaneUA}. This
697	<	insensivity to $L_i$ fluctuations allows reliable measurement of $G$'s
698	<	without the necessity of extremely cautious equilibration process.
691	>	during equilibrating the liquid phase. Due to the stiffness of the
692	>	crystalline Au structure, $L_x$ and $L_y$ would not change noticeably
693	>	after equilibration. Although $L_z$ could fluctuate $\sim$1\% after a
694	>	system is fully equilibrated in the NPT ensemble, this fluctuation, as
695	>	well as those of $L_x$ and $L_y$ (which is significantly smaller),
696	>	would not be magnified on the calculated $G$'s, as shown in Table
697	>	\ref{AuThiolHexaneUA}. This insensivity to $L_i$ fluctuations allows
698	>	reliable measurement of $G$'s without the necessity of extremely
699	>	cautious equilibration process.
700
701		As stated in our computational details, the spacing filled with
702		solvent molecules can be chosen within a range. This allows some
#	Line 422 | Line 708 | Our NIVS algorithm allows change of unphysical thermal
708		smaller system size would be preferable, given that the liquid phase
709		structure is not affected.
710
711	+	\subsubsection{Effects of applied flux}
712		Our NIVS algorithm allows change of unphysical thermal flux both in
713		direction and in quantity. This feature extends our investigation of
714		interfacial thermal conductance. However, the magnitude of this
#	Line 437 | Line 724 | $J_z$ is listed in Table \ref{AuThiolHexaneUA} and [RE
724		the thermal flux across the interface. For our simulations, we denote
725		$J_z$ to be positive when the physical thermal flux is from the liquid
726		to metal, and negative vice versa. The $G$'s measured under different
727	<	$J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These
728	<	results do not suggest that $G$ is dependent on $J_z$ within this flux
729	<	range. The linear response of flux to thermal gradient simplifies our
730	<	investigations in that we can rely on $G$ measurement with only a
731	<	couple $J_z$'s and do not need to test a large series of fluxes.
727	>	$J_z$ is listed in Table \ref{AuThiolHexaneUA} and
728	>	\ref{AuThiolToluene}. These results do not suggest that $G$ is
729	>	dependent on $J_z$ within this flux range. The linear response of flux
730	>	to thermal gradient simplifies our investigations in that we can rely
731	>	on $G$ measurement with only a couple $J_z$'s and do not need to test
732	>	a large series of fluxes.
733
446	–	%ADD MORE TO TABLE
734		\begin{table*}
735		\begin{minipage}{\linewidth}
736		\begin{center}
737		\caption{Computed interfacial thermal conductivity ($G$ and
738	<	$G^\prime$) values for the Au/butanethiol/hexane interface
739	<	with united-atom model and different capping agent coverage
740	<	and solvent molecule numbers at different temperatures using a
741	<	range of energy fluxes.}
742	<
743	<	\begin{tabular}{cccccc}
738	>	$G^\prime$) values for the 100\% covered Au-butanethiol/hexane
739	>	interfaces with UA model and different hexane molecule numbers
740	>	at different temperatures using a range of energy
741	>	fluxes. Error estimates indicated in parenthesis.}
742	>
743	>	\begin{tabular}{ccccccc}
744		\hline\hline
745	<	Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\
746	<	coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) &
745	>	$\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
746	>	$J_z$ & $G$ & $G^\prime$ \\
747	>	(K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
748		\multicolumn{2}{c}{(MW/m$^2$/K)} \\
749		\hline
750	<	0.0   & 200 & 200 & 0.96 & 43.3 & 42.7 \\
751	<	&     &     & 1.91 & 45.7 & 42.9 \\
752	<	&     & 166 & 0.96 & 43.1 & 53.4 \\
753	<	88.9  & 200 & 166 & 1.94 & 172  & 108  \\
754	<	100.0 & 250 & 200 & 0.96 & 81.8 & 67.0 \\
755	<	&     & 166 & 0.98 & 79.0 & 62.9 \\
756	<	&     &     & 1.44 & 76.2 & 64.8 \\
757	<	& 200 & 200 & 1.92 & 129  & 87.3 \\
758	<	&     &     & 1.93 & 131  & 77.5 \\
759	<	&     & 166 & 0.97 & 115  & 69.3 \\
760	<	&     &     & 1.94 & 125  & 87.1 \\
750	>	200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
751	>	& 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
752	>	&     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
753	>	&     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
754	>	&     &     &       &  1.91 & 139(10)   & 101(10)   \\
755	>	&     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
756	>	& 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
757	>	&     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
758	>	&     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
759	>	&     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
760	>	\hline
761	>	250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
762	>	&     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
763	>	& 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
764	>	&     & No  & 0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
765	>	&     &     &       &  1.44 & 76.2(5.0) & 64.8(3.8) \\
766	>	&     &     &       & -0.95 & 56.4(2.5) & 54.4(1.1) \\
767	>	&     &     &       & -1.85 & 47.8(1.1) & 53.5(1.5) \\
768		\hline\hline
769		\end{tabular}
770		\label{AuThiolHexaneUA}
#	Line 477 | Line 772 | Furthermore, we also attempted to increase system aver
772		\end{minipage}
773		\end{table*}
774
775	+	\subsubsection{Effects due to average temperature}
776	+
777		Furthermore, we also attempted to increase system average temperatures
778		to above 200K. These simulations are first equilibrated in the NPT
779		ensemble under normal pressure. As stated above, the TraPPE-UA model
#	Line 485 | Line 782 | butanethiol as well.[MAY NEED FIGURE] And this reduced
782		temperature is higher than 250K. Additionally, the equilibrated liquid
783		hexane density under 250K becomes lower than experimental value. This
784		expanded liquid phase leads to lower contact between hexane and
785	<	butanethiol as well.[MAY NEED FIGURE] And this reduced contact would
785	>	butanethiol as well.[MAY NEED SLAB DENSITY FIGURE]
786	>	And this reduced contact would
787		probably be accountable for a lower interfacial thermal conductance,
788		as shown in Table \ref{AuThiolHexaneUA}.
789
#	Line 500 | Line 798 | in that higher degree of contact could yield increased
798		important role in the thermal transport process across the interface
799		in that higher degree of contact could yield increased conductance.
800
503	–	[ADD SIGNS AND ERROR ESTIMATE TO TABLE]
801		\begin{table*}
802		\begin{minipage}{\linewidth}
803		\begin{center}
804		\caption{Computed interfacial thermal conductivity ($G$ and
805	<	$G^\prime$) values for the Au/butanethiol/toluene interface at
806	<	different temperatures using a range of energy fluxes.}
805	>	$G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
806	>	interface at different temperatures using a range of energy
807	>	fluxes. Error estimates indicated in parenthesis.}
808
809	<	\begin{tabular}{cccc}
809	>	\begin{tabular}{ccccc}
810		\hline\hline
811	<	$\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
812	<	(K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
811	>	$\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
812	>	(K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
813		\hline
814	<	200 & 1.86 & 180 & 135 \\
815	<	& 2.15 & 204 & 113 \\
816	<	& 3.93 & 175 & 114 \\
817	<	300 & 1.91 & 143 & 125 \\
818	<	& 4.19 & 134 & 113 \\
814	>	200 & 0.933 &  2.15 & 204(12) & 113(12) \\
815	>	&       & -1.86 & 180(3)  & 135(21) \\
816	>	&       & -3.93 & 176(5)  & 113(12) \\
817	>	\hline
818	>	300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
819	>	&       & -4.19 & 135(9)  & 113(12) \\
820		\hline\hline
821		\end{tabular}
822		\label{AuThiolToluene}
#	Line 545 | Line 844 | undertaken at $<T>\sim$200K.
844		reconstructions could eliminate the original $x$ and $y$ dimensional
845		homogeneity, measurement of $G$ is more difficult to conduct under
846		higher temperatures. Therefore, most of our measurements are
847	<	undertaken at $<T>\sim$200K.
847	>	undertaken at $\langle T\rangle\sim$200K.
848
849		However, when the surface is not completely covered by butanethiols,
850		the simulated system is more resistent to the reconstruction
851	<	above. Our Au-butanethiol/toluene system did not see this phenomena
852	<	even at $<T>\sim$300K. The Au(111) surfaces have a 90\% coverage of
853	<	butanethiols and have empty three-fold sites. These empty sites could
854	<	help prevent surface reconstruction in that they provide other means
855	<	of capping agent relaxation. It is observed that butanethiols can
856	<	migrate to their neighbor empty sites during a simulation. Therefore,
857	<	we were able to obtain $G$'s for these interfaces even at a relatively
858	<	high temperature without being affected by surface reconstructions.
851	>	above. Our Au-butanethiol/toluene system had the Au(111) surfaces 90\%
852	>	covered by butanethiols, but did not see this above phenomena even at
853	>	$\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
854	>	capping agents could help prevent surface reconstruction in that they
855	>	provide other means of capping agent relaxation. It is observed that
856	>	butanethiols can migrate to their neighbor empty sites during a
857	>	simulation. Therefore, we were able to obtain $G$'s for these
858	>	interfaces even at a relatively high temperature without being
859	>	affected by surface reconstructions.
860
561	–	\subsection{Influence of Capping Agent Coverage on $G$}
562	–	To investigate the influence of butanethiol coverage on interfacial
563	–	thermal conductance, a series of different coverage Au-butanethiol
564	–	surfaces is prepared and solvated with various organic
565	–	molecules. These systems are then equilibrated and their interfacial
566	–	thermal conductivity are measured with our NIVS algorithm. Table
567	–	\ref{tlnUhxnUhxnD} lists these results for direct comparison between
568	–	different coverages of butanethiol.
861
862	<	With high coverage of butanethiol on the gold surface,
571	<	the interfacial thermal conductance is enhanced
572	<	significantly. Interestingly, a slightly lower butanethiol coverage
573	<	leads to a moderately higher conductivity. This is probably due to
574	<	more solvent/capping agent contact when butanethiol molecules are
575	<	not densely packed, which enhances the interactions between the two
576	<	phases and lowers the thermal transfer barrier of this interface.
577	<	[COMPARE TO AU/WATER IN PAPER]
862	>	\section{Discussion}
863
864	+	\subsection{Capping agent acts as a vibrational coupler between solid
865	+	and solvent phases}
866	+	To investigate the mechanism of interfacial thermal conductance, the
867	+	vibrational power spectrum was computed. Power spectra were taken for
868	+	individual components in different simulations. To obtain these
869	+	spectra, simulations were run after equilibration, in the NVE
870	+	ensemble, and without a thermal gradient. Snapshots of configurations
871	+	were collected at a frequency that is higher than that of the fastest
872	+	vibrations occuring in the simulations. With these configurations, the
873	+	velocity auto-correlation functions can be computed:
874	+	\begin{equation}
875	+	C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
876	+	\label{vCorr}
877	+	\end{equation}
878	+	The power spectrum is constructed via a Fourier transform of the
879	+	symmetrized velocity autocorrelation function,
880	+	\begin{equation}
881	+	\hat{f}(\omega) =
882	+	\int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
883	+	\label{fourier}
884	+	\end{equation}
885
886	<	significant conductance enhancement compared to the gold/water
887	<	interface without capping agent and agree with available experimental
888	<	data. This indicates that the metal-metal potential, though not
889	<	predicting an accurate bulk metal thermal conductivity, does not
890	<	greatly interfere with the simulation of the thermal conductance
891	<	behavior across a non-metal interface.
892	<	The results show that the two definitions used for $G$ yield
893	<	comparable values, though $G^\prime$ tends to be smaller.
886	>	From Figure \ref{coverage}, one can see the significance of the
887	>	presence of capping agents. Even when a fraction of the Au(111)
888	>	surface sites are covered with butanethiols, the conductivity would
889	>	see an enhancement by at least a factor of 3. This indicates the
890	>	important role cappping agent is playing for thermal transport
891	>	phenomena on metal / organic solvent surfaces.
892	>
893	>	Interestingly, as one could observe from our results, the maximum
894	>	conductance enhancement (largest $G$) happens while the surfaces are
895	>	about 75\% covered with butanethiols. This again indicates that
896	>	solvent-capping agent contact has an important role of the thermal
897	>	transport process. Slightly lower butanethiol coverage allows small
898	>	gaps between butanethiols to form. And these gaps could be filled with
899	>	solvent molecules, which acts like ``heat conductors'' on the
900	>	surface. The higher degree of interaction between these solvent
901	>	molecules and capping agents increases the enhancement effect and thus
902	>	produces a higher $G$ than densely packed butanethiol arrays. However,
903	>	once this maximum conductance enhancement is reached, $G$ decreases
904	>	when butanethiol coverage continues to decrease. Each capping agent
905	>	molecule reaches its maximum capacity for thermal
906	>	conductance. Therefore, even higher solvent-capping agent contact
907	>	would not offset this effect. Eventually, when butanethiol coverage
908	>	continues to decrease, solvent-capping agent contact actually
909	>	decreases with the disappearing of butanethiol molecules. In this
910	>	case, $G$ decrease could not be offset but instead accelerated. [MAY NEED
911	>	SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
912
913	+	A comparison of the results obtained from differenet organic solvents
914	+	can also provide useful information of the interfacial thermal
915	+	transport process. The deuterated hexane (UA) results do not appear to
916	+	be much different from those of normal hexane (UA), given that
917	+	butanethiol (UA) is non-deuterated for both solvents. These UA model
918	+	studies, even though eliminating C-H vibration samplings, still have
919	+	C-C vibrational frequencies different from each other. However, these
920	+	differences in the infrared range do not seem to produce an observable
921	+	difference for the results of $G$ (Figure \ref{uahxnua}).
922
923	<	\begin{table*}
924	<	\begin{minipage}{\linewidth}
925	<	\begin{center}
926	<	\caption{Computed interfacial thermal conductivity ($G$ and
927	<	$G^\prime$) values for the Au/butanethiol/hexane interface
928	<	with united-atom model and different capping agent coverage
929	<	and solvent molecule numbers at different temperatures using a
930	<	range of energy fluxes.}
598	<
599	<	\begin{tabular}{cccccc}
600	<	\hline\hline
601	<	Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\
602	<	coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) &
603	<	\multicolumn{2}{c}{(MW/m$^2$/K)} \\
604	<	\hline
605	<	0.0   & 200 & 200 & 0.96 & 43.3 & 42.7 \\
606	<	&     &     & 1.91 & 45.7 & 42.9 \\
607	<	&     & 166 & 0.96 & 43.1 & 53.4 \\
608	<	88.9  & 200 & 166 & 1.94 & 172  & 108  \\
609	<	100.0 & 250 & 200 & 0.96 & 81.8 & 67.0 \\
610	<	&     & 166 & 0.98 & 79.0 & 62.9 \\
611	<	&     &     & 1.44 & 76.2 & 64.8 \\
612	<	& 200 & 200 & 1.92 & 129  & 87.3 \\
613	<	&     &     & 1.93 & 131  & 77.5 \\
614	<	&     & 166 & 0.97 & 115  & 69.3 \\
615	<	&     &     & 1.94 & 125  & 87.1 \\
616	<	\hline\hline
617	<	\end{tabular}
618	<	\label{tlnUhxnUhxnD}
619	<	\end{center}
620	<	\end{minipage}
621	<	\end{table*}
923	>	\begin{figure}
924	>	\includegraphics[width=\linewidth]{uahxnua}
925	>	\caption{Vibrational spectra obtained for normal (upper) and
926	>	deuterated (lower) hexane in Au-butanethiol/hexane
927	>	systems. Butanethiol spectra are shown as reference. Both hexane and
928	>	butanethiol were using United-Atom models.}
929	>	\label{uahxnua}
930	>	\end{figure}
931
932	<	\subsection{Influence of Chosen Molecule Model on $G$}
933	<	[MAY COMBINE W MECHANISM STUDY]
932	>	Furthermore, results for rigid body toluene solvent, as well as other
933	>	UA-hexane solvents, are reasonable within the general experimental
934	>	ranges\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406}. This
935	>	suggests that explicit hydrogen might not be a required factor for
936	>	modeling thermal transport phenomena of systems such as
937	>	Au-thiol/organic solvent.
938
939	<	For the all-atom model, the liquid hexane phase was not stable under NPT
940	<	conditions. Therefore, the simulation length scale parameters are
941	<	adopted from previous equilibration results of the united-atom model
942	<	at 200K. Table \ref{AuThiolHexaneAA} shows the results of these
943	<	simulations. The conductivity values calculated with full capping
944	<	agent coverage are substantially larger than observed in the
945	<	united-atom model, and is even higher than predicted by
946	<	experiments. It is possible that our parameters for metal-non-metal
947	<	particle interactions lead to an overestimate of the interfacial
948	<	thermal conductivity, although the active C-H vibrations in the
949	<	all-atom model (which should not be appreciably populated at normal
950	<	temperatures) could also account for this high conductivity. The major
638	<	thermal transfer barrier of Au/butanethiol/hexane interface is between
639	<	the liquid phase and the capping agent, so extra degrees of freedom
640	<	such as the C-H vibrations could enhance heat exchange between these
641	<	two phases and result in a much higher conductivity.
939	>	However, results for Au-butanethiol/toluene do not show an identical
940	>	trend with those for Au-butanethiol/hexane in that $G$ remains at
941	>	approximately the same magnitue when butanethiol coverage differs from
942	>	25\% to 75\%. This might be rooted in the molecule shape difference
943	>	for planar toluene and chain-like {\it n}-hexane. Due to this
944	>	difference, toluene molecules have more difficulty in occupying
945	>	relatively small gaps among capping agents when their coverage is not
946	>	too low. Therefore, the solvent-capping agent contact may keep
947	>	increasing until the capping agent coverage reaches a relatively low
948	>	level. This becomes an offset for decreasing butanethiol molecules on
949	>	its effect to the process of interfacial thermal transport. Thus, one
950	>	can see a plateau of $G$ vs. butanethiol coverage in our results.
951
952	+	\subsection{Influence of Chosen Molecule Model on $G$}
953	+	In addition to UA solvent/capping agent models, AA models are included
954	+	in our simulations as well. Besides simulations of the same (UA or AA)
955	+	model for solvent and capping agent, different models can be applied
956	+	to different components. Furthermore, regardless of models chosen,
957	+	either the solvent or the capping agent can be deuterated, similar to
958	+	the previous section. Table \ref{modelTest} summarizes the results of
959	+	these studies.
960	+
961		\begin{table*}
962		\begin{minipage}{\linewidth}
963		\begin{center}
964
965		\caption{Computed interfacial thermal conductivity ($G$ and
966	<	$G^\prime$) values for the Au/butanethiol/hexane interface
967	<	with all-atom model and different capping agent coverage at
968	<	200K using a range of energy fluxes.}
966	>	$G^\prime$) values for interfaces using various models for
967	>	solvent and capping agent (or without capping agent) at
968	>	$\langle T\rangle\sim$200K. (D stands for deuterated solvent
969	>	or capping agent molecules; ``Avg.'' denotes results that are
970	>	averages of simulations under different $J_z$'s. Error
971	>	estimates indicated in parenthesis.)}
972
973	<	\begin{tabular}{cccc}
973	>	\begin{tabular}{llccc}
974		\hline\hline
975	<	Thiol & $J_z$ & $G$ & $G^\prime$ \\
976	<	coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
975	>	Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
976	>	(or bare surface) & model & (GW/m$^2$) &
977	>	\multicolumn{2}{c}{(MW/m$^2$/K)} \\
978		\hline
979	<	0.0   & 0.95 & 28.5 & 27.2 \\
980	<	& 1.88 & 30.3 & 28.9 \\
981	<	100.0 & 2.87 & 551  & 294  \\
982	<	& 3.81 & 494  & 193  \\
979	>	UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
980	>	& UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
981	>	& AA hexane    & Avg. & 131(6)    & 122(10)   \\
982	>	& UA toluene   & 1.96 & 187(16)   & 151(11)   \\
983	>	& AA toluene   & 1.89 & 200(36)   & 149(53)   \\
984	>	\hline
985	>	AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
986	>	& AA hexane    & Avg. & 442(14)   & 356(31)   \\
987	>	& AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
988	>	& UA toluene   & 1.98 & 125(25)   & 97(60)    \\
989	>	& AA toluene   & 3.79 & 487(56)   & 290(42)   \\
990	>	\hline
991	>	AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
992	>	& AA hexane    & 1.92 & 243(29)   & 191(11)   \\
993	>	& AA toluene   & 1.93 & 364(36)   & 322(67)   \\
994	>	\hline
995	>	bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
996	>	& UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
997	>	& AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
998	>	& UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
999		\hline\hline
1000		\end{tabular}
1001	<	\label{AuThiolHexaneAA}
1001	>	\label{modelTest}
1002		\end{center}
1003		\end{minipage}
1004		\end{table*}
1005
1006	+	To facilitate direct comparison, the same system with differnt models
1007	+	for different components uses the same length scale for their
1008	+	simulation cells. Without the presence of capping agent, using
1009	+	different models for hexane yields similar results for both $G$ and
1010	+	$G^\prime$, and these two definitions agree with eath other very
1011	+	well. This indicates very weak interaction between the metal and the
1012	+	solvent, and is a typical case for acoustic impedance mismatch between
1013	+	these two phases.
1014
1015	<	\subsection{Mechanism of Interfacial Thermal Conductance Enhancement
1016	<	by Capping Agent}
1017	<	[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL]
1015	>	As for Au(111) surfaces completely covered by butanethiols, the choice
1016	>	of models for capping agent and solvent could impact the measurement
1017	>	of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
1018	>	interfaces, using AA model for both butanethiol and hexane yields
1019	>	substantially higher conductivity values than using UA model for at
1020	>	least one component of the solvent and capping agent, which exceeds
1021	>	the general range of experimental measurement results. This is
1022	>	probably due to the classically treated C-H vibrations in the AA
1023	>	model, which should not be appreciably populated at normal
1024	>	temperatures. In comparison, once either the hexanes or the
1025	>	butanethiols are deuterated, one can see a significantly lower $G$ and
1026	>	$G^\prime$. In either of these cases, the C-H(D) vibrational overlap
1027	>	between the solvent and the capping agent is removed (Figure
1028	>	\ref{aahxntln}). Conclusively, the improperly treated C-H vibration in
1029	>	the AA model produced over-predicted results accordingly. Compared to
1030	>	the AA model, the UA model yields more reasonable results with higher
1031	>	computational efficiency.
1032
1033	+	\begin{figure}
1034	+	\includegraphics[width=\linewidth]{aahxntln}
1035	+	\caption{Spectra obtained for All-Atom model Au-butanethil/solvent
1036	+	systems. When butanethiol is deuterated (lower left), its
1037	+	vibrational overlap with hexane would decrease significantly,
1038	+	compared with normal butanethiol (upper left). However, this
1039	+	dramatic change does not apply to toluene as much (right).}
1040	+	\label{aahxntln}
1041	+	\end{figure}
1042
1043	<	%subsubsection{Vibrational spectrum study on conductance mechanism}
1044	<	To investigate the mechanism of this interfacial thermal conductance,
1045	<	the vibrational spectra of various gold systems were obtained and are
1046	<	shown as in the upper panel of Fig. \ref{vibration}. To obtain these
1047	<	spectra, one first runs a simulation in the NVE ensemble and collects
1048	<	snapshots of configurations; these configurations are used to compute
1049	<	the velocity auto-correlation functions, which is used to construct a
1050	<	power spectrum via a Fourier transform. The gold surfaces covered by
1051	<	butanethiol molecules exhibit an additional peak observed at a
683	<	frequency of $\sim$170cm$^{-1}$, which is attributed to the vibration
684	<	of the S-Au bond. This vibration enables efficient thermal transport
685	<	from surface Au atoms to the capping agents. Simultaneously, as shown
686	<	in the lower panel of Fig. \ref{vibration}, the large overlap of the
687	<	vibration spectra of butanethiol and hexane in the all-atom model,
688	<	including the C-H vibration, also suggests high thermal exchange
689	<	efficiency. The combination of these two effects produces the drastic
690	<	interfacial thermal conductance enhancement in the all-atom model.
1043	>	However, for Au-butanethiol/toluene interfaces, having the AA
1044	>	butanethiol deuterated did not yield a significant change in the
1045	>	measurement results. Compared to the C-H vibrational overlap between
1046	>	hexane and butanethiol, both of which have alkyl chains, that overlap
1047	>	between toluene and butanethiol is not so significant and thus does
1048	>	not have as much contribution to the heat exchange
1049	>	process. Conversely, extra degrees of freedom such as the C-H
1050	>	vibrations could yield higher heat exchange rate between these two
1051	>	phases and result in a much higher conductivity.
1052
1053	+	Although the QSC model for Au is known to predict an overly low value
1054	+	for bulk metal gold conductivity\cite{kuang:164101}, our computational
1055	+	results for $G$ and $G^\prime$ do not seem to be affected by this
1056	+	drawback of the model for metal. Instead, our results suggest that the
1057	+	modeling of interfacial thermal transport behavior relies mainly on
1058	+	the accuracy of the interaction descriptions between components
1059	+	occupying the interfaces.
1060	+
1061	+	\subsection{Role of Capping Agent in Interfacial Thermal Conductance}
1062	+	The vibrational spectra for gold slabs in different environments are
1063	+	shown as in Figure \ref{specAu}. Regardless of the presence of
1064	+	solvent, the gold surfaces covered by butanethiol molecules, compared
1065	+	to bare gold surfaces, exhibit an additional peak observed at the
1066	+	frequency of $\sim$170cm$^{-1}$, which is attributed to the S-Au
1067	+	bonding vibration. This vibration enables efficient thermal transport
1068	+	from surface Au layer to the capping agents. Therefore, in our
1069	+	simulations, the Au/S interfaces do not appear major heat barriers
1070	+	compared to the butanethiol / solvent interfaces.
1071	+
1072	+	\subsubsection{Overlap of power spectrum}
1073	+	Simultaneously, the vibrational overlap between butanethiol and
1074	+	organic solvents suggests higher thermal exchange efficiency between
1075	+	these two components. Even exessively high heat transport was observed
1076	+	when All-Atom models were used and C-H vibrations were treated
1077	+	classically. Compared to metal and organic liquid phase, the heat
1078	+	transfer efficiency between butanethiol and organic solvents is closer
1079	+	to that within bulk liquid phase.
1080	+
1081	+	Furthermore, our observation validated previous
1082	+	results\cite{hase:2010} that the intramolecular heat transport of
1083	+	alkylthiols is highly effecient. As a combinational effects of these
1084	+	phenomena, butanethiol acts as a channel to expedite thermal transport
1085	+	process. The acoustic impedance mismatch between the metal and the
1086	+	liquid phase can be effectively reduced with the presence of suitable
1087	+	capping agents.
1088	+
1089		\begin{figure}
1090		\includegraphics[width=\linewidth]{vibration}
1091		\caption{Vibrational spectra obtained for gold in different
1092	<	environments (upper panel) and for Au/thiol/hexane simulation in
1093	<	all-atom model (lower panel).}
697	<	\label{vibration}
1092	>	environments.}
1093	>	\label{specAu}
1094		\end{figure}
699	–	% MAY NEED TO CONVERT TO JPEG
1095
1096	+	[MAY ADD COMPARISON OF AU SLAB WIDTHS BUT NOT MUCH TO TALK ABOUT...]
1097	+
1098		\section{Conclusions}
1099	+	The NIVS algorithm we developed has been applied to simulations of
1100	+	Au-butanethiol surfaces with organic solvents. This algorithm allows
1101	+	effective unphysical thermal flux transferred between the metal and
1102	+	the liquid phase. With the flux applied, we were able to measure the
1103	+	corresponding thermal gradient and to obtain interfacial thermal
1104	+	conductivities. Under steady states, single trajectory simulation
1105	+	would be enough for accurate measurement. This would be advantageous
1106	+	compared to transient state simulations, which need multiple
1107	+	trajectories to produce reliable average results.
1108
1109	+	Our simulations have seen significant conductance enhancement with the
1110	+	presence of capping agent, compared to the bare gold / liquid
1111	+	interfaces. The acoustic impedance mismatch between the metal and the
1112	+	liquid phase is effectively eliminated by proper capping
1113	+	agent. Furthermore, the coverage precentage of the capping agent plays
1114	+	an important role in the interfacial thermal transport
1115	+	process. Moderately lower coverages allow higher contact between
1116	+	capping agent and solvent, and thus could further enhance the heat
1117	+	transfer process.
1118
1119	<	[NECESSITY TO STUDY THERMAL CONDUCTANCE IN NANOCRYSTAL STRUCTURE]\cite{vlugt:cpc2007154}
1119	>	Our measurement results, particularly of the UA models, agree with
1120	>	available experimental data. This indicates that our force field
1121	>	parameters have a nice description of the interactions between the
1122	>	particles at the interfaces. AA models tend to overestimate the
1123	>	interfacial thermal conductance in that the classically treated C-H
1124	>	vibration would be overly sampled. Compared to the AA models, the UA
1125	>	models have higher computational efficiency with satisfactory
1126	>	accuracy, and thus are preferable in interfacial thermal transport
1127	>	modelings. Of the two definitions for $G$, the discrete form
1128	>	(Eq. \ref{discreteG}) was easier to use and gives out relatively
1129	>	consistent results, while the derivative form (Eq. \ref{derivativeG})
1130	>	is not as versatile. Although $G^\prime$ gives out comparable results
1131	>	and follows similar trend with $G$ when measuring close to fully
1132	>	covered or bare surfaces, the spatial resolution of $T$ profile is
1133	>	limited for accurate computation of derivatives data.
1134
1135	+	Vlugt {\it et al.} has investigated the surface thiol structures for
1136	+	nanocrystal gold and pointed out that they differs from those of the
1137	+	Au(111) surface\cite{landman:1998,vlugt:cpc2007154}. This difference
1138	+	might lead to change of interfacial thermal transport behavior as
1139	+	well. To investigate this problem, an effective means to introduce
1140	+	thermal flux and measure the corresponding thermal gradient is
1141	+	desirable for simulating structures with spherical symmetry.
1142	+
1143		\section{Acknowledgments}
1144		Support for this project was provided by the National Science
1145		Foundation under grant CHE-0848243. Computational time was provided by
1146		the Center for Research Computing (CRC) at the University of Notre
1147	<	Dame. \newpage
1147	>	Dame.
1148	>	\newpage
1149
1150		\bibliography{interfacial}
1151

Diff Legend

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