| 52 |
|
|
| 53 |
|
\begin{abstract} |
| 54 |
|
|
| 55 |
< |
We report on simulations of heat conduction from Au(111) / hexane interfaces in which the surface has been protected by a mix of short and long chain alkanethiolate ligands. A variant of reverse non-equilibrium molecular dynamics (RNEMD) was used to create a thermal flux between the metal and solvent, and thermal conductance was computed using the resulting thermal profiles of the interface. We find a non-monotonic dependence of the interfacial thermal conductance on the fraction of long chains present at the interface, and correlate this behavior to solvent ordering and escape from the thiolate layer immediately in contact with the metal surface. |
| 56 |
< |
Our results suggest that a mixed vibrational transfer / convection model is |
| 57 |
< |
necessary to explain the features of heat transfer at this interface. The |
| 58 |
< |
alignment of the solvent chains with the ordered ligand allows rapid transfer |
| 59 |
< |
of energy to the trapped solvent and is the dominant feature for ordered |
| 60 |
< |
ligand layers. Diffusion of the vibrationally excited solvent into the bulk |
| 61 |
< |
also plays a significant role when the ligands are less tightly packed. |
| 55 |
> |
We report on simulations of heat conduction through Au(111) / hexane |
| 56 |
> |
interfaces in which the surface has been protected by a mixture of |
| 57 |
> |
short and long chain alkanethiolate ligands. Reverse |
| 58 |
> |
non-equilibrium molecular dynamics (RNEMD) was used to create a |
| 59 |
> |
thermal flux between the metal and solvent, and thermal conductance |
| 60 |
> |
was computed using the resulting thermal profiles near the |
| 61 |
> |
interface. We find a non-monotonic dependence of the interfacial |
| 62 |
> |
thermal conductance on the fraction of long chains present at the |
| 63 |
> |
interface, and correlate this behavior to both solvent ordering and |
| 64 |
> |
the rate of solvent escape from the thiolate layer immediately in |
| 65 |
> |
contact with the metal surface. Our results suggest that a mixed |
| 66 |
> |
vibrational transfer / convection model is necessary to explain the |
| 67 |
> |
features of heat transfer at this interface. The alignment of the |
| 68 |
> |
solvent chains with the ordered ligand allows rapid transfer of |
| 69 |
> |
energy to the trapped solvent and is the dominant feature for |
| 70 |
> |
ordered ligand layers. Diffusion of the vibrationally excited |
| 71 |
> |
solvent into the bulk also plays a significant role when the ligands |
| 72 |
> |
are less tightly packed. |
| 73 |
|
|
| 74 |
|
\end{abstract} |
| 75 |
|
|
| 83 |
|
\section{Introduction} |
| 84 |
|
|
| 85 |
|
The structural details of the interfaces of metal nanoparticles |
| 86 |
< |
determines how energy flows between these particles and their |
| 86 |
> |
determine how energy flows between these particles and their |
| 87 |
|
surroundings. Understanding this energy flow is essential in designing |
| 88 |
|
and functionalizing metallic nanoparticles for use in plasmonic photothermal |
| 89 |
|
therapies,\cite{Jain:2007ux,Petrova:2007ad,Gnyawali:2008lp,Mazzaglia:2008to,Huff:2007ye,Larson:2007hw} |
| 94 |
|
interfacial thermal conductance, $G$, which can be somewhat difficult |
| 95 |
|
to determine experimentally.\cite{Wilson:2002uq,Plech:2005kx} |
| 96 |
|
|
| 97 |
< |
Metallic particles have also been proposed for use in highly efficient |
| 97 |
> |
Metallic particles have also been proposed for use in efficient |
| 98 |
|
thermal-transfer fluids, although careful experiments by Eapen {\it et |
| 99 |
|
al.} have shown that metal-particle-based nanofluids have thermal |
| 100 |
|
conductivities that match Maxwell predictions.\cite{Eapen:2007th} The |
| 130 |
|
eliminate this barrier |
| 131 |
|
($G\rightarrow\infty$).\cite{doi:10.1021/la904855s} |
| 132 |
|
|
| 133 |
< |
In previous simulations, we applied a variant of reverse |
| 134 |
< |
non-equilibrium molecular dynamics (RNEMD) to calculate the |
| 135 |
< |
interfacial thermal conductance at a metal / organic solvent interface |
| 136 |
< |
that had been chemically protected by butanethiolate groups. Our |
| 133 |
> |
Recently, Hase and coworkers employed Non-Equilibrium Molecular |
| 134 |
> |
Dynamics (NEMD) simulations to study thermal transport from hot |
| 135 |
> |
Au(111) substrate to a self-assembled monolayer of alkylthiol with |
| 136 |
> |
relatively long chain (8-20 carbon atoms).\cite{hase:2010,hase:2011} |
| 137 |
> |
These simulations explained many of the features of the experiments of |
| 138 |
> |
Wang {\it et al.} However, ensemble averaged measurements for heat |
| 139 |
> |
conductance of interfaces between the capping monolayer on Au and a |
| 140 |
> |
solvent phase have yet to be studied with their approach. In previous |
| 141 |
> |
simulations, our group applied a variant of reverse non-equilibrium |
| 142 |
> |
molecular dynamics (RNEMD) to calculate the interfacial thermal |
| 143 |
> |
conductance at a metal / organic solvent interface that had been |
| 144 |
> |
chemically protected by butanethiolate groups.\cite{kuang:AuThl} Our |
| 145 |
|
calculations suggested an explanation for the very large thermal |
| 146 |
|
conductivity at alkanethiol-capped metal surfaces when compared with |
| 147 |
|
bare metal/solvent interfaces. Specifically, the chemical bond |
| 153 |
|
A notable result of the previous simulations was the non-monotonic |
| 154 |
|
dependence of $G$ on the fractional coverage of the metal surface by |
| 155 |
|
the chemical protecting group. Gaps in surface coverage allow the |
| 156 |
< |
solvent molecules to come into direct contact with ligands that have been |
| 157 |
< |
heated by the metal surface, absorb thermal energy from the ligands, |
| 158 |
< |
and then diffuse away. Quantifying the role of surface coverage is |
| 159 |
< |
difficult as the ligands have lateral mobility on the surface and can |
| 160 |
< |
aggregate to form domains on the timescale of the simulation. |
| 156 |
> |
solvent molecules to come into direct contact with ligands that have |
| 157 |
> |
been heated by the metal surface, absorb thermal energy from the |
| 158 |
> |
ligands, and then diffuse away. Quantifying the role of overall |
| 159 |
> |
surface coverage was difficult because the ligands have lateral |
| 160 |
> |
mobility on the surface and can aggregate to form domains on the |
| 161 |
> |
timescale of the simulation. |
| 162 |
|
|
| 163 |
< |
To isolate this effect while avoiding lateral mobility of the |
| 164 |
< |
surface ligands, the current work involves mixed-chain-length |
| 165 |
< |
monolayers in which the length mismatch between long and short chains |
| 166 |
< |
is sufficient to accommodate solvent molecules. These completely |
| 167 |
< |
covered (but mixed-chain) surfaces are significantly less prone to |
| 168 |
< |
surface domain formation, and allow us to further investigate the |
| 169 |
< |
mechanism of heat transport to the solvent. A thermal flux is created |
| 170 |
< |
using velocity shearing and scaling reverse non-equilibrium molecular |
| 171 |
< |
dynamics (VSS-RNEMD), and the resulting temperature profiles are |
| 172 |
< |
analyzed to yield information about the interfacial thermal |
| 173 |
< |
conductance. |
| 163 |
> |
To isolate the effect of ligand/solvent coupling while avoiding |
| 164 |
> |
lateral mobility of the surface ligands, the current work utilizes |
| 165 |
> |
monolayers of mixed chain-lengths in which the length mismatch between |
| 166 |
> |
long and short chains is sufficient to accommodate solvent |
| 167 |
> |
molecules. These completely covered (but mixed-chain) surfaces are |
| 168 |
> |
significantly less prone to surface domain formation, and allow us to |
| 169 |
> |
further investigate the mechanism of heat transport to the solvent. A |
| 170 |
> |
thermal flux is created using velocity shearing and scaling reverse |
| 171 |
> |
non-equilibrium molecular dynamics (VSS-RNEMD), and the resulting |
| 172 |
> |
temperature profiles are analyzed to yield information about the |
| 173 |
> |
interfacial thermal conductance. |
| 174 |
|
|
| 175 |
|
|
| 176 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 180 |
|
|
| 181 |
|
There are many ways to compute bulk transport properties from |
| 182 |
|
classical molecular dynamics simulations. Equilibrium Molecular |
| 183 |
< |
Dynamics (EMD) simulations can be used by computing relevant time |
| 184 |
< |
correlation functions and assuming that linear response theory |
| 183 |
> |
Dynamics (EMD) simulations can be used to compute the relevant time |
| 184 |
> |
correlation functions and transport coefficients can be calculated |
| 185 |
> |
assuming that linear response theory |
| 186 |
|
holds.\cite{PhysRevB.37.5677,MASSOBRIO:1984bl,PhysRev.119.1,Viscardy:2007rp,che:6888,kinaci:014106} |
| 187 |
|
For some transport properties (notably thermal conductivity), EMD |
| 188 |
|
approaches are subject to noise and poor convergence of the relevant |
| 233 |
|
large.\cite{Maginn:2010} |
| 234 |
|
|
| 235 |
|
The most useful alternative RNEMD approach developed so far utilizes a |
| 236 |
< |
series of simultaneous velocity shearing and scaling exchanges between |
| 236 |
> |
series of simultaneous velocity shearing and scaling (VSS) exchanges between |
| 237 |
|
the two slabs.\cite{Kuang2012} This method provides a set of |
| 238 |
|
conservation constraints while simultaneously creating a desired flux |
| 239 |
|
between the two slabs. Satisfying the constraint equations ensures |
| 260 |
|
the center of mass velocities in the $C$ and $H$ slabs, respectively. |
| 261 |
|
Within the two slabs, particles receive incremental changes or a |
| 262 |
|
``shear'' to their velocities. The amount of shear is governed by the |
| 263 |
< |
imposed momentum flux, $j_z(\mathbf{p})$ |
| 263 |
> |
imposed momentum flux, $\mathbf{j}_z(\mathbf{p})$ |
| 264 |
|
\begin{eqnarray} |
| 265 |
|
\mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\ |
| 266 |
|
\mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2} |
| 267 |
|
\end{eqnarray} |
| 268 |
< |
where $M_{\{c,h\}}$ is the total mass of particles within each slab and $\Delta t$ |
| 269 |
< |
is the interval between two separate operations. |
| 268 |
> |
where $M_{\{c,h\}}$ is the total mass of particles within each of the |
| 269 |
> |
slabs and $\Delta t$ is the interval between two separate operations. |
| 270 |
|
|
| 271 |
|
To simultaneously impose a thermal flux ($J_z$) between the slabs we |
| 272 |
|
use energy conservation constraints, |
| 287 |
|
$j_z(\mathbf{p})$, and thermal flux, $J_z$, values. Since the VSS |
| 288 |
|
operations do not change the kinetic energy due to orientational |
| 289 |
|
degrees of freedom or the potential energy of a system, configurations |
| 290 |
< |
after the VSS move have exactly the same energy ({\it and} linear |
| 290 |
> |
after the VSS move have exactly the same energy (and linear |
| 291 |
|
momentum) as before the move. |
| 292 |
|
|
| 293 |
|
As the simulation progresses, the VSS moves are performed on a regular |
| 309 |
|
The VSS-RNEMD approach is versatile in that it may be used to |
| 310 |
|
implement both thermal and shear transport either separately or |
| 311 |
|
simultaneously. Perturbations of velocities away from the ideal |
| 312 |
< |
Maxwell-Boltzmann distributions are minimal, as is thermal anisotropy. This ability to generate simultaneous thermal and |
| 313 |
< |
shear fluxes has been previously utilized to map out the shear |
| 314 |
< |
viscosity of SPC/E water over a wide range of temperatures (90~K) with |
| 315 |
< |
a single 1 ns simulation.\cite{Kuang2012} |
| 312 |
> |
Maxwell-Boltzmann distributions are minimal, as is thermal anisotropy. |
| 313 |
> |
This ability to generate simultaneous thermal and shear fluxes has |
| 314 |
> |
been previously utilized to map out the shear viscosity of SPC/E water |
| 315 |
> |
over a wide range of temperatures (90~K) with a single 1 ns |
| 316 |
> |
simulation.\cite{Kuang2012} |
| 317 |
|
|
| 318 |
|
\begin{figure} |
| 319 |
|
\includegraphics[width=\linewidth]{figures/rnemd} |
| 351 |
|
or thermal flux and watching the gradient response of the |
| 352 |
|
material. This means that the {\it interfacial} conductance is a |
| 353 |
|
second derivative property which is subject to significant noise and |
| 354 |
< |
therefore requires extensive sampling. We have been able to |
| 355 |
< |
demonstrate the use of the second derivative approach to compute |
| 356 |
< |
interfacial conductance at chemically-modified metal / solvent |
| 357 |
< |
interfaces. However, a definition of the interfacial conductance in |
| 358 |
< |
terms of a temperature difference ($\Delta T$) measured at $z_0$, |
| 354 |
> |
therefore requires extensive sampling. Previous work has demonstrated |
| 355 |
> |
the use of the second derivative approach to compute interfacial |
| 356 |
> |
conductance at chemically-modified metal / solvent interfaces. |
| 357 |
> |
However, a definition of the interfacial conductance in terms of a |
| 358 |
> |
temperature difference ($\Delta T$) measured at $z_0$, |
| 359 |
|
\begin{displaymath} |
| 360 |
|
G = \frac{J_z}{\Delta T_{z_0}}, |
| 361 |
|
\end{displaymath} |
| 472 |
|
|
| 473 |
|
We have implemented the VSS-RNEMD algorithm in OpenMD, our group |
| 474 |
|
molecular dynamics code.\cite{openmd} An 1188 atom gold slab was |
| 475 |
< |
equilibrated at 1 atm and 200 K. The periodic box was then expanded |
| 476 |
< |
in the $z$ direction to expose two Au(111) faces on either side of the |
| 477 |
< |
11-atom thick slab. |
| 475 |
> |
prepared and equilibrated at 1 atm and 200 K. The periodic box was |
| 476 |
> |
then expanded in the $z$ direction to expose two Au(111) faces on |
| 477 |
> |
either side of the 11-layer slab. |
| 478 |
|
|
| 479 |
|
A full monolayer of thiolates (1/3 the number of surface gold atoms) |
| 480 |
|
was placed on three-fold hollow sites on the gold interfaces. The |
| 495 |
|
|
| 496 |
|
The simulation box $z$-dimension was set to roughly double the length |
| 497 |
|
of the gold / thiolate block. Hexane solvent molecules were placed in |
| 498 |
< |
the vacant portion of the box using the packmol algorithm. Figure \ref{fig:timelapse} shows two of the mixed chain length systems both before and after the RNEMD simulation. |
| 498 |
> |
the vacant portion of the box using the packmol |
| 499 |
> |
algorithm.\cite{packmol} Figure \ref{fig:timelapse} shows two of the |
| 500 |
> |
mixed chain length interfaces both before and after the RNEMD simulation. |
| 501 |
|
|
| 502 |
|
\begin{figure} |
| 503 |
|
\includegraphics[width=\linewidth]{figures/timelapse} |
| 521 |
|
ensemble before proceeding with the VSS-RNEMD and data collection |
| 522 |
|
stages. |
| 523 |
|
|
| 524 |
< |
A kinetic energy flux was applied using VSS-RNEMD. The total |
| 525 |
< |
simulation time was 3 nanoseconds, with velocity scaling occurring |
| 526 |
< |
every 10 femtoseconds. The ``hot'' slab was centered in the gold and the |
| 527 |
< |
``cold'' slab was placed in the center of the solvent region. The entire |
| 528 |
< |
system had a (time-averaged) temperature of 220 K, with a temperature |
| 529 |
< |
difference between the hot and cold slabs of approximately 30 K. The |
| 530 |
< |
average temperature and kinetic energy flux were selected to avoid |
| 531 |
< |
solvent freezing (or glass formation) and to prevent the thiolates from |
| 532 |
< |
burying in the gold slab. The Au-S interaction has a deep potential |
| 533 |
< |
energy well, which allows sulfur atoms to embed into the gold slab at |
| 534 |
< |
temperatures above 250 K. Simulation conditions were chosen to avoid |
| 535 |
< |
both of these undesirable situations. |
| 524 |
> |
A kinetic energy flux was applied using VSS-RNEMD during a data |
| 525 |
> |
collection period of 3 nanoseconds, with velocity scaling moves |
| 526 |
> |
occurring every 10 femtoseconds. The ``hot'' slab was centered in the |
| 527 |
> |
gold and the ``cold'' slab was placed in the center of the solvent |
| 528 |
> |
region. The entire system had a (time-averaged) temperature of 220 K, |
| 529 |
> |
with a temperature difference between the hot and cold slabs of |
| 530 |
> |
approximately 30 K. The average temperature and kinetic energy flux |
| 531 |
> |
were selected to avoid solvent freezing (or glass formation) and to |
| 532 |
> |
prevent the thiolates from burying in the gold slab. The Au-S |
| 533 |
> |
interaction has a deep potential energy well, which allows sulfur |
| 534 |
> |
atoms to embed into the gold slab at temperatures above 250 K. |
| 535 |
> |
Simulation conditions were chosen to avoid both of these |
| 536 |
> |
situations. |
| 537 |
|
|
| 538 |
|
Temperature profiles of the system were created by dividing the box |
| 539 |
|
into $\sim$ 3 \AA \, bins along the $z$ axis and recording the average |
| 681 |
|
% **CONCLUSIONS** |
| 682 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 683 |
|
\section{Conclusions} |
| 684 |
< |
Our results suggest that a mixed vibrational transfer / convection model is necessary to explain the features of heat transfer at this interface. The alignment of the solvent chains with the ordered ligand allows rapid transfer of energy to the trapped solvent and becomes the dominant feature for ordered ligand layers. Diffusion of the vibrationally excited solvent into the bulk also plays a significant role when the ligands are less tightly packed. |
| 684 |
> |
Our results suggest that a mixed vibrational transfer / convection |
| 685 |
> |
model may be necessary to explain the features of heat transfer at |
| 686 |
> |
this interface. The alignment of the solvent chains with the ordered |
| 687 |
> |
ligand allows rapid transfer of energy to the trapped solvent and |
| 688 |
> |
becomes the dominant feature for ordered ligand layers. Diffusion of |
| 689 |
> |
the vibrationally excited solvent into the bulk also plays a |
| 690 |
> |
significant role when the ligands are less tightly packed. |
| 691 |
|
|
| 692 |
|
In the language of earlier continuum approaches to interfacial |
| 693 |
|
conductance,\cite{RevModPhys.61.605} the alignment of the chains is an |
