| 22 |
|
\setlength{\abovecaptionskip}{20 pt} |
| 23 |
|
\setlength{\belowcaptionskip}{30 pt} |
| 24 |
|
|
| 25 |
< |
%\renewcommand\citemid{\ } % no comma in optional referenc note |
| 25 |
> |
%\renewcommand\citemid{\ } % no comma in optional reference note |
| 26 |
|
\bibpunct{[}{]}{,}{s}{}{;} |
| 27 |
|
\bibliographystyle{aip} |
| 28 |
|
|
| 71 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 72 |
|
|
| 73 |
|
\section{Introduction} |
| 74 |
< |
|
| 74 |
> |
[BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM] |
| 75 |
|
Interfacial thermal conductance is extensively studied both |
| 76 |
|
experimentally and computationally, and systems with interfaces |
| 77 |
|
present are generally heterogeneous. Although interfaces are commonly |
| 97 |
|
|
| 98 |
|
\section{Methodology} |
| 99 |
|
\subsection{Algorithm} |
| 100 |
+ |
[BACKGROUND FOR MD METHODS] |
| 101 |
|
There have been many algorithms for computing thermal conductivity |
| 102 |
|
using molecular dynamics simulations. However, interfacial conductance |
| 103 |
|
is at least an order of magnitude smaller. This would make the |
| 132 |
|
algorithm conserves momenta and energy and does not depend on an |
| 133 |
|
external thermostat. |
| 134 |
|
|
| 135 |
< |
(wondering how much detail of algorithm should be put here...) |
| 135 |
> |
\subsection{Defining Interfacial Thermal Conductivity $G$} |
| 136 |
> |
For interfaces with a relatively low interfacial conductance, the bulk |
| 137 |
> |
regions on either side of an interface rapidly come to a state in |
| 138 |
> |
which the two phases have relatively homogeneous (but distinct) |
| 139 |
> |
temperatures. The interfacial thermal conductivity $G$ can therefore |
| 140 |
> |
be approximated as: |
| 141 |
> |
\begin{equation} |
| 142 |
> |
G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle - |
| 143 |
> |
\langle T_\mathrm{cold}\rangle \right)} |
| 144 |
> |
\label{lowG} |
| 145 |
> |
\end{equation} |
| 146 |
> |
where ${E_{total}}$ is the imposed non-physical kinetic energy |
| 147 |
> |
transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle |
| 148 |
> |
T_\mathrm{cold}\rangle}$ are the average observed temperature of the |
| 149 |
> |
two separated phases. |
| 150 |
|
|
| 151 |
< |
\subsection{Force Field Parameters} |
| 152 |
< |
Our simulation systems consists of metal gold lattice slab solvated by |
| 138 |
< |
organic solvents. In order to study the role of capping agents in |
| 139 |
< |
interfacial thermal conductance, butanethiol is chosen to cover gold |
| 140 |
< |
surfaces in comparison to no capping agent present. |
| 151 |
> |
When the interfacial conductance is {\it not} small, two ways can be |
| 152 |
> |
used to define $G$. |
| 153 |
|
|
| 154 |
< |
The Au-Au interactions in metal lattice slab is described by the |
| 155 |
< |
quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC |
| 156 |
< |
potentials include zero-point quantum corrections and are |
| 157 |
< |
reparametrized for accurate surface energies compared to the |
| 146 |
< |
Sutton-Chen potentials\cite{Chen90}. |
| 147 |
< |
|
| 148 |
< |
Straight chain {\it n}-hexane and aromatic toluene are respectively |
| 149 |
< |
used as solvents. For hexane, both United-Atom\cite{TraPPE-UA.alkanes} |
| 150 |
< |
and All-Atom\cite{OPLSAA} force fields are used for comparison; for |
| 151 |
< |
toluene, United-Atom\cite{TraPPE-UA.alkylbenzenes} force fields are |
| 152 |
< |
used with rigid body constraints applied. (maybe needs more details |
| 153 |
< |
about rigid body) |
| 154 |
< |
|
| 155 |
< |
Buatnethiol molecules are used as capping agent for some of our |
| 156 |
< |
simulations. United-Atom\cite{TraPPE-UA.thiols} and All-Atom models |
| 157 |
< |
are respectively used corresponding to the force field type of |
| 158 |
< |
solvent. |
| 159 |
< |
|
| 160 |
< |
To describe the interactions between metal Au and non-metal capping |
| 161 |
< |
agent and solvent, we refer to Vlugt\cite{vlugt:cpc2007154} and derive |
| 162 |
< |
other interactions which are not parametrized in their work. (can add |
| 163 |
< |
hautman and klein's paper here and more discussion; need to put |
| 164 |
< |
aromatic-metal interaction approximation here) |
| 165 |
< |
|
| 166 |
< |
[TABULATED FORCE FIELD PARAMETERS NEEDED] |
| 167 |
< |
|
| 168 |
< |
\section{Computational Details} |
| 169 |
< |
\subsection{System Geometry} |
| 170 |
< |
Our simulation systems consists of a lattice Au slab with the (111) |
| 171 |
< |
surface perpendicular to the $z$-axis, and a solvent layer between the |
| 172 |
< |
periodic Au slabs along the $z$-axis. To set up the interfacial |
| 173 |
< |
system, the Au slab is first equilibrated without solvent under room |
| 174 |
< |
pressure and a desired temperature. After the metal slab is |
| 175 |
< |
equilibrated, United-Atom or All-Atom butanethiols are replicated on |
| 176 |
< |
the Au surface, each occupying the (??) among three Au atoms, and is |
| 177 |
< |
equilibrated under NVT ensemble. According to (CITATION), the maximal |
| 178 |
< |
thiol capacity on Au surface is $1/3$ of the total number of surface |
| 179 |
< |
Au atoms. |
| 180 |
< |
|
| 181 |
< |
\cite{packmol} |
| 182 |
< |
|
| 183 |
< |
\subsection{Simulation Parameters} |
| 184 |
< |
|
| 185 |
< |
When the interfacial conductance is {\it not} small, there are two |
| 186 |
< |
ways to define $G$. If we assume the temperature is discretely |
| 187 |
< |
different on two sides of the interface, $G$ can be calculated with |
| 188 |
< |
the thermal flux applied $J$ and the temperature difference measured |
| 189 |
< |
$\Delta T$ as: |
| 154 |
> |
One way is to assume the temperature is discretely different on two |
| 155 |
> |
sides of the interface, $G$ can be calculated with the thermal flux |
| 156 |
> |
applied $J$ and the maximum temperature difference measured along the |
| 157 |
> |
thermal gradient max($\Delta T$), which occurs at the interface, as: |
| 158 |
|
\begin{equation} |
| 159 |
|
G=\frac{J}{\Delta T} |
| 160 |
|
\label{discreteG} |
| 161 |
|
\end{equation} |
| 162 |
< |
We can as well assume a continuous temperature profile along the |
| 163 |
< |
thermal gradient axis $z$ and define $G$ as the change of bulk thermal |
| 164 |
< |
conductivity $\lambda$ at a defined interfacial point: |
| 162 |
> |
|
| 163 |
> |
The other approach is to assume a continuous temperature profile along |
| 164 |
> |
the thermal gradient axis (e.g. $z$) and define $G$ at the point where |
| 165 |
> |
the magnitude of thermal conductivity $\lambda$ change reach its |
| 166 |
> |
maximum, given that $\lambda$ is well-defined throughout the space: |
| 167 |
|
\begin{equation} |
| 168 |
|
G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big| |
| 169 |
|
= \Big|\frac{\partial}{\partial z}\left(-J_z\Big/ |
| 170 |
|
\left(\frac{\partial T}{\partial z}\right)\right)\Big| |
| 171 |
< |
= J_z\Big|\frac{\partial^2 T}{\partial z^2}\Big| |
| 171 |
> |
= |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big| |
| 172 |
|
\Big/\left(\frac{\partial T}{\partial z}\right)^2 |
| 173 |
|
\label{derivativeG} |
| 174 |
|
\end{equation} |
| 175 |
+ |
|
| 176 |
|
With the temperature profile obtained from simulations, one is able to |
| 177 |
|
approximate the first and second derivatives of $T$ with finite |
| 178 |
|
difference method and thus calculate $G^\prime$. |
| 179 |
|
|
| 180 |
|
In what follows, both definitions are used for calculation and comparison. |
| 181 |
|
|
| 182 |
+ |
[IMPOSE G DEFINITION INTO OUR SYSTEMS] |
| 183 |
+ |
To facilitate the use of the above definitions in calculating $G$ and |
| 184 |
+ |
$G^\prime$, we have a metal slab with its (111) surfaces perpendicular |
| 185 |
+ |
to the $z$-axis of our simulation cells. With or withour capping |
| 186 |
+ |
agents on the surfaces, the metal slab is solvated with organic |
| 187 |
+ |
solvents, as illustrated in Figure \ref{demoPic}. |
| 188 |
+ |
|
| 189 |
+ |
\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} |
| 208 |
+ |
\includegraphics[width=\linewidth]{gradT} |
| 209 |
+ |
\caption{The 1st and 2nd derivatives of temperature profile can be |
| 210 |
+ |
obtained with finite difference approximation.} |
| 211 |
+ |
\label{gradT} |
| 212 |
+ |
\end{figure} |
| 213 |
+ |
|
| 214 |
+ |
\section{Computational Details} |
| 215 |
+ |
\subsection{System Geometry} |
| 216 |
+ |
In our simulations, Au is used to construct a metal slab with bare |
| 217 |
+ |
(111) surface perpendicular to the $z$-axis. Different slab thickness |
| 218 |
+ |
(layer numbers of Au) are simulated. This metal slab is first |
| 219 |
+ |
equilibrated under normal pressure (1 atm) and a desired |
| 220 |
+ |
temperature. After equilibration, butanethiol is used as the capping |
| 221 |
+ |
agent molecule to cover the bare Au (111) surfaces evenly. The sulfur |
| 222 |
+ |
atoms in the butanethiol molecules would occupy the three-fold sites |
| 223 |
+ |
of the surfaces, and the maximal butanethiol capacity on Au surface is |
| 224 |
+ |
$1/3$ of the total number of surface Au atoms[CITATION]. A series of |
| 225 |
+ |
different coverage surfaces is investigated in order to study the |
| 226 |
+ |
relation between coverage and conductance. |
| 227 |
+ |
|
| 228 |
+ |
[COVERAGE DISCRIPTION] However, since the interactions between surface |
| 229 |
+ |
Au and butanethiol is non-bonded, the capping agent molecules are |
| 230 |
+ |
allowed to migrate to an empty neighbor three-fold site during a |
| 231 |
+ |
simulation. Therefore, the initial configuration would not severely |
| 232 |
+ |
affect the sampling of a variety of configurations of the same |
| 233 |
+ |
coverage, and the final conductance measurement would be an average |
| 234 |
+ |
effect of these configurations explored in the simulations. [MAY NEED FIGURES] |
| 235 |
+ |
|
| 236 |
+ |
After the modified Au-butanethiol surface systems are equilibrated |
| 237 |
+ |
under canonical ensemble, Packmol\cite{packmol} is used to pack |
| 238 |
+ |
organic solvent molecules in the previously vacuum part of the |
| 239 |
+ |
simulation cells, which guarantees that short range repulsive |
| 240 |
+ |
interactions do not disrupt the simulations. Two solvents are |
| 241 |
+ |
investigated, one which has little vibrational overlap with the |
| 242 |
+ |
alkanethiol and plane-like shape (toluene), and one which has similar |
| 243 |
+ |
vibrational frequencies and chain-like shape ({\it n}-hexane). The |
| 244 |
+ |
spacing filled by solvent molecules, i.e. the gap between periodically |
| 245 |
+ |
repeated Au-butanethiol surfaces should be carefully chosen so that it |
| 246 |
+ |
would not be too short to affect the liquid phase structure, nor too |
| 247 |
+ |
long, leading to over cooling (freezing) or heating (boiling) when a |
| 248 |
+ |
thermal flux is applied. In our simulations, this spacing is usually |
| 249 |
+ |
$35 \sim 60$\AA. |
| 250 |
+ |
|
| 251 |
+ |
The initial configurations generated by Packmol are further |
| 252 |
+ |
equilibrated with the $x$ and $y$ dimensions fixed, only allowing |
| 253 |
+ |
length scale change in $z$ dimension. This is to ensure that the |
| 254 |
+ |
equilibration of liquid phase does not affect the metal crystal |
| 255 |
+ |
structure in $x$ and $y$ dimensions. Further equilibration are run |
| 256 |
+ |
under NVT and then NVE ensembles. |
| 257 |
+ |
|
| 258 |
+ |
After the systems reach equilibrium, NIVS is implemented to impose a |
| 259 |
+ |
periodic unphysical thermal flux between the metal and the liquid |
| 260 |
+ |
phase. Most of our simulations are under an average temperature of |
| 261 |
+ |
$\sim$200K. Therefore, this flux usually comes from the metal to the |
| 262 |
+ |
liquid so that the liquid has a higher temperature and would not |
| 263 |
+ |
freeze due to excessively low temperature. This induced temperature |
| 264 |
+ |
gradient is stablized and the simulation cell is devided evenly into |
| 265 |
+ |
N slabs along the $z$-axis and the temperatures of each slab are |
| 266 |
+ |
recorded. When the slab width $d$ of each slab is the same, the |
| 267 |
+ |
derivatives of $T$ with respect to slab number $n$ can be directly |
| 268 |
+ |
used for $G^\prime$ calculations: |
| 269 |
+ |
\begin{equation} |
| 270 |
+ |
G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big| |
| 271 |
+ |
\Big/\left(\frac{\partial T}{\partial z}\right)^2 |
| 272 |
+ |
= |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big| |
| 273 |
+ |
\Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2 |
| 274 |
+ |
= |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big| |
| 275 |
+ |
\Big/\left(\frac{\partial T}{\partial n}\right)^2 |
| 276 |
+ |
\label{derivativeG2} |
| 277 |
+ |
\end{equation} |
| 278 |
+ |
|
| 279 |
+ |
\subsection{Force Field Parameters} |
| 280 |
+ |
Our simulations include various components. Therefore, force field |
| 281 |
+ |
parameter descriptions are needed for interactions both between the |
| 282 |
+ |
same type of particles and between particles of different species. |
| 283 |
+ |
|
| 284 |
+ |
The Au-Au interactions in metal lattice slab is described by the |
| 285 |
+ |
quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC |
| 286 |
+ |
potentials include zero-point quantum corrections and are |
| 287 |
+ |
reparametrized for accurate surface energies compared to the |
| 288 |
+ |
Sutton-Chen potentials\cite{Chen90}. |
| 289 |
+ |
|
| 290 |
+ |
For both solvent molecules, straight chain {\it n}-hexane and aromatic |
| 291 |
+ |
toluene, United-Atom (UA) and All-Atom (AA) models are used |
| 292 |
+ |
respectively. The TraPPE-UA |
| 293 |
+ |
parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used |
| 294 |
+ |
for our UA solvent molecules. In these models, pseudo-atoms are |
| 295 |
+ |
located at the carbon centers for alkyl groups. By eliminating |
| 296 |
+ |
explicit hydrogen atoms, these models are simple and computationally |
| 297 |
+ |
efficient, while maintains good accuracy. However, the TraPPE-UA for |
| 298 |
+ |
alkanes is known to predict a lower boiling point than experimental |
| 299 |
+ |
values. Considering that after an unphysical thermal flux is applied |
| 300 |
+ |
to a system, the temperature of ``hot'' area in the liquid phase would be |
| 301 |
+ |
significantly higher than the average, to prevent over heating and |
| 302 |
+ |
boiling of the liquid phase, the average temperature in our |
| 303 |
+ |
simulations should be much lower than the liquid boiling point. [NEED MORE DISCUSSION] |
| 304 |
+ |
For UA-toluene model, rigid body constraints are applied, so that the |
| 305 |
+ |
benzene ring and the methyl-C(aromatic) bond are kept rigid. This |
| 306 |
+ |
would save computational time.[MORE DETAILS NEEDED] |
| 307 |
+ |
|
| 308 |
+ |
Besides the TraPPE-UA models, AA models for both organic solvents are |
| 309 |
+ |
included in our studies as well. For hexane, the OPLS |
| 310 |
+ |
all-atom\cite{OPLSAA} force field is used. [MORE DETAILS] |
| 311 |
+ |
For toluene, the United Force Field developed by Rapp\'{e} {\it et |
| 312 |
+ |
al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS] |
| 313 |
+ |
|
| 314 |
+ |
The capping agent in our simulations, the butanethiol molecules can |
| 315 |
+ |
either use UA or AA model. The TraPPE-UA force fields includes |
| 316 |
+ |
parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used in |
| 317 |
+ |
our simulations corresponding to our TraPPE-UA models for solvent. |
| 318 |
+ |
and All-Atom models [NEED CITATIONS] |
| 319 |
+ |
However, the model choice (UA or AA) of capping agent can be different |
| 320 |
+ |
from the solvent. Regardless of model choice, the force field |
| 321 |
+ |
parameters for interactions between capping agent and solvent can be |
| 322 |
+ |
derived using Lorentz-Berthelot Mixing Rule. |
| 323 |
+ |
|
| 324 |
+ |
To describe the interactions between metal Au and non-metal capping |
| 325 |
+ |
agent and solvent, we refer to Vlugt\cite{vlugt:cpc2007154} and derive |
| 326 |
+ |
other interactions which are not yet finely parametrized. [can add |
| 327 |
+ |
hautman and klein's paper here and more discussion; need to put |
| 328 |
+ |
aromatic-metal interaction approximation here]\cite{doi:10.1021/jp034405s} |
| 329 |
+ |
|
| 330 |
+ |
[TABULATED FORCE FIELD PARAMETERS NEEDED] |
| 331 |
+ |
|
| 332 |
+ |
|
| 333 |
+ |
[SURFACE RECONSTRUCTION PREVENTS SIMULATION TEMP TO GO HIGHER] |
| 334 |
+ |
|
| 335 |
+ |
|
| 336 |
|
\section{Results} |
| 337 |
+ |
[REARRANGEMENT NEEDED] |
| 338 |
|
\subsection{Toluene Solvent} |
| 339 |
|
|
| 340 |
< |
The simulations follow a protocol similar to the previous gold/water |
| 215 |
< |
interfacial systems. The results (Table \ref{AuThiolToluene}) show a |
| 340 |
> |
The results (Table \ref{AuThiolToluene}) show a |
| 341 |
|
significant conductance enhancement compared to the gold/water |
| 342 |
|
interface without capping agent and agree with available experimental |
| 343 |
|
data. This indicates that the metal-metal potential, though not |
