ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/chainLength/GoldThiolsPaper.tex
Revision: 3850
Committed: Fri Dec 21 19:29:36 2012 UTC (12 years, 7 months ago) by gezelter
Content type: application/x-tex
File size: 35671 byte(s)
Log Message:
edits

File Contents

# Content
1 \documentclass[11pt]{article}
2 \usepackage{amsmath}
3 \usepackage{amssymb}
4 \usepackage{times}
5 \usepackage{mathptm}
6 \usepackage{setspace}
7 \usepackage{endfloat}
8 \usepackage{caption}
9 \usepackage{tabularx}
10 \usepackage{longtable}
11 \usepackage{graphicx}
12 \usepackage{multirow}
13 \usepackage{multicol}
14 \usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions
15 \usepackage[square, comma, sort&compress]{natbib}
16 \usepackage{url}
17 \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
18 \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
19 9.0in \textwidth 6.5in \brokenpenalty=10000
20
21 % double space list of tables and figures
22 %\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
23 \setlength{\abovecaptionskip}{20 pt}
24 \setlength{\belowcaptionskip}{30 pt}
25
26 \bibpunct{}{}{,}{s}{}{;}
27
28 \citestyle{nature}
29 \bibliographystyle{achemso}
30
31 \begin{document}
32
33 \newcolumntype{A}{p{1.5in}}
34 \newcolumntype{B}{p{0.75in}}
35
36 \title{Simulations of heat conduction at thiolate-capped gold
37 surfaces: The role of chain length and solvent penetration}
38
39 \author{Kelsey M. Stocker and J. Daniel
40 Gezelter\footnote{Corresponding author. \ Electronic mail:
41 gezelter@nd.edu} \\
42 251 Nieuwland Science Hall, \\
43 Department of Chemistry and Biochemistry,\\
44 University of Notre Dame\\
45 Notre Dame, Indiana 46556}
46
47 \date{\today}
48
49 \maketitle
50
51 \begin{doublespace}
52
53 \begin{abstract}
54
55 \end{abstract}
56
57 \newpage
58
59 %\narrowtext
60
61 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62 % **INTRODUCTION**
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 \section{Introduction}
65
66 The structural details of the interfaces of metal nanoparticles
67 determines how energy flows between these particles and their
68 surroundings. Understanding this energy flow is essential in designing
69 and functionalizing metallic nanoparticles for plasmonic photothermal
70 therapies,\cite{Jain:2007ux,Petrova:2007ad,Gnyawali:2008lp,Mazzaglia:2008to,Huff:2007ye,Larson:2007hw}
71 which rely on the ability of metallic nanoparticles to absorb light in
72 the near-IR, a portion of the spectrum in which living tissue is very
73 nearly transparent. The relevant physical property controlling the
74 transfer of this energy as heat into the surrounding tissue is the
75 interfacial thermal conductance, $G$, which can be somewhat difficult
76 to determine experimentally.\cite{Wilson:2002uq,Plech:2005kx}
77
78 Metallic particles have also been proposed for use in highly efficient
79 thermal-transfer fluids, although careful experiments by Eapen {\it et
80 al.} have shown that metal-particle-based nanofluids have thermal
81 conductivities that match Maxwell predictions.\cite{Eapen:2007th} The
82 likely cause of previously reported non-Maxwell
83 behavior\cite{Eastman:2001wb,Keblinski:2002bx,Lee:1999ct,Xue:2003ya,Xue:2004oa}
84 is percolation networks of nanoparticles exchanging energy via the
85 solvent,\cite{Eapen:2007mw} so it is important to get a detailed
86 molecular picture of particle-ligand and ligand-solvent interactions
87 in order to understand the thermal behavior of complex fluids. To
88 date, there have been few reported values (either from theory or
89 experiment) for $G$ for ligand-protected nanoparticles embedded in
90 liquids, and there is a significant gap in knowledge about how
91 chemically distinct ligands or protecting groups will affect heat
92 transport from the particles.
93
94 Experimentally, the thermal properties of various nanostructured
95 interfaces have been investigated by a number of groups. Cahill and
96 coworkers studied thermal transport from metal nanoparticle/fluid
97 interfaces, epitaxial TiN/single crystal oxide interfaces, and
98 hydrophilic and hydrophobic interfaces between water and solids with
99 different self-assembled
100 monolayers.\cite{cahill:793,Wilson:2002uq,PhysRevB.67.054302,doi:10.1021/jp048375k,PhysRevLett.96.186101}
101 Wang {\it et al.} studied heat transport through long-chain
102 hydrocarbon monolayers on gold substrate at the individual molecular
103 level,\cite{Wang10082007} Schmidt {\it et al.} studied the role of
104 cetyltrimethylammonium bromide (CTAB) on the thermal transport between
105 gold nanorods and solvent,\cite{doi:10.1021/jp8051888} and Juv\'e {\it
106 et al.} studied the cooling dynamics, which is controlled by thermal
107 interface resistance of glass-embedded metal
108 nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
109 normally considered barriers for heat transport, Alper {\it et al.}
110 have suggested that specific ligands (capping agents) could completely
111 eliminate this barrier
112 ($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}
113
114 In previous simulations, we applied a variant of reverse
115 non-equilibrium molecular dynamics (RNEMD) to calculate the
116 interfacial thermal conductance at a metal / organic solvent interface
117 that had been chemically protected by butanethiolate groups. Our
118 calculations suggested an explanation for the very large thermal
119 conductivity at alkanethiol-capped metal surfaces when compared with
120 bare metal/solvent interfaces. Specifically, the chemical bond
121 between the metal and the ligand introduces a vibrational overlap that
122 is not present without the protecting group, and the overlap between
123 the vibrational spectra (metal to ligand, ligand to solvent) provides
124 a mechanism for rapid thermal transport across the interface.
125
126 A notable result of the previous simulations was the non-monotonic
127 dependence of $G$ on the fractional coverage of the metal surface by
128 the chemical protecting group. Gaps in surface coverage allow the
129 solvent molecules come into direct contact with ligands that had been
130 heated by the metal surface, absorb thermal energy from the ligands,
131 and then diffuse away. Quantifying the role of surface coverage is
132 difficult as the ligands have lateral mobility on the surface and can
133 aggregate to form domains on the timescale of the simulation.
134
135 To isolate this effect without worrying about lateral mobility of the
136 surface ligands, the current work involves mixed-chain-length
137 monolayers in which the length mismatch between long and short chains
138 is sufficient to accomodate solvent molecules. These completely
139 covered (but mixed-chain) surfaces are significantly less prone to
140 surface domain formation, and allow us to further investigate the
141 mechanism of heat transport to the solvent. A thermal flux is created
142 using velocity shearing and scaling reverse non-equilibrium molecular
143 dynamics (VSS-RNEMD), and the resulting temperature profiles are
144 analyzed to yield information about the interfacial thermal
145 conductance.
146
147
148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 % **METHODOLOGY**
150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151 \section{Methodology}
152
153 There are many ways to compute bulk transport properties from
154 classical molecular dynamics simulations. Equilibrium Molecular
155 Dynamics (EMD) simulations can be used by computing relevant time
156 correlation functions and assuming that linear response theory
157 holds.\cite{PhysRevB.37.5677,MASSOBRIO:1984bl,PhysRev.119.1,Viscardy:2007rp,che:6888,kinaci:014106}
158 For some transport properties (notably thermal conductivity), EMD
159 approaches are subject to noise and poor convergence of the relevant
160 correlation functions. Traditional Non-equilibrium Molecular Dynamics
161 (NEMD) methods impose a gradient (e.g. thermal or momentum) on a
162 simulation.\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2002dp,PhysRevA.34.1449,JiangHao_jp802942v}
163 However, the resulting flux is often difficult to
164 measure. Furthermore, problems arise for NEMD simulations of
165 heterogeneous systems, such as phase-phase boundaries or interfaces,
166 where the type of gradient to enforce at the boundary between
167 materials is unclear.
168
169 {\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt
170 a different approach in that an unphysical {\it flux} is imposed
171 between different regions or ``slabs'' of the simulation
172 box.\cite{MullerPlathe:1997xw,ISI:000080382700030,Kuang2010} The
173 system responds by developing a temperature or momentum {\it gradient}
174 between the two regions. Since the amount of the applied flux is known
175 exactly, and the measurement of a gradient is generally less
176 complicated, imposed-flux methods typically take shorter simulation
177 times to obtain converged results for transport properties. The
178 corresponding temperature or velocity gradients which develop in
179 response to the applied flux are then related (via linear response
180 theory) to the transport coefficient of interest. These methods are
181 quite efficient, in that they do not need many trajectories to provide
182 information about transport properties. To date, they have been
183 utilized in computing thermal and mechanical transfer of both
184 homogeneous
185 liquids~\cite{MullerPlathe:1997xw,ISI:000080382700030,Maginn:2010} as
186 well as heterogeneous
187 systems.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl}
188
189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190 % VSS-RNEMD
191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
192 \subsection{VSS-RNEMD}
193 The original ``swapping'' approaches by M\"{u}ller-Plathe {\it et
194 al.}\cite{ISI:000080382700030,MullerPlathe:1997xw} can be understood
195 as a sequence of imaginary elastic collisions between particles in
196 regions separated by half of the simulation cell. In each collision,
197 the entire momentum vectors of both particles may be exchanged to
198 generate a thermal flux. Alternatively, a single component of the
199 momentum vectors may be exchanged to generate a shear flux. This
200 algorithm turns out to be quite useful in many simulations of bulk
201 liquids. However, the M\"{u}ller-Plathe swapping approach perturbs the
202 system away from ideal Maxwell-Boltzmann distributions, and this has
203 undesirable side-effects when the applied flux becomes
204 large.\cite{Maginn:2010}
205
206 The most useful alternative RNEMD approach developed so far utilizes a
207 series of simultaneous velocity shearing and scaling exchanges between
208 the two slabs.\cite{Kuang2012} This method provides a set of
209 conservation constraints while simultaneously creating a desired flux
210 between the two slabs. Satisfying the constraint equations ensures
211 that the new configurations are sampled from the same NVE ensemble.
212
213 The VSS moves are applied periodically to scale and shift the particle
214 velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and
215 $C$) which are separated by half of the simulation box,
216 \begin{displaymath}
217 \begin{array}{rclcl}
218
219 & \underline{\mathrm{shearing}} & &
220 \underline{~~~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~~~} \\ \\
221 \mathbf{v}_i \leftarrow &
222 \mathbf{a}_c\ & + & c\cdot\left(\mathbf{v}_i - \langle\mathbf{v}_c
223 \rangle\right) + \langle\mathbf{v}_c\rangle \\
224 \mathbf{v}_j \leftarrow &
225 \mathbf{a}_h & + & h\cdot\left(\mathbf{v}_j - \langle\mathbf{v}_h
226 \rangle\right) + \langle\mathbf{v}_h\rangle .
227
228 \end{array}
229 \end{displaymath}
230 Here $\langle\mathbf{v}_c\rangle$ and $\langle\mathbf{v}_h\rangle$ are
231 the center of mass velocities in the $C$ and $H$ slabs, respectively.
232 Within the two slabs, particles receive incremental changes or a
233 ``shear'' to their velocities. The amount of shear is governed by the
234 imposed momentum flux, $j_z(\mathbf{p})$
235 \begin{eqnarray}
236 \mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\
237 \mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2}
238 \end{eqnarray}
239 where $M_{\{c,h\}}$ is total mass of particles within each slab and $\Delta t$
240 is the interval between two separate operations.
241
242 To simultaneously impose a thermal flux ($J_z$) between the slabs we
243 use energy conservation constraints,
244 \begin{eqnarray}
245 K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\mathbf{v}_c
246 \rangle^2) + \frac{1}{2}M_c (\langle \mathbf{v}_c \rangle + \mathbf{a}_c)^2 \label{vss3}\\
247 K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\mathbf{v}_h
248 \rangle^2) + \frac{1}{2}M_h (\langle \mathbf{v}_h \rangle +
249 \mathbf{a}_h)^2 \label{vss4}.
250 \label{constraint}
251 \end{eqnarray}
252 Simultaneous solution of these quadratic formulae for the scaling
253 coefficients, $c$ and $h$, will ensure that the simulation samples from
254 the original microcanonical (NVE) ensemble. Here $K_{\{c,h\}}$ is the
255 instantaneous translational kinetic energy of each slab. At each time
256 interval, it is a simple matter to solve for $c$, $h$, $\mathbf{a}_c$,
257 and $\mathbf{a}_h$, subject to the imposed momentum flux,
258 $j_z(\mathbf{p})$, and thermal flux, $J_z$ values. Since the VSS
259 operations do not change the kinetic energy due to orientational
260 degrees of freedom or the potential energy of a system, configurations
261 after the VSS move have exactly the same energy ({\it and} linear
262 momentum) as before the move.
263
264 As the simulation progresses, the VSS moves are performed on a regular
265 basis, and the system develops a thermal or velocity gradient in
266 response to the applied flux. Using the slope of the temperature or
267 velocity gradient, it is quite simple to obtain of thermal
268 conductivity ($\lambda$),
269 \begin{equation}
270 J_z = -\lambda \frac{\partial T}{\partial z},
271 \end{equation}
272 and shear viscosity ($\eta$),
273 \begin{equation}
274 j_z(p_x) = -\eta \frac{\partial \langle v_x \rangle}{\partial z}.
275 \end{equation}
276 Here, the quantities on the left hand side are the imposed flux
277 values, while the slopes are obtained from linear fits to the
278 gradients that develop in the simulated system.
279
280 The VSS-RNEMD approach is versatile in that it may be used to
281 implement both thermal and shear transport either separately or
282 simultaneously. Perturbations of velocities away from the ideal
283 Maxwell-Boltzmann distributions are minimal, and thermal anisotropy is
284 kept to a minimum. This ability to generate simultaneous thermal and
285 shear fluxes has been previously utilized to map out the shear
286 viscosity of SPC/E water over a wide range of temperatures (90~K) with
287 a single 1 ns simulation.\cite{Kuang2012}
288
289 \begin{figure}
290 \includegraphics[width=\linewidth]{figures/rnemd}
291 \caption{The VSS-RNEMD approach imposes unphysical transfer of
292 linear momentum or kinetic energy between a ``hot'' slab and a
293 ``cold'' slab in the simulation box. The system responds to this
294 imposed flux by generating velocity or temperature gradients. The
295 slope of the gradients can then be used to compute transport
296 properties (e.g. shear viscosity or thermal conductivity).}
297 \label{fig:rnemd}
298 \end{figure}
299
300 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301 % INTERFACIAL CONDUCTANCE
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303 \subsection{Reverse Non-Equilibrium Molecular Dynamics approaches
304 to interfacial transport}
305
306 Interfaces between dissimilar materials have transport properties
307 which can be defined as derivatives of the standard transport
308 coefficients in a direction normal to the interface. For example, the
309 {\it interfacial} thermal conductance ($G$) can be thought of as the
310 change in the thermal conductivity ($\lambda$) across the boundary
311 between materials:
312 \begin{align}
313 G &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\
314 &= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/
315 \left(\frac{\partial T}{\partial z}\right)_{z_0}^2.
316 \label{derivativeG}
317 \end{align}
318 where $z_0$ is the location of the interface between two materials and
319 $\mathbf{\hat{n}}$ is a unit vector normal to the interface (assumed
320 to be the $z$ direction from here on). RNEMD simulations, and
321 particularly the VSS-RNEMD approach, function by applying a momentum
322 or thermal flux and watching the gradient response of the
323 material. This means that the {\it interfacial} conductance is a
324 second derivative property which is subject to significant noise and
325 therefore requires extensive sampling. We have been able to
326 demonstrate the use of the second derivative approach to compute
327 interfacial conductance at chemically-modified metal / solvent
328 interfaces. However, a definition of the interfacial conductance in
329 terms of a temperature difference ($\Delta T$) measured at $z_0$,
330 \begin{displaymath}
331 G = \frac{J_z}{\Delta T_{z_0}},
332 \end{displaymath}
333 is useful once the RNEMD approach has generated a stable temperature
334 gap across the interface.
335
336 \begin{figure}
337 \includegraphics[width=\linewidth]{figures/resistor_series}
338 \caption{The inverse of the interfacial thermal conductance, $G$, is
339 the Kapitza resistance, $R_K$. Because the gold / thiolate/
340 solvent interface extends a significant distance from the metal
341 surface, the interfacial resistance $R_K$ can be computed by
342 summing a series of temperature drops between adjacent temperature
343 bins along the $z$ axis.}
344 \label{fig:resistor_series}
345 \end{figure}
346
347 In the particular case we are studying here, there are two interfaces
348 involved in the transfer of heat from the gold slab to the solvent:
349 the metal/thiolate interface and the thiolate/solvent interface. We
350 could treat the temperature on each side of an interface as discrete,
351 making the interfacial conductance the inverse of the Kaptiza
352 resistance, or $G = \frac{1}{R_k}$. To model the total conductance
353 across multiple interfaces, it is useful to think of the interfaces as
354 a set of resistors wired in series. The total resistance is then
355 additive, $R_{total} = \sum_i R_{i}$ and the interfacial conductance
356 is the inverse of the total resistance, or $G = \frac{1}{\sum_i
357 R_i}$). In the interfacial region, we treat each bin in the
358 VSS-RNEMD temperature profile as a resistor with resistance
359 $\frac{T_{2}-T_{1}}{J_z}$, $\frac{T_{3}-T_{2}}{J_z}$, etc. The sum of
360 the set of resistors which spans the gold/thiolate interface, thiolate
361 chains, and thiolate/solvent interface simplifies to
362 \begin{equation}
363 R_{K} = \frac{T_{n}-T_{1}}{J_z},
364 \label{eq:finalG}
365 \end{equation}
366 or the temperature difference between the gold side of the
367 gold/thiolate interface and the solvent side of the thiolate/solvent
368 interface over the applied flux.
369
370
371 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
372 % **COMPUTATIONAL DETAILS**
373 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
374 \section{Computational Details}
375
376 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
377 % FORCE-FIELD PARAMETERS
378 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
379 \subsection{Force-Field Parameters}
380
381 Our simulations include a number of chemically distinct components.
382 Figure \ref{fig:structures} demonstrates the sites defined for both
383 the {\it n}-hexane and alkanethiolate ligands present in our
384 simulations.
385
386 \begin{figure}
387 \includegraphics[width=\linewidth]{figures/structures}
388 \caption{Topologies of the thiolate capping agents and solvent
389 utilized in the simulations. The chemically-distinct sites (S,
390 \ce{CH2}, and \ce{CH3}) are treated as united atoms. Most
391 parameters are taken from references \bibpunct{}{}{,}{n}{}{,}
392 \protect\cite{TraPPE-UA.alkanes} and
393 \protect\cite{TraPPE-UA.thiols}. Cross-interactions with the Au
394 atoms were adapted from references
395 \protect\cite{landman:1998},~\protect\cite{vlugt:cpc2007154},~and
396 \protect\cite{hautman:4994}.}
397 \label{fig:structures}
398 \end{figure}
399
400 The Au-Au interactions in the metal lattice slab are described by the
401 quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
402 potentials include zero-point quantum corrections and are
403 reparametrized for accurate surface energies compared to the
404 Sutton-Chen potentials.\cite{Chen90}
405
406 For the {\it n}-hexane solvent molecules, the TraPPE-UA
407 parameters\cite{TraPPE-UA.alkanes} were utilized. In this model,
408 sites are located at the carbon centers for alkyl groups. Bonding
409 interactions, including bond stretches and bends and torsions, were
410 used for intra-molecular sites closer than 3 bonds. For non-bonded
411 interactions, Lennard-Jones potentials were used. We have previously
412 utilized both united atom (UA) and all-atom (AA) force fields for
413 thermal conductivity,\cite{kuang:AuThl,Kuang2012} and since the united
414 atom force fields cannot populate the high-frequency modes that are
415 present in AA force fields, they appear to work better for modeling
416 thermal conductivity. The TraPPE-UA model for alkanes is known to
417 predict a slightly lower boiling point than experimental values. This
418 is one of the reasons we used a lower average temperature (200K) for
419 our simulations.
420
421 The TraPPE-UA force field includes parameters for thiol
422 molecules\cite{TraPPE-UA.thiols} which were used for the
423 alkanethiolate molecules in our simulations. To derive suitable
424 parameters for butanethiol adsorbed on Au(111) surfaces, we adopted
425 the S parameters from Luedtke and Landman\cite{landman:1998} and
426 modified the parameters for the CTS atom to maintain charge neutrality
427 in the molecule.
428
429 To describe the interactions between metal (Au) and non-metal atoms,
430 we refer to an adsorption study of alkyl thiols on gold surfaces by
431 Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
432 Lennard-Jones form of potential parameters for the interaction between
433 Au and pseudo-atoms CH$_x$ and S based on a well-established and
434 widely-used effective potential of Hautman and Klein for the Au(111)
435 surface.\cite{hautman:4994} As our simulations require the gold slab
436 to be flexible to accommodate thermal excitation, the pair-wise form
437 of potentials they developed was used for our study.
438
439 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
440 % SIMULATION PROTOCOL
441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
442 \subsection{Simulation Protocol}
443
444 We have implemented the VSS-RNEMD algorithm in OpenMD, our group
445 molecular dynamics code.\cite{openmd} An 1188 atom gold slab was
446 equilibrated at 1 atm and 200 K. The periodic box was then expanded
447 in the $z$-direction to expose two Au(111) faces on either side of the
448 11-atom thick slab.
449
450 A full monolayer of thiolates (1/3 the number of surface gold atoms)
451 was placed on three-fold hollow sites on the gold interfaces. The
452 effect of thiolate binding sites on the thermal conductance was tested
453 by placing thiolates at both fcc and hcp hollow sites. No appreciable
454 difference in the temperature profile was noted due to the location of
455 thiolate binding.
456
457 To test the role of thiolate chain length on interfacial thermal
458 conductance, full coverages of each of five chain lengths were tested:
459 butanethiolate (C$_4$), hexanethiolate (C$_6$), octanethiolate
460 (C$_8$), decanethiolate (C$_{10}$), and dodecanethiolate
461 (C$_{12}$). To test the effect of mixed chain lengths, full coverages
462 of C$_4$ / C$_{10}$ and C$_4$ / C$_{12}$ mixtures were created in
463 short/long chain percentages of 25/75, 50/50, 62.5/37.5, 75/25, and
464 87.5/12.5. The short and long chains were placed on the surface hollow
465 sites in a random configuration.
466
467 The simulation box $z$-dimension was set to roughly double the length
468 of the gold/thiolate block. Hexane solvent molecules were placed in
469 the vacant portion of the box using the packmol algorithm.
470
471 The system was equilibrated to 220 K in the canonical (NVT) ensemble,
472 allowing the thiolates and solvent to warm gradually. Pressure
473 correction to 1 atm was done using an isobaric-isothermal (NPT)
474 integrator that allowed expansion or contraction only in the $z$
475 direction, maintaining the crystalline structure of the gold as close
476 to the bulk result as possible. The diagonal elements of the pressure
477 tensor were monitored during the pressure equilibration stage. If the
478 $xx$ and/or $yy$ elements had a mean above zero throughout the
479 simulation -- indicating residual surface tension in the plane of the
480 gold slab -- an additional short NPT equilibration step was performed
481 allowing all box dimensions to change. Once the pressure was stable
482 at 1 atm, a final equilibration stage was performed at constant
483 temperature. All systems were equilibrated in the microcanonical (NVE)
484 ensemble before proceeding with the VSS-RNEMD and data collection
485 stages.
486
487 A kinetic energy flux was applied using VSS-RNEMD. The total
488 simulation time was 3 nanoseconds, with velocity scaling occurring
489 every 10 femtoseconds. The hot slab was centered in the gold and the
490 cold slab was placed in the center of the solvent region. The entire
491 system has a (time-averaged) temperature of 220 K, with a temperature
492 difference between the hot and cold slabs of approximately 30 K. The
493 average temperature and kinetic energy flux were selected to prevent
494 solvent freezing (or glass formation) and to not allow the thiolates
495 to bury in the gold slab. The Au-S interaction has a deep potential
496 energy well, which allows sulfur atoms to embed into the gold slab at
497 temperatures above 250 K. Simulation conditions were chosen to avoid
498 both of these undesirable situations.
499
500 Temperature profiles of the system were created by dividing the box
501 into $\sim$ 3 \AA \, bins along the $z$ axis and recording the average
502 temperature of each bin.
503
504
505
506 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
507 % **RESULTS**
508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
509 \section{Results}
510
511 The solvent, hexane, a straight chain flexible alkane, is structurally
512 similar to the thiolate alkane tails, and previous work has shown that
513 UA models of hexane and butanethiolate have a high degree of
514 vibrational overlap.\cite{kuang:AuThl} This overlap provides a
515 mechanism for thermal energy transfer from the thiolates to the
516 solvent.
517
518 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
519 % CHAIN LENGTH
520 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
521 \subsection{Effect of Chain Length}
522
523 We examined full coverages of five chain lengths, n = 4, 6, 8, 10, 12. As shown in table \ref{table:chainlengthG}, the interfacial conductance is constant as a function of chain length, indicating that the length of the thiolate alkyl chains does not play a significant role in the transport of heat across the gold/thiolate and thiolate/solvent interfaces. In all cases, the hexane solvent was unable to penetrate into the thiolate layer, leading to a persistent 2-4 \AA \, gap between the solvent region and the thiolates. However, while the identity of the alkyl thiolate capping agent has little effect on the interfacial thermal conductance, the presence of a full monolayer of capping agent provides a two-fold increase in the G value relative to a bare gold surface.
524 \begin{longtable}{p{4cm} p{3cm}}
525 \caption{Computed interfacial thermal conductance (G) values for bare gold and 100\% coverages of various thiolate alkyl chain lengths.}
526 \\
527 \centering {\bf Chain Length (n)} & \centering\arraybackslash {\bf G (MW/m$^2$/K)} \\ \hline
528 \endhead
529 \hline
530 \endfoot
531 \centering bare gold & \centering\arraybackslash 30.2 \\
532 \centering 4 & \centering\arraybackslash 59.4 \\
533 \centering 6 & \centering\arraybackslash 60.2 \\
534 \centering 8 & \centering\arraybackslash 61.0 \\
535 \centering 10 & \centering\arraybackslash 58.2 \\
536 \centering 12 & \centering\arraybackslash 58.8
537 \label{table:chainlengthG}
538 \end{longtable}
539
540 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
541 % MIXED CHAINS
542 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
543 \subsection{Effect of Mixed Chain Lengths}
544
545 Previous work demonstrated that for butanethiolate monolayers on a Au(111) surface, the interfacial conductance was a non-monotonic function of the percent coverage. This is believed to be due to enhanced solvent-thiolate coupling through greater penetration of solvent molecules into the thiolate layer. At lower coverages, hexane solvent can more penetrate into the interfacial layer by fitting into gaps between the thiolates. However, a side effect of low coverages is surface aggregation of the thiolates. To simulate the effect of low coverages while preventing thiolate domain formation, we maintain 100\% thiolate coverage while varying the proportions of short (butanethiolate, n = 4) and long (decanethiolate, n = 10 or dodecanethiolate, n = 12) thiolate alkyl chains. We show in figure \ref{fig:Gstacks} that the interfacial conductance of full coverage mixed-chain systems is a non-monotonic function of the fraction of long chains.
546 \begin{figure}
547 \includegraphics[width=\linewidth]{figures/Gstacks}
548 \caption{Interfacial thermal conductivity of mixed-chains has a non-monotonic dependence on the fraction of long chains (lower panels). At low long chain coverage, the solvent escape rate dominates the heat transfer process, while the solvent-thiolate orientational order dominates at high long chain coverage (upper panels).}
549 \label{fig:Gstacks}
550 \end{figure}
551 Mixtures of butanethiolate/decanethiolate (n = 4, 10) have a peak interfacial conductance for 25\%/75\% short/long chains. At this short/long chain fraction, the orientational ordering of the solvent molecules and chains becomes a dominant effect while the solvent escape rate sharply decreases. Mixtures of butanethiolate/dodecanethiolate (n = 4, 12) have a peak interfacial conductance for 12.5\%/87.5\% short/long chains. The solvent-thiolate orientational ordering reaches its maximum value at this short/long chain ratio. Long chain fractions of over $0.5$ for the n = 4, 12 system are very well ordered, but this effect is tempered by the exceptionally slow solvent escape rate.
552
553 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
554 % **DISCUSSION**
555 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
556 \section{Discussion}
557
558 In the mixed chain-length simulations, solvent molecules can become
559 temporarily trapped or entangled with the thiolate chains. Their
560 residence in close proximity to the higher temperature environment
561 close to the surface allows them to carry heat away from the surface
562 quite efficiently. There are two aspects of this behavior that are
563 relevant to thermal conductance of the interface: the residence time
564 of solvent molecules in the thiolate layer, and the orientational
565 ordering of the C-C chains as a mechanism for transferring vibrational
566 energy to these entrapped solvent molecules. These competing effects are responsible for the non-monotonic dependence of the interfacial thermal conductance on the fraction of long chain capping agents.
567
568 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
569 % RESIDENCE TIME
570 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
571 \subsection{Mobility of solvent in the interfacial layer}
572
573 We use a simple survival correlation function, $C(t)$, to quantify the
574 residence time of a solvent molecule in the long thiolate chain
575 layer. This function correlates the identity of all hexane molecules
576 within the $z$-coordinate range of the thiolate layer at two separate
577 times. If the solvent molecule is present at both times, the
578 configuration contributes a $1$, while the absence of the molecule at
579 the later time indicates that the solvent molecule has migrated into
580 the bulk, and this configuration contributes a $0$. A steep decay in
581 $C(t)$ indicates a high turnover rate of solvent molecules from the
582 chain region to the bulk. We define the escape rate for trapped solvent molecules at the interface as
583 \begin{equation}
584 k_{escape} = \left( \int_0^T C(t) dt \right)^{-1}
585 \label{eq:mobility}
586 \end{equation}
587 where T is the length of the simulation. This is a direct measure of
588 the rate at which solvent molecules entangled in the thiolate layer
589 can escape into the bulk. As $k_{escape} \rightarrow 0$, the
590 solvent has become permanently trapped in the thiolate layer. In
591 figure \ref{fig:Gstacks} we show that interfacial solvent mobility
592 decreases as the percentage of long thiolate chains increases.
593
594 \begin{figure}
595 \includegraphics[width=\linewidth]{figures/timelapse}
596 \caption{Images of 25\%/75\% C$_4$/C$_{12}$ (top panel) and 75\%/25\% C$_4$/C$_{12}$ (bottom panel) systems at the beginning and end of a 3 ns simulation. Initial interfacial solvent molecules are colored light blue. The solvent-thiolate disorder and interfacial solvent turnover of the 25\%/75\% short/long system stand in stark contrast to the highly ordered but entrapped interfacial solvent of the 75\%/25\% system.}
597 \label{fig:timelapse}
598 \end{figure}
599
600 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
601 % ORDER PARAMETER
602 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
603
604 \subsection{Vibrational coupling via orientational ordering}
605
606 As the fraction of long-chain thiolates becomes large, the entrapped
607 solvent molecules must find specific orientations relative to the mean
608 orientation of the thiolate chains. This configuration allows for
609 efficient thermal energy exchange between the thiolate alkyl chain and
610 the solvent molecules.
611
612 To quantify this cooperative
613 ordering, we computed the orientational order parameters and director
614 axes for both the thiolate chains and for the entrapped solvent. The
615 director axis can be easily obtained by diagonalization of the order
616 parameter tensor,
617 \begin{equation}
618 \mathsf{Q}_{\alpha \beta} = \frac{1}{2 N} \sum_{i=1}^{N} \left( 3 \mathbf{e}_{i
619 \alpha} \mathbf{e}_{i \beta} - \delta_{\alpha \beta} \right)
620 \end{equation}
621 where $\mathbf{e}_{i \alpha}$ was the $\alpha = x,y,z$ component of
622 the unit vector $\mathbf{e}_{i}$ along the long axis of molecule $i$.
623 For both kinds of molecules, the $\mathbf{e}$ vector is defined using
624 the terminal atoms of the chains.
625
626 The largest eigenvalue of $\overleftrightarrow{\mathsf{Q}}$ is
627 traditionally used to obtain orientational order parameter, while the
628 eigenvector corresponding to the order parameter yields the director
629 axis ($\mathbf{d}(t)$) which defines the average direction of
630 molecular alignment at any time $t$. The overlap between the director
631 axes of the thiolates and the entrapped solvent is time-averaged,
632 \begin{equation}
633 \bar{d} = \langle \mathbf{d}_{thiolates} \left( t \right) \cdot
634 \mathbf{d}_{solvent} \left( t \right) \rangle_t
635 \label{eq:orientation}
636 \end{equation}
637 and reported in figure \ref{fig:Gstacks}.
638
639 Once the solvent molecules have picked up thermal energy from the
640 thiolates, they carry heat away from the gold as they diffuse back
641 into the bulk solvent. When the percentage of long chains decreases,
642 the tails of the long chains are much more disordered and do not
643 provide structured channels for the solvent to fill.
644
645 Although the alignment of the chains with the entrapped solvent is one
646 possible mechanism for the non-monotonic increase in the conductance
647 as a function
648
649
650 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
651 % **CONCLUSIONS**
652 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
653 \section{Conclusions}
654
655
656
657 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
658 % **ACKNOWLEDGMENTS**
659 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
660 \section*{Acknowledgments}
661 Support for this project was provided by the
662 National Science Foundation under grant CHE-0848243. Computational
663 time was provided by the Center for Research Computing (CRC) at the
664 University of Notre Dame.
665
666 \newpage
667
668 \bibliography{thiolsRNEMD}
669
670 \end{doublespace}
671 \end{document}
672
673
674
675
676
677
678
679

Properties

Name Value
svn:executable *