PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 5108 K:/Prolecc&Pres/ClRC/DateDue.indd MOLECULAR DYNAMICS SIMULATION OF THERMAL ENERGY TRANSPORT ACROSS MATERIAL INTERFACES By Tengfei Luo A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY MECHANICAL ENGINEERING 2009 MOLECUI 801 fit them Tnernal simulated IEMD) ar {self-ass mfliecule Gwen-Kl hem EM W si coverag Conducu Abstract MOLECULAR DYNAMICS SIMULATION OF THERMAL ENERGY TRANPSORT ACROSS MATERIAL INTERFACES By Tengfei Luo Both ab-initio and classical molecular dynamics (MD) were used to study the thermal energy transport phenomena across nano-scele material interfaces. Thermal equilibration in semiconductor ultra-thin layered superlattices was simulated by ab-initio MD with density functional theory (DFT). Equilibrium MD (EMD) and Non-equilibrium MD (NEMD) simulations were performed on Au-SAM (self-assembly monolayer)-Au junctions with alkanedithiols being the SAM ‘ molecules. The in-plane thermal conductivities were calculated using EMD with Green-Kubo method. The out-of-plane thermal conductances were calculated in both EMD and NEMD simulations. Au substrate thickness effect, temperature effect, simulated normal pressure effect, molecular chain length effect, molecule coverage effect and molecule-substrate bonding strength effect on thermal conductivity/conductance were studied. Vibration density of states (VDOS) was calculated, and the mechanism of thermal energy transport across the material junctions was analyzed. The calculated thermal conductance at high temperatures agrees well with available experimental data. The temperature dependence of thermal conductance has a similar trend to experimental 3958M itefiac slim am SF lamina was us As-S it were fl simulat energy tamed lileria: W81: Wu: mom a'G'oss West Midi Elem: observations. SAM molecular coverage was found to be important on the interfacial thermal conductance. Analysis of the junction response to a heat pulse showed that the Au-SAM interface resistance was much larger than the substrate and SAM resistances. The results showed that the Au-SAM interface resistance dominated thermal energy transport across the junction The DFT ab-initio method was used to study the bondings of thiols on As-terminated GaAs (001) surfaces. As-S interactions were simulated by the Morse potential, and the parameters were fitted to an energy hypersurface obtained from DFT calculations. NEMD simulations were then performed on GaAs-SAM-GaAs junctions to study thermal energy transport across thioI-GaAs interfaces. NEMD simulations were also carried out to study thermal energy transport across different graphene-polymer interfaces. The results of this study will be useful for the current molecular electronics ‘ industry in which thermal dissipation is a critical problem to be resolved. It is concluded that the interfacial resistance is the barrier for thermal transport across molecule-solid junctions. As a result, methods to facilitate thermal transport across the interfaces, such as depositing denser SAM, forming stronger molecule—solid bonds, choosing materials with better vibration coupling, are to be considered in the emerging technology of the manufacturing of molecular electronics. m JfS them the fi Acknowledgements At this moment of finishing the dissertation, I would like to thank my advisor, Dr. John R. Lloyd for his exceptional and patient guidance, warm-hearted encouragement, and firm support. His broadness and depth of knowledge, and passion to explore the unknown fields embodied a spirit to seek excellence. I particularly thank Dr. Lloyd for giving me the chance to participate in this exciting project, and his dedication to my academic development. I thank professors in my Ph. D guidance committee: Dr. S. D. Mahanti, whose in-depth understanding of phonon transport and DFT help me to understand thermal energy transport from a microscopic view; Dr. N. Priejzev, from whom I learnt so much of molecular modeling; Dr. F. Jaberi, whose heat conduction course lead to me the world of the heat transfer; Dr. A. Engeda, who taught me thermodynamics - the foundation of thermal physics. I gratefully acknowledge national science foundation award ID 0522594 for the financial support. I want to express my sincere appreciation to my wife, Yuanyuan. Her understanding and confidence in me support me to work. hard with perseverance and self-discipline. Finally, I wish to especially thank my parents for their constant love and encouragement. iv LIST OF llST 0F Homeric CWT] 1.1 l 1.2 1.3 CHAPT 2.1 . Table of Contents LIST OF TABLES ................................................................................................ viii LIST OF FIGURES ............................................................................................... x Nomenclature ..................................................................................................... xiv CHAPTER 1 INTRODUCTION ............................................................................. 1 1.1 Motivation ................................................................................................. 1 1.2 Research Background .............................................................................. 6 1.3 Specific Research Problems and Objectives ............................................ 8 CHAPTER 2 THEORIES AND METHODOLOGIES ............................................ 12 2.1 Ab—initio Molecular Dynamics .................................................................. 12 2.1.1 Bom-Oppenheimer Molecular Dynamics ....................................... 12 2.1.2 Car-Parrinello Molecular Dynamics (CPMD) ................................. 16 2.1.3 Density Functional Theory ............................................................. 21 2.1.4 Plane Wave Basis ......................................................................... 22 2.1.5 Pseudopotentials ........................................................................... 23 2.1.6 Limits of ab—initio MD ..................................................................... 24 2.2 Classical Molecular Dynamics ................................................................ 24 2.3 Molecular Dynamics Simulation in Then'nal Energy Transport ................ 26 CHAPTER 3 AB-INITIO MOLECULAR DYNAMICS SIMULATION OF THERMAL ENERGY TRANSPORT ACROSS MATERIAL INTERFACES ............................ 32 3.1 Backgromd and Objective ...................................................................... 32 3.2 System and Simulation ........................................................................... 34 3.3 Results and Discussion .......................................................................... 41 3.3.1 64-atom, Silicon-Silicon ................................................................. 41 3.3.2 64-atom, Gennanium-Germanium ................................................. 52 3.3.3 64-atom, Silicon-Germanium ......................................................... 54 3.4 Summary and Conclusion ....................................................................... 66 CHAPTER 4 Equilibrium Molecular Dynamics Calculation of Thermal Conductivity/conductance of Au—SAM-Au Junctions ........................................... 69 4.1 Introduction of Molecular Electronics ...................................................... 69 4.2 Equilibrium Molecular Dynamics and Green-Kubo Method .................... 73 4.3 System, Simulation and Potential ........................................................... 74 4.4 I 4.5 CHAPT THERII 5.1 5.2 4.4 Results and Discussion .......................................................................... 87 4.4.1 Defining the Value of Thermal Conductivity ................................... 87 4.4.2 Boundary Conditions ..................................................................... 97 4.4.3 Finite Size Effect ............................................................................ 97 4.4.4 Substrate Thickness Effect ............................................................ 98 4.4.5 Temperature Effect ...................................................................... 102 4.4.6 Simulated Normal Pressure Effect ............................................... 108 4.4.7 Chain Length Effect ..................................................................... 110 4.4.8 SAM Molecule Coverage Effect ................................................... 112 4.4.9 Vibration Coupling Analysis ......................................................... 113 4.5 Summary and Conclusion ..................................................................... 115 CHAPTER 5 NON-EQUILIBRIUM MOLECULAR DYNAMICS SIMULATION OF THERMAL ENERGY TRANSPORT IN AU-SAM-AU JUNCTIONS ................... 118 5.1 Non-Equilibrium Molecular Dynamics and Simulation System ............. 118 5.2 Results and Discussion ........................................................................ 121 5.2.1 Finite Size Effect ......................... . ................................................. 121 5.2.2 Temperature Effect ...................................................................... 134 5.2.3 Simulated Pressure Effect ........................................................... 142 5.2.4 Coverage Effect ........................................................................... 146 5.2.5 Au-SAM Bond Strength Effect ...................................................... 147 5.2.6 Heat pulse ................................................................................... 152 5.2.7 Further VDOS Analysis ................................................................ 155 5.3 Summary and Conclusion ..................................................................... 159 CHAPTER 6 DENSITY FUNCTIONAL THEORY STUDY OF THIOL ADSORPTION ON GAAS(001) SURFACE ....................................................... 162 6.1 Introduction ........................................................................................... 162 6.2 System and Simulation ......................................................................... 163 6.3 Results and Discussion ........................................................................ 166 6.3.1 GaAs 8qu Property Reproduction ............................................... 166 6.3.2 GaAs(001) Surface Reconfiguration ............................................ 167 6.3.3 Thiol adsorption on GaAs(001) Surface ...................................... 170 6.3.4 Empirical Potential Fitting ............................................................ 174 6.4 Summary .............................................................................................. 176 CHAPTER 7 NEMD STUDY OF THERMAL ENERGY TRANSPORT IN GAAS-SAM—GAAS JUNCTIONS ...................................................................... 178 7.1 Introduction ........................................................................................... 178 7.2 Results and Discussions ....................................................................... 178 7.2.1 Potential for GaAs Substrate ....................................................... 178 7.2.2 Thiol Adsorption on GaAs(001) Surface ...................................... 181 7.2.3 NEMD Study of SAM-GaAs Interfacial Thermal Transport .......... 182 7.3 Summary .............................................................................................. 199 CHAPTER 8 NEMD STUDY OF THERMAL ENERGY TRANSPORT ACROSS GRAPHENE-POLYMER INTERFACES ............................................................ 201 8.1 Introduction ........................................................................................... 201 8.2 Simulation and System ......................................................................... 202 8.3 Results and Discussions ....................................................................... 204 8.4 Summary .............................................................................................. 211 CHAPTER 9 CONCLUSIONS ........................................................................... 212 CHAPTER 10 FUTURE DIRECTIONS ............................................................. 219 10.1 QMMM in Thermal Energy Transport across Material Interfaces ........ 219 10.2 Electron-Phonon (e-ph) Interaction in Thermal Energy Transport ....... 219 10.3 Electron Contribution to Thermal Energy Transport ............................ 220 10.4 A More Sophisticated Potential Model forAu Substrate ...................... 221 10.5 A Full Anharrnonic Potential for SAM Molecules ................................. 221 BIBLIOGRAPHY ............................................................................................... 222 vii iaiele . Table . iabée . lab-29 Table Table Taele Tabie Table Table Table Table Table Table Tabie Table 4.1 Table 4.2 Table 4.3 Table 4.4 Table 4.5 Table 4.6 Table 5.1 Table 5.2 Table 6.1 Table 6.2 Table 6.3 Table 6.4 Table 6.5 Table 7.1 Table 7.2 LIST OF TABLES Potential Functions with Parameters Used in the Simulations... Finite Size Effect on Thermal Conductivities Substrate Thickness Effect on Thermal Conductivities and Thermal Conductance Simulated Pressure Dependence of Thermal Conductivities/conductance in x-, y- and z-directions .............. Conductivities in x—, Chain Length Effect on Thermal y-directions at 100K...... .. .. Molecule Coverage Dependence of Thermal Conductivities and Thermal Conductance... Thermal Conductance of Systems with the Constant TemperaturesMethod Thermal Conductance Systems with Heat Fluxes Imposed ...... Calculated bulk properties using different exchange-correlation functionals... Energetics of Reconfigured GaAs Surfaces............ Characteristic Dimensions of As-terminated GaAs(001) Bond lengths and binding energies of thiol-on-GaAs... .... Functional and parameters forAs-S Morse Potential... 8qu properties predicted using empirical potentials... LJ parameters for different atoms... 76 99 109 110 113 122 124 167 169 170 173 176 179 181 Table E Table 8.1 LJ parameters for polymer and graphene simulations... .. 203 Figure Flare Figure Figure Figure Flgre Flgure Figure FQUE Figure FIQUTE FEUm Figure Figure FQUK figure Figure 1.1 Figure 2.2 Figure 2.3 Figure 2.4 Figure 2.5 Figure 2.6 Figure 3.1 F igure.3.2 Figure 3.3 Figure 3.4 Figure 3.5 Figure 3.6 Figure 3.7 Figure 3.8 Figure 4.1 Figure 4.2 LIST OF FIGURES Images in this dissertation are presented in color Moore's Law[1] Typical energies in a CPMD run [20] A typical model used in nonequilibrium method to calculate Temperature profile in 2 direction [35] .. Linearfit ofthe nonlinear temperature profile [35] Heat correlation function decay with the time [35] Material combinationsofthin layers Periodic boundary condition in 81-69 system... Temperature evolutions in the two contacting thin layers Atomic vibration power spectra (Si-Si) .................................. Temperature evolutions in contacting thin layers (Ge-Ge)... Atomic vibration power spectra (Ge-Ge)... .. Temperature evolutions in the two contacting thin layers Atomic vibration power spectra (Si-Ge)... Simulated Au-SAM-Au systems of different sizes... .. Procedures of preparing the Au-SAM-Au simulation system 20 27 28 29 31 43 52 53 55 81 83 figure . Figure Figure Figure Flgure Figure Figure qum FQUK figure Flgllrr Figure Figure 4.3 Figure 4.4 Figure 4.5 Figure 4.6 Figure 4.7 Figure 4.8 Figure 4.9 Figure 4.10 Figure 4.11 Figure 5.1 Figure 5.2 Figure 5.3 Figure 5.4 Figure 5.5 Absorption sites and tilt directions of SAM molecules on Au(111)surface [81] A Typical Normalized Heat Current Auto—Correlation (HCAC)Function Different thermal conductivity profiles obtained from the Temperature Dependence of ln-plane Thermal Conductivities of Systems with Free Boundary Condition in z-direction Temperature Dependence of Thermal Conductivities and Conductance of Systems with P80 in z-direction... Out-of-plane Thermal Conductance vs. Simulation Cell Temperature Dependence of Au-SAM-Au Junction Thermal Conductance with Different Alkanedithiol Chain Lengths ..... Experimental Data of Temperature Dependence of Au-SAM-GaAs Junction Thermal Conductance with DifferentAlkanedithiol Chain Lengths [84]... Normalized Vibrational Power Spectra of Au Substrate, SulfurHeadsand Carbon atoms ATypical Set-up ofthe Simulation System...... Steady state temperature profiles...... Heatfluxes... Au-SAM interfacial thermal conductances at different mean temperatures... vooss of surface Au (1" column), s (2"d column), 01(3rd column), C2(4"‘ column) and C4(5"' column) atoms at 87 92 104 105 109 111 111 115 120 128 131 134 Figure 5. meS wa5 me5 EMS I meE Figure l mel HWW Figure Figure Figure Figure Figure temperatures of 350K (top row), 250K (2" row), 150K (3rd row) and 50K (bottom row)... 138 Figure 5.6 Interfacial thermal conductance versus simulation supercell 144 thicknesses... Figure 5.7 Coverage dependence of Au-SAM interfacial thermal 146 conductanceat100K Figure 5.8 Binding Energy Dependence ofAu-SAM Interfacial Thermal 149 Figure 5.9 VDOS’s of the surface Au atoms and the S atoms for the 150 original binding energy and 130% binding energy... .. Figure 5.10 Snapshots ofTemperature Profiles... 153 Figure 5.11 VDOS’s of the Au atoms at different locations of the 156 substrate Figure 5.12 Frequency dependent temperatures of the LF and IF 159 modesoftheCatoms in the SAM molecules Figure6.1 An8-layerGaAssystem 165 Figure6.2 GaAssubstratewittherminated onlhebottom 165 Figure6.3 ThiolonGaAs(001)substrate 166 Figure 6.4 As-terminated GaAs(001)surfaces 168 Figure 6.5 S—CH3 adsorption on As-terminated GaAs(001) surface 171 Figure 6.6 S—(CH2)2-CH3 adsorption on As-terminated GaAs(001) 172 surfacetypel Figure 6.7 Electron Isodensity Surface corresponding to 0.02 a.u. of a 174 S-(CH2)3—CH3 bonded on the GaAs(001) surface type I ..... Figure 1 Figure Figure Figure Fig-ere Figure Figure FIQUFE Figure Figure 6.8 Figure 7.1 Figure 7.2 Figure 7.3 Figure 7.4 Figure 7.5 Figure 7.6 Figure 8.1 Figure 8.2 Figure 8.3 Figure 9.1 S-Asbondenergy asafunctionofbond Optimized As-tarminated type I GaAs(001) suface... Thiol adsorbed on a GaAs(001) surface... System and temperature profile of NEMD on crystalline System and temperature profile of NEMD on a VDOSofatomsatdifferentz-coordinates........................ Participation ratio of (a) GaAs-SAM-GaAs junction (14 atoms), (b) SAM molecule (10 atoms) and (c) pure carbon segments (8atoms) Pure amorphous polymer system and me NEMD temperatureprofile...................................................... Single graphene sheet-polymer interface system and the steady state temperature profile Graphite-polymer interfaces and temperature profile... Regime Map for Thermal Energy Transport...... 175 180 182 184 187 193 197 205 210 216 As At AT AT OSC S ‘S N '6' q is. NOMENCLATURE element a Net kinetic energy change Velocity scaling interval Temperature difference Fluctuation of temperature Fictitious mass of orbital degrees of freedom Relative error Total wavefunction nuclear wavefunction electronic wavefunction frequency Cross sectional area perpendicular to the heat flux direction Dar (m) D (0) a,ne D (60) a ,eq amplitude factor of nuclear wavefunction Vibration density of state at frequency a) VDOS’s from non-equilibrium runs VDOS’s from equilibrium runs Absolute error convenient phase factor Total energy of atom i mean kinetic energy over all the atoms Force vector force on atom i due to its neighbor j from the pair potential Thermal conductance 5" .F f Ki 110') Planck’s constant Hamiltonian Hamiltonian of electron system heat current vector Heat flux Thermal conductivity electronic fictitious kinetic energy Car-Parrinello Lagrangian Electron mass Nuclear mass dimensionality electronic one-body density xvi N rs“ ‘ 2 part cs“ I k. .31 03‘ ‘ N h 72(0) Participation ratio Position vector of electrons Position vector of atom i Position vector of nuclei Phase of nuclear wavefunction Frequency dependent temperature Temperature Average temperature Mean temperature velocity volume xvii 1.1 110 at sorr ransis‘ phenor CHAPTER 1 INTRODUCTION 1 .1 Motivation "..The first microprocessor only had 22 hundred transistors. We are looking at something a million times that complex in the next generations—a billion transistors. \M'iat that gives us in the way of flexibility to design products is phenomenal." —Gordon E. Moore Figure 1.1 Moore's Law [1] 1970 bumnsjsnmns uxoooxxxrooo ldkyoreflslgavv /{./w uxrooqooo 4// If .I/ / .- i’/ ////,a Looqooo i‘/’ ,1 / //”‘ uxooo {/1/ I 1970 1980 1990 2000 2010 In 19 popularlyl doublesal ltlslawT iedlodevl Silie 310280m 281111 by The ultra- 7330-91191 geometric Studied e; tampon, SYSIETTIS: allenn‘on In . ”Men: SWfan ””1- This 00ndUlrlii In 1965, Intel co-founder Gordon Moore saw the future. His prediction, now popularfy known as Moore's Law, states that the number of transistors on a chip doubles about every two years. Engineers we astonishingly clever at maintaining this law. This anticipated trend towards miniaturization of electronic devices has led to device featu‘es in the nanometer range. Silicon integration suppliers have been producing 45nm CPUs since 2007 [1] and 28nm silicon-on—insulator transistors are predicted to reach a gate length of 28 nm by the year 2009 [2]. Together with the trend towards miniaturization, are the ultra-fast transportation phenomena. The time scales of interest in nano—engineering ranges from femtoseconds to nanoseconds. In such confined geometrical and time scales, charge transport and mass transport have been studied extensively for many years. On the other hand, nanoscele thermal energy transport, which has great influence on the performance of most miniaturized systems such as CPU chips and molecular electronics, has not gained as much attention as other transport phenomena. In the current electronics industry, one of the major challenges in high-density Si-based electronics is power dissipation. Thermal management is slowly emerging as a potential barrier for pushing transistor gate lengths below 50 nm. This is despite the fact that single crystal Si has a room-temperature thermal conductivity of about 140 Wlm-K, which is much higher than other relevant semiconductors (SiGe —- 10 WIm-K; GaAs - 50 WIm-K; InP — 70 WIm-K; lnGaAs — 10113 Hews fats belle COM 68p: \ “an: ”801 10 WIm-K). Besides the well developed Si-besed electronics, organic electronics [3] are steadily emerging as a promising candidate for low-cost information processing, storage, display, and communication devices. The aforementioned thermal management problem in Si—based electronics is aggavated by the fact that amorphous polymers have thermal conductivities in the range of 0.1— 0.217 /m-K , which is much lower than other most inorganic semiconductors and metals. However, several recent studies [4,5] have found that the thermal energy transport in chain-like molecules is ultra-efficient along the molecule chains. To better utilize the high thermal conductivity along the molecule chains, techniques are needed to manufacture these molecules so that they line-up in appropriate order. Fortunately, some chain-like molecules usually have the self-assembly feature which facilitates manipulation at the nanoscele. Since the materials of the popular electronics have large thermal conductivities, thermal energy transport across interfaces in these devices is especially important, and they often play the role of a barrier to thermal energy dissipation. As a result, understanding and predicting thermal transport across nanoscele material interfaces becomes essential to further the advance of nanoengineering. 12 Reseal Them obse'ied I using the a 31.955 an lansmssr wrprelier 0811311131. reflection 2 because 11 order of lf arid Pohl 910110115 ll Talsmissi Gallery 01 1.2 Research Background Thermal resistance at the interface between solid Cu and liquid He was first observed by Kapitza in 1941 [6] and subsequently explained by Khalatnikov [7] using the acoustic mismatch model (AMM) of phonon reflection and transmission across an interface. Much of the subsequent work has focused on phonon transmission across solid-solid interfaces, which was summarized in a comprehensive review by Swartz and Pohl [8] and by Cahill et. al. [9]. In particular, it was shown that the assumption of specular (or non-diffuse) phonon reflection and transmission breaks down for high-frequency phonons (>100 GHz) because the chemical and physical roughness of an interface can become on the order of the wavelength of the dominant heat-carrying phonons. This led Swartz and Pohl to propose a diffuse mismatch model [4] (DMM), which assumes that phonons lose their memory after reaching the interface, and that the probability of transmission to either side of the interface depends on the ratio of the phonon density of states. However, the DMM seems to continuously over-predict the solid-solid interfacial thermal resistance. For solid-molecular interfaces, because the phonon densities of states between the metal and the molecule are vastly different, the DMM is expected to predict a large thermal resistance. It must be noted that no one has systematically applied DMM to metal-molecule interfaces, although some predictions [10] without accounting for molecular vibration spectra fall within a factor of 2. It I can be nanesc imerad simulal atoms pcienb Simula' 0011031 particle fibula This is hand; flare L“Tiler. interla Tans; Valid It has become increasingly clear that the classical MD simulation technique can be a very valuable computation tool in studying energy transport at the nanoscele. Nanoscale particle-particle, molecule-molecule, and phonon-phonon interactions have been studied [11-16] using classical MD. The MD computation simulates real phenomena by solving the equations of motion for an ensemble of atoms that exert forces on each while in the presence of the intermolecular potential fields surrounding them. With the use of large numbers of atoms in the simulation, physical properties can be determined using ensemble averaging. As opposed to energy transport models that treat the energy transport units as particles and as a result use the Boltzmann transport equation, there is no need to calculate the relaxation time of various scattering events in the MD simulation. This is a major advantage with the molecular dynamics simulation. On the other hand, MD simulations can be combined with Boltzmann equations, since relaxation time can be extracted from the generated trajectories. It must also be pointed out that one of the biggest challenges in fundamental understanding of material interfaces is behavior of the atomic bond at the interfaces. The atomic bond at the interface also greatly influences the thennal transport efficiency across the boundaries. However, there is yet a universally valid potential to simulate the interfacial atomic bonds that can be used in classical MD simulations. As a result, a precise characteristiztion of the interface bonding and the development of a set of accurate empirical interatomic potentials 519' 0191'; re. bat its 1 for Ie-i: 1115. int 1.3 in] are critical to any interface related studies. The empirical potential functions are often determined by fitting parameters to reproduce one or several properties of a certain system from experimental measurements or ab-initio calculations. An alternative route to study the interface problem without the need of developing of an empirical potential is the ab-initio method. Ab—initio MD computes the forces acting on the nuclei from electronic structure calculations that are performed “on-the-fly" as the atomic trajectory is generated [17]. Despite its accuracy and independence of empirical potentials, ab-initio method suffers from the problem of only being able to simulate very limited numbers of atoms because the "on—the-fly' calculations are very computing resource demanding. In heat conduction of solid state materials, there exist two kinds of heat carriers: phonons and electrons. Phonons are the quanta of the lattice vibrational field [18]. Phonons dominate the energy transport in semiconductors and insulators. \M'iile in metals, electrons are important energy transporters. In this study, the phonon part of thermal energy transport across material interfaces is studied by both ab—initio and classical MD. 1.3 Specific Research Problems and Objectives There are several nanoscele interfacial thermal transport problems studied in this project. (1) ab—initio MD simulation of themial equilibration in ultra-thin layered semiconductor superlattices. Since the ab-initio MD can model any atomic interaction without the need of empirical potentials, and such a method has never been used in thermal energy transport studies, the objective of this part is to take the first step to use the empirical potential-free ab—initio MD to simulate a simple phenomenon, thermal energy equilibration, in silicon and/or germanium semiconductor superlattices with layer thickness of several angstroms. (2) Equilibrium MD (EMD) calculation of thermal conductivity and thermal conductance of Au-SAM-Au junctions. SAM-on-Au substrates is one of the most studied molecular electronics systems. A set of well established empirical potentials is available for classical MD simulations. In this part, the alkanedithiols (’5 " ((3172),. - S ") are used as the SAM molecules. The objective is to use classical EMD to calculate the thermal conductivity and thermal conductance of the Au-SAM-Au junctions utilizing the Green-Kubo method. (3) Non-equilibrium MD (NEMD) study of thermal energy transport in Au—SAM-Au junctions. Different from the EMD which relies on the fluctuation-dissipation theory to obtain thermal transport properties, NEMD mimics the experimental setup by applying a constant heat bath to the system, and uses the resultant thermal profile and heat current to calculate the thermal conductance according to the Fourier’s Law of heat conduction. The objective of this part is to use NEMD to calculate err elec lnzet ban Slml ab-i Gal W‘l‘l d8: NE Stu ll‘li thennal conductance, and to study the physical mechanism of thermal energy transport across the AuSAM interfaces by analyzing the observed phenomena. (4) DFT characterization of thief-GaAs bondings and parameterization of an empirical potential to simulate the thiol-GaAs bondings. SAM can be grown on GaAs substrate to form another kind of molecular electronic system. Unlike the Au—SAM interface, there is no well established interatomic potential for GaAs.SAM interfacial bondings. To study the thermal transport across such interfaces, high quality potentials need to be developed to simulate the SAM-GaAs interactions. The objective of this part is to use an ab—initio method which employs DFT to characterize SAM molecule absorption on GaAs [001] surfaces, and to calculate the bonding energies. The Morse potential which models the SAM-GaAs covalent bond is to be parameterized for classical MD simulation. (5) NEMD study of thermal energy transport in GaAs-SAM-GaAs junctions. Having obtained the empirical potential for the SAM-GaAs interactions, classical MD simulations can be performed. The objective of this part is to use NEMD to calculate the thermal conductance of GaAs-SAM interfaces, and to study the physics of thermal energy transport across such interfaces. (6) NEMD study of thermal energy transport across graphene-polymer interfaces. As graphene-polymer nanocomposite materials become increasingly popular, 10 the need of understanding the thermal energy transport across graphene-polymer interface is urgent. The objective of this part is to provide some preliminary results from NEMD simulations of thermal energy transport across graphene-polymer interfaces. 11 ZIAl inter densl‘ tom 1 2.1.1 21.1 Wail Tlllcl deb Rn ham CHAPTER 2 THEORIES AND METHODOLOGIES 2.1 Ab-initio Molecular Dynamics In this section, brief descriptions of ab-initio theory, approximation, and important equations are presented. The most widely used ab-initio method - density function theory is also introduced. The contents are abstracted mainly from ref. [17] which has a more detailed description. 2.1.1 Born-Oppenheimer Molecular Dynamics 2.1.1.1 Theory As we know, any material can be described to consist of nuclei and electrons. Because electrons travel thousands of times faster than nuclei, electrons can follow any nuclear motion with almost infinitesimal time elapse. As a result, the electronic and nuclear motion can be decoupled since the electrons can respond instantaneously to any change in nuclear coordinates. The many-body wavefunction of the electron-nucleus system can be written as the product of a nuclear and an electronic wavefunction. Since the electronic wavefunction only depends on the instantaneous nuclear configuration, and not on time, its behavior can be described with the time-independent SchrOdinger equation. On the other hand, the nuclei are massive enough to be treated as classical particles, 12 responding to the electronic forces according to Newton's second law. Such a treatment, which decouples the electronic and nuclear degrees of freedom, is called the Born-Oppenheimer approximation. The quantum molecular dynamics method that utilizes the Born-Oppenheimer approximation is called Bom—Oppenheimer molecular dynamics (BOMD). In a BOMD, the time-independent electronic structure is calculated in each MD step as a function of the nuclear coordinates. In the following sections, important equations of BOMD are derived briefly based on the content of reference [17]. 2.1.1.2 Derivation of Equations The time-dependent SchrOdinger equation governing an electron-nucleus system is . a -. ~ -' -* lh-a7q)({n’i}’{r"’l};t) = Hq)({re,i}’{rn,l};t) (2.1) —o re,1 a r ",1 are the electronic degrees of freedom and nuclear degrees of freedom; (I) is the total wavefunction of nuclear-electronic system; and H is the system’s Hamiltonian which can be expressed as h2 h2 H=— V2— V.2+V F. , -' :ZMI I 22m 1 n—e({ e,1} {rn,]}) e ii2 2 _. . =—Z 2M Vi +H.({r.,,1,tr.,,}) (22’ I I 13 were iiuelea eibsys (I) where V"-.. is sum of the electron-electron, electron-nuclear, and nuclear-nuclear Coulomb interactions, He is the Hamiltonian of the electron m subsystem, M I is the nuclear mass and e is the electron mass. Separate the contributions from nuclei and electrons to the wavefunction ({i,,-}.{?,,,,-};t); <1>({F.,.},{r;,,.};t) s 1140;...-};t);t({ii.,.-};t)epr-Ii2 Idt'E.(t')l (2.3) where Ee is a convenient phase factor; V1 is the electronic wavefunction and Z is the nuclear wavefunction. Inserting equation (2.3) into equation (2.1) and (2.2) yields .§_I,I_/___ it that— :2 2 Vizl/l + { IdfiZ*K-eZ}I/’ (2.4) e i §£=—Z hz V?;(+{JdFI//*H I//}Z 0t ,- 2m ' e e Eq. (2.4) describes the electronic degrees of freedom and eq. (2.5) describes the (2.5) e nuclear degrees offreedom. From these two equations, we can see that given the nuclear structue, the electronic degrees of freedom can be calculated, and given 14 the electronic structure, the nuclear degrees of freedom can be calculated. Rewrite the nuclear wavefunction as 2015,33) = B({i';,.-};t)eXP[iS({?;,r};t)] (2.6) where B is an amplitude factor and S the phase. Substitute (2.6) into equation (2.5) and take the classical extreme in which it —> O , then we get 63+2: 1 07,592+ jdéw‘H.r/=o 317 , 2M, (2.7) We write this equation as a Hamilton-Jacobi formulation as .. E+H({rn,1}a{VIS})=O (23) with the classical Hamilton function H({?;,,},{E})=T({E})+V({f;,r}) (2.9) It appears that 131 5 V15 (2.10) O a... and from Newtonian equation of motion P1 = “V1 V( {’11 I ) we find B = MIrnJ (t) = TV: Idf'eW He‘l/ (2.11) According to the Born-Oppenheimer approximation, which states that electronic si’Stem is at the ground—state at each instance of fixed nuclear positions, equation (2.11 ) Wes 15 Equaic eeCUOi define haslo fispl 21.21 ZIJL‘ diflan 30.11 Howe Sebiy lien CWhp Was 1 The 5 much He W0 >} (2.12) Equation (2.12) together with time-independent Schrddinger equation for an Mr... (0 = —V, 13mm electronic system EoI/jo = HeWo (2.13) define the BOMD. In each molecular dynamics step, the minimum of (H e > has to be reached. Based on the BOMD, the Car-Parrinello MD can be derived. This procedure is carried on in the next section. 2.1.2 Car-Parrinello Molecular Dynamics (CPMD) 2.1.2.1 Motivation OF CPMD ln BOMD, as we discussed in the last section, there are no electron dynamics involved in solving the Bani-Oppenheimer equations (2.12) and (2.13), so, these equations can be integrated on the time scale given by nuclear motion. However, this means that the electronic structure problem has to be solved self-consistently at each molecular dynamics step which introduces a relatively high computational workload. A non-obvious approach to cut down the computational expenses of DFT-MD, which includes the electrons in a single state, was proposed by Car and Parrinello in 1985 [19]. The CPMD takes advantage of the smooth time-evolution of the dynamically evolving electronic subsystem as much as possible while it makes an acceptable compromise concerning the 16 Teri-b 0‘ election evolve ' art‘s av 11.2.2 Tl quanlu nucleat energy lliis gr Woes 911990 ”Iran the lo FESpe Tome r6396 length of the time step which is in the nuclear motion time scale. In the CPMD, the electronic wavefunctions are calculated as time-dependent quantities which evolve with time elapse. However, the evolving electronic wavefunction never drifts away from the ground-state. 2.1.2.2 Car-Parrinello Lagrangian and Equation of Motion The basic idea of the Car-Parrinello approach can be viewed to exploit the quantum-mechanical adiabatic time-scale separation of fast electronic and slow nuclear motion by transforming that into classical—mechanical adiabatic energy-scale separation in the framework of dynamics theory. In order to achieve this goal the two-component quantum/classical problem is mapped onto a two-component purely classical problem with two separate energy scales at the expense of loosing the explicit time-dependence of the quantum subsystem dynamics [17]. If a suitable Lagrangian can be constructed for the nuclear-electronic system, the force on the nuclei can be calculated from the derivative of a Lagrangian with respect to the nuclear positions. Similar to the nuclear degrees of freedom, the force on the electronic orbitals can be calculated using a functional derivative with respect to the orbitals. The Lagrangian Car and Parrinello postulated was [19] 17 "where canb A5181 .ll 1 -: 1 . . L =Z'2'Mrr:,1+Z§#r-<¢rI¢r-> 5? I J —+constraints (2,14) C___V__J potenfialcnetgy where M- are the “fictitious masses” of orbital degrees of freedom. The Lagrange equations of motion for both nuclear positions and the orbitals L can be derived from eq. (2.14): d 6ch _ 6ch git-651;“, _ 05.; (2°15) d 5ch _ 5ch — (2.16) dt 6r): 5(1)? After some manipulation, the Newtonian equations of motion can be obtained 1: a 0 . M Ir", 1 (t) = - 6? ( He V10) + 6" {constramts} (217) 72,1 n,1 He 6 . I/l >+ ,, {constraints} 0 M (2.18) I .. 5 fli¢i (t) z: — 5¢* (W0 1 Similar to the thermodynamics temperature which is proportional to the 18 "‘13 as: rrlr my 11th Thu: 42 nuclear kinetic energy ( Z M I r 11.1 ), a ‘fictitious temperature' that is related to the 1 electronic orbital degrees of freedom can be defined to be proportional to the orbital ‘klnetic energy’ (Z-ui <¢i I Q) ). According to the Bom-Oppenheimer l approximation, the electrons are always in the ground-state which has the minimum energy. As a result, the ‘fictitious temperature’ associated with the electronic subsystem is always kept around the minimum energy m1n{¢i}(i//0 He [20] for an example). From the figure, we can see that the electronic fictitious W0) , which is called ‘cold electrons’ (see Figure 2.2 kinetic energy, Kr , is kept at a value far less than the nuclei Hamiltonian H] . Thus, a ground-state wavefunction optimized for the initial configuration of the nuclei will stay close to its ground state during time evolution if it is kept at a sufficiently low temperature. 19 tit—Lirgzy (hurt 1‘00) [5 Figure Electro Fr -7. 16 " ' H ’5 -7. 18 - <1) ._ H 7. 16 i—> , H1 is ..C '7.18 '- V —7. 16 - :>. on EKS E —7. 18 LIE—1 311035 . l ' " ' i Ki ' I . . ' . I l ' ' . I i ' ' ' I 0 7 140 147 252 259 time (10‘3 atu) Figure 2.2 Typical energies in a CPMD run [20] ( K, is fictitious kinetic energy of memomcmmnmsmmembm,lfls eummommmwwemyrunmmowpnen For practical calculations which simulate long physical times, the temperatures of the cold electronic subsystems are often maintained at a very low value compared to the thennodynamice temperature using Nose-Hoover thermostats [21,22]. In such a way, the ‘cold electrons’ will stay on the Bom-Oppenheimer surface, and will not exchange energy with the high temperature nuclei. This is possible if the power spectra of the electronic and the nuclear dynamics do not have substantial overlap in the frequency domain so that energy transfer from the “hot nuclei” to the ‘cold electrons' becomes practically impossible on the relevant time scales. 20 2.1.3l cti-tar." V 1 sysle! 00131: when WhICI Sails slmp Wave 011m 353C 2.1.3 Density Functional Theory Now that equation (2.12) for BOMD and equation (2.17) for CPMD are obtained, the next problem is how to calculate the forces on the nuclei, He V1 “3’“!ch W0 >} . The total ground-state energy of the interacting system of electrons with classical nuclei fixed at positions {rnJ} can be obtained through the minimum of the Kohn-Sham energy He . __ . KS IBEIKV’O Wo>}-‘gf?E [{8}] (2.19) where the Kohn—Sham energy E“ is represented as [23] EKS[{¢.}]=T.[{¢.~}]+ jar Vm(?)n(F) +£- [df V,,(r’)n(i~‘)+E,,[n] (2.20) which is an explicit functional of the set of auxiliary functions {(0.17)} that satisfy the orthononnality relation (r). 1%) =51,- . This is a dramatic simplification since the minimization with respect to all possible many-body wavefunctions {W} is replaced by a minimization with respect to a set of orthonormal one-particle functions, the Kohn-Sham orbitals {(01} . The associated electronic one-body density or charge density 21 ls obi Hare elem exd‘i lsal ace 2 210') = 2f. | ¢.(r)l (2.21) is obtained from a single Slater determinant built from the occupied orbitals, where {it} are integer occupation numbers. The first term in the KohnSham functional equation (2.20) is the kinetic energy of a non-interacting reference system consisting of the same number of electrons exposed to the same external potential as in the fully interacting system. The second term comes from the fixed external potential. The third term is the Hatree energy, is. the classical electrostatic energy of two charge clouds which stem from the electronic density. The last term is the contribution from the exchange-correlation functional, and this term is the most intricate one. This term is a lumped term of exchange and correlation effects. (More detailed descriptions of these terms can be found in reference [17]). There are many schemes to approximate the exchange-correlation function. The most popular two are the Local Density Approximation (LDA) [23] and Generalized Gradient Approximation (GGA) [24]. 2.1.4 Plane Wave Basis In many current density functional theory codes, plane waves are used as the basissetinordertorepresentthe periodicpartoftheorbitals (A. The method makes use of the reality that the ubiquitous periodicity of the lattice produces a periodic potential and thus imposes the same periodicity on the density. Plane waves are originless functions, i.e. they do not depend on the positions of the nuclei. This facilitates force calculations. It also implies that plane waves are an unbiased basis set. By increasing the energy cutoff Ea“ the quality of the basis can be improved. Another feature of the plane wave basis is that derivatives in real-space are simply multiplications in G-space, and both spaces can be efficiently connected via Fast Fourier Transforms (Ff-T). Thus, one can easily evaluate operators in the space in which they are diagonal. 2.1.5 Pseudopotentials In order to minimize the size of the plane wave basis necessary for the calculation, core electrons are replaced by pseudopotentials. Pseudopotentials are required to conectly represent the long range interactions of the core and to produce pseudo-wavefunction solutions that approach the full wavefunction outside a core radius rc . Inside this radius the pseudopotential and the wavefunction should be as smooth as possible, in order to allow for a small plane wave cutoff. In addition, it is desired that a pseudopotential is transferable, this means that the same pseudopotential can be used in calculations for different system under different conditions with comparable accuracy. 23 2.1.6 Limits of w-initio MD Although ab-initio MD is accurate because of the 'on-the-fiy' feature, and it can theoretically be applied to any system without suffering the transferability problem, it requires a large amount of computational resources and can only model a very limited number of atoms (up to several hmdreds) for very limited physical timescale (up to tens of picoseconds). Such a limitation makes the classical MD attractive in many research areas where large systems and long physical times need to be simulated. 2.2 Classical Molecular Dynamics Classical MD simulation is the modern realization of an essentially old-fashioned idea in science; namely, the behavior of a system can be computed if we have a set of initial conditions plus forces of interaction [25]. This classical simulation method has long been used and is well developed as a tool in statistical mechanics and chemistry [26]. In principal, the MD method can be applied to all phases of gas, liquid and solid and to interfaces of these three phases. MD is a computational method that simulates the real motion of every single atom that builds up the system of interest by employing empirical potential fields, and then solving Newton’s equation of motion. In fact, Newton’s equation of motion can be derived directly from the ab—initio equation of motion discussed in section 2.1. 24 Fit nuclei éniegra celenil claret Teailsli lire sir Slmula beset Cale-u. I0 lie. the C, Is brr and E From equation (2.11), the nuclei equation of motion can be expressed as "V1 V1921» = M1521 (t) = F (2.22) According to the Bom—Oppenheimer approximation, V is only a function of nuclei coordinates. In classical MD, the electronic degrees of freedom are integrated out beforehand, and empirical functions are employed to express the potential V . Depending on the nature of electronic structure or simply bonding characteristics, the choice of suitable potential functions is extremely important for realistic simulations. The accuracy of the potential function limits the accuracy of the simulation, and thus the accuracy of any information to be extracted from the simulation. Fortunately well-tested empirical potentials are available for most of the materials of interest in this study. These “predefined“ potentials are either based on experimental data or on independent dz-initio electronic stmcture calculations. These predefined potentials are the solutions to the question of how to describe or in practice how to approximate the interatomic interactions which is the core problem of any molecular dynamics scheme. Typically, the full interaction is broken up into two-body, three—body and many-body contributions, long-range and short-range terms, which have to be represented in suitable functional forms. However, the need to devise a “fixed model potential” implies serious drawbacks. First, the predefined potentials suffer from the transforrnability 25 gamer (roped hector diener ebebo arise lntrac be a transit atren lager tab», 510018: 13111 abdrn direct] °Ilhe phOno Ills Se main], problem. The potential functions are always fitted to mproduce one or several properties of a certain system under certain conditions, but when these fitted functions are used to calculate other properties, or the simulation is carried on different conditions, the result becomes questionable [27,28]. Secondly, if the electronic stmcture, and thus the bonding pattern, change qualitatively in the course of the simulation, a ‘fixed model potential” can not represent the interactions for all time while the simulation proceeds. These disadvantages affect the accuracy of simulations involving chemical interactions or phase transformation. Fortunately, in this research, our simulations do not involve the aforementioned phenomena. Moreover, classical MD is a good tool to simulate large scale systems for a relatively long physical duration (up to nanosecond scale), under which conditions the thermal energy transport phenomena can be studied. 2.3 Molecular Dynamics Simulation in Thermal Energy Transport Classical MD has become more and more popular in studies of nano-scale and micro-scale thennal transport. Research using classical MD is mainly in two directions: one is using the MD trajectory to investigate the physical mechanisms of thermal transport. Studies in this category often focus on characteristics of phonons, the main heat carriers in semiconductors, and their interactions [29-31]; the second direction is predicting the thermal transport properties. Such works mainly focus on the thermal conductivity calculation [32-34]. 26 There are generally two kinds of molecular dynamics simulations in calculating thermal conductivity. One is NEMD, which is also called the direct method. This method is analogous to the experimental measurement of thennal conductivity. One approach is to impose a steady heat flux across the direction of interest of the modeling system. A typical example of this method is shown in Figure 2.3 [35]. Z=-Lz/ 2 Z=-Lz/ 4 2:0 Z=Lz/4 Z=Lz/2 I l I I I J2 J2 Jz ‘— ‘I'AE ———* -A8 a— s—b <—-——> 5 / 5 \ l 2 Figure 2.3 A typical model used in nonequilibrium method to calculate thermal conductivity [35]. To impose a steady heat bath, thermal energy A3 is added or subtracted inathinslabofthickness 5 centeredat z=1th/4 ateachMDtimestep. Each particle velocity in the source/sink region is scaled by the same factor such that the resulting net kinetic energy is increased or decreased by an amount A5 . Equilibrium between this added kinetic energy and the system’s potential energy is expected to be reached within a typical vibrational period (<<1ps). When the system achieves steady state, the heat flux is given by J2 = A5 / 2AAt . 27 The 1er picos MC 119116 caTrier The resulting temperature gradient is then calculated from the data of the atomic velocities (see Figure 2.4 [35]). Usually a long simulation time, about tens of picoseconds, is needed to achieve the stationary temperature profile by using this direct method. 4x411288 Si T=500K Temperature (K) 8 450 T I I T I 1 T 0 20 40 60 80 100 120 140 z(nm) Figure 2.4 Temperature profile in Z direction [35] It can be seen in Figure 2.4 that there is nonlinear temperature profile near the heat source and sink. To cope with this nonlinearity, a linear fit is always carried out, see Figure 2.5 [35]. 28 Tiler is us liar) sud me; § e i: 0' ..5 1° 12.4K 40nm Temperature (K) & 8 4704 450 0 1'0 2'0 3'0 40 5'0 61) (a) 2 (nm) Figure 2.5 Linear fit of the nonlinear temperature profile [35] Then, Fourier’s law of heat conduction, AT J2 = ”k E (2.23) is used to calculate the effective thermal conductivity k. The nonequilibrium simulation is also an effective tool to study the heat transfer across interfaces [29,36]. However, this method has some drawbacks such as the requirement of a large number of atoms, limitation of the phonon mean free path and the need to impose unphysically large temperature gradients. Another method to calculate thermal conductivity employing MD is the so-called equilibrium approach. This approach calculates the thermal conductivity using the Green-Kubo method. This method does not have those disadvantages of nonequilibrium method. Instead of imposing a temperature difference or a heat 29 flux, the Green-Kubo method relies on small statistical temperature fluctuations to drive instantaneous heat fluxes. The thermal conductivity of the system to the disturbance can be calculated from the instantaneous heat flux autocorrelation function according to the Green-Kubo formula [37] (2.24) : ndBTm2 where n is the dimensionality; k3 is the Boltzmann’s constant; V is the o:i‘cit < .7(0)j(t) > system’s volume; T M is the system’s average temperature; (2.25) is the heat current autocorrelation function; and the heat current vector is expressed as .. d _, J :71; .- E3; (2.26) where E, is the total energy of atom i which consists of kinetic and potential —-> energies; ’1' is the position vector of atom i. For pair-wise potentials, equation (2.26) can be rewritten as J=ZEr +— 271E?” (2.27) 2i ,j, j¢i -. where F, is the force on atom i due to its neighbor j from the pair potential. 30 A problem in the Gmen-Kubo method is often caused by the very long simulation time needed to let the heat current correlation function decay to zero (see Figure 2.6 [351) 0.2 6X6X6 SI T=1 000K 0.15 -0.05 r 1 0 50 100 1 50 200 r (95) Figure 2.6 Heat conelation function decay with the time [35]. The usual method to deal with the slow decay problem is to fit the heat current autoconelation function to an exponential function and integrate it Several studies show that such an exponential fit gives good results [32,38]. However, such an exponential fit is not always suitable. McGaughey et. al. [39] found the fitting is not suitable for complex silica structures, and used a running average to smooth the noisy integration profile to define a converged region to attain the thermal conductivity values. 31 AB-II 3.1 Be Amer ultra: pier/ll molic “pier flel pro; We the pie. are Illr 01.1 CHAPTER 3 AB-INITIO MOLECULAR DYNAMICS SIMULATION OF THERMAL ENERGY TRANSPORT ACROSS MATERIAL INTERFACES 3.1 Background and Objective There are many different approaches to study nano-scale thermal transport Among them, classical molecular dynamics (MD) has the potential to simulate ultra-small and ultra-fast transport phenomena [32,33]. As discussed in the previous chapter, classical MD is a computational method that simulates the real motion of every single atom involved in thermal energy transport by using a “predefined' empirical potential and then solving Newton’s equation of motion. After a trajectory of the atoms’ motion is obtained through the simulation, properties of the system can be calculated. Thus, if one can find an empirical potential field model accurate enough to describe the interaction of the atoms in the system of interest, any property related to the motion of the atoms can be predicted through a classical MD simulation. Such empirical potential functions are often determined by fitting parameters to reproduce one or several properties of a certain system under prescribed conditions. When these fitted potential functions are used to calculate other properties or when the simulation is carried out on different conditions, the results become questionable. Not only does 32 classical MD suffer from this transferability problem, it also fails to take electrons into consideration because the “predefined' potential has already integrated out the electronic degrees of freedom. This feature prevents its application on thermal energy or electronic transport simulation in metallic materials where the electrons play an important role in transport phenomena. The reign of traditional molecular dynamics and electronic structure methods was greatly extended by the family of techniques that is called "ab-initio molecular dynamics” [17]. The basic idea undertying this ab-initio molecular dynamics methodistocomputetheforcesactingonflrenudeifiomelecfionicshucture calculations that are performed “on-the-fly" as the molecular dynamics trajectory is generated [25]. In this method, both electronic and nuclear dogmas of freedom are considered, and any system consisting of any atom can be simulated without suffering from the transferability problem. There has not been any work using ab-initio molecular dynamics to model thermal energy transport directly. In this pat of the project, an ab-initio MD approach employing Density Functional Theory (DFT) is used to simulate the thermal energy equilibration in different types of multi-thin—Iayers in different temperatures. This should be considered the first step to use ab—initio method to model lattice thermal energy transport (phonon transport). This work did not involve explicit electron transport properties which can be calculated using ab-initio method separately by calculate the electronic structures of a system. The 33 modeling of phonon transport, electron transport and even phonon-electron interaction are expected to be performed in future works within the reign of ab—initio method. The results show reasonable physics and thus suggest the possibility that the ab-initio molecular dynamics simulation is appropriate for the simulation of nanoscele thermal energy transport. 3.2 System and Simulation In this part of the project, the CPMD package version 3.11[40] was used. This CPMD program is a plane wave/pseudopotential implementation of Density Functional Theory, particularfy designed for ab-initio molecular dynamics. Visual Molecular Dynamics (VMD) [41] is used to visualize the structure of systems. In the thermal energy transport simulation, different systems were studied. The simulated systems contain 64 atoms. In terms of material, Silicon—Silicon (Si-Si) with a lattice constant of 5.431 A , Germanium—Gennanium (Ge-Ge) with a lattice constant of 5.658 A and Silicon-Germanium (Si-Ge) thin layer combinations were studied. For the Si-Ge system, the simulation cell size was chosen as twice the Germanium lattice constant (i.e. 11.316A ), in which case the silicon density is slightly smaller than its bulk density, and the germanium density is the same as its bulk density. A geometry optimization with periodic boundary condition showed that such a structure is stable. The reasons we chose Silicon and Germanium for simulation are: 1) they have similar crystal structure and 2) they are semiconductors which have well defined band gaps. Thus we do not need to consider excited states which result in much more computational time. Examples of the thin layer structures are shown in Figure 3.1. 35 Figure 3.1 Material combinations of thin layers (a.Si-Si; b.Si-Ge; c.Ge-Ge) 36 (a) 37 Figure 3.1 cont. (b) 38 Figure 3.1 cont. 39 Because both bulk Silicon and Germanium have Face-Centered-Cubic (FCC) crystal structures, the three different combinations of materials were set as FCC lattices. Geometric optimizations were perfomled and showed that FCC crystal arrangements are stable in these three combinations. Periodic boundary conditions were employed in both systems and in all three (x-, y—, z-) directions as seen in Figure 3.2. The system’s periodicity in the z-direction suggests that each layer would have thermal communication with the neighboring two layers. ‘ vs ¢ 1 o ‘ v‘ a ‘ '0' r,- ‘ r I s V 1,0 a s Figure 3.2 Periodic boundary condition in Si—Ge system The pseudopotential used for Silicon was the Trouiller—Martins norm-conserving pseudo-potential [42]. A nonn-conserving pseudopotential means that outside the core, the real and pseudo wavefunctions generate the same charge density. For Germanium, the Stumpf-Gonze—Schaeffler pseudopotential, which is also noun-conserving, was used. These pseudopotential are well establist and give enoug'i accuacy while not very calculation demanding. To deal with the exchange-correlation function in the Kohn—Sham energy (see equation (2.20)), the local density approximation (LDA) [43] was employed for both materials. The procedure of the simulation is as follows. First, the neighboring two layers are set to different temperatures using Nose-Hoover thermostats [21,22]. Then, all the thermostats are removed to let the thermal energy equilibrate. Finally, the temperature evolutions of the two layers are visualized, and the vibration power spectra are calculated. One thing worth mentioning is that the temperature here is defined according to the kinetic theory: _ 2r. T 3,93 (3.1) where Ex is the mean kinetic energy over all the atoms, and thus temperature T is used as a label indicating how much kinetic energy does the system has. 3.3 Results and Discussion 3.3.1 64-atom, Silicon-Silicon For this system, the contacting Silicon thin layers, as shown in Figure 3.1 (a), 41 are set to 100K and 50K in case (a), 300K and 100K in case (b), 400K and 100K in case (c), and the temperatures of the two layers are maintained by thermostats for 10000, 25000 and 30000 time steps, respectively. Then, thermostats are released to allow the energy to equilibrate, see Figure 3.3. Then the vibration power spectra of the two layers are calculated as seen in Figure 3.4. For this system which contains 64 Silicon atoms, there is no material mismatch. The thermal energy transport is as would be expected in the bulk material. The equilibration time required was very short (around 2000 steps 200 ff ). The final temperatures of each layer were equal to the average set temperatures. From the vibration coupling point of view, see Figure 3.4, the vibration spectra of the two thin layers overlapped almost perfectly, which resulted in fast thermal communication and thus short equilibration time. In Silicon, phonons are the thennal carriers. The thicknesses of the layers are around 0.5nm which is much less than the phonon mean free path. Phonons travel ballistically from one layer to another. We see that the time required to transport different amounts of thermal energy are almost the same. This verifies the ballistic transport of the phonons. 42 Figure 3.3 Temperature evolutions in the two contacting thin layers (Si-Si). Dotted lines refers to high temperature Si and solid lines represent low temperature Si 43 Te m pe ratureg 8 y; (a) 100K vs 50K .1 ...-nine:- a.u.-... ..r... ...... . OI.‘ .- IHUH-nuo-o~' - Dense-”...? - High Temp. Si Low Temp. Si ......a...:..:. l 1441 l .0 3.3.3-...... ......... o .... ... a.u..c ...-‘9 es ...... "...“; a o c ....... '- ...: ...... . l . I 2. se...........:............ .. ... . ... ....\ u-e .. .... Sa- ........ 2...:- o 3... so... o‘fiang-uaeu-uf: .a ..... :- u$ ..... ......o:=. gen-office’s ........ ...:...:.........::e.::....e .... ...-.....zssssss. b-——-.::~.-e ...-no... . - a.u....Juuun—zouoou ...: soooeoe..........::.......1 ..........uz§e..... ocean . ... .sss-‘ss::~u.....e..... .......... .- ‘ifiueasaesee...............o....eo . a so . a .c . J...§..:...:.:. rises..." ” 3...... .. x Ioo'OoaeE-leee::o::e . so none 0&5}...-.Huuuunnuuunuqszssz: ...::..:x.....: . .. 1:31.?- ausfi cu: . 5.33:3"Hints? .....s‘s s---s---~ .............. s ............. ..ee.:-=.a:se s ...-......oee:...-u“unflnnnuorawovdfof'Joooo- a a.-tee-us.-en...u.....u...n¢......-.-.- u.an...a...sensuous-.....oeoneae . 0 (areas-...... ..........:.n...-e..uuunnunnueunnnu—n—h: ....- --------------- unh-\~\\§“ cos—rareeuafe::::...-u} ......... e:- l .oaaoasuuvonnnvuuuuufi....-o; ...-......e..=zo alcove-....u-ue. ...-... coaavafcooeoeo......o—n:...::-eeon-eunuaonuheouoounnu-oe-sssse :- i)- oooOoosssawppnwunu .... ...-2:... ...... ...-J... ....... ............:..e.o.u :59: P b p _ . - ...:zva...............~.r:: fill r- r Bun—Id Eil- I m m w 1 ”.382.th JJ 15000 200m 10000 Time step (4 a.u. per step) ”wk-JHWLQQEQF Figure 3.3 cont. (1:) 300K vs 100K one... :- .... a e - Pia-...... .- a . .- .l l. I n- o- . r I I ll’ah 55" u. a ..a o o e :3. and! act- -.-DID.I.I'I. l ...-....a ...-Jo; In. I. 43.: :1. ... 2.53s .1141 IIIIIIIIIIII High Temp. Si Low Temp. Si v «al.-nut. i“ . ...L.......o.. . 033*7n. .3: .o issues"? . . .......:..=.... m.................. ......::...:......fiwu...::.u.u.n. ..."...feu... so... 2.33:: l issues-$5}............................-.......-.a.u.-.... e . e . es... ....... ... a . . . u . . . . u u e e . o.— o J O I ‘ ‘ ............ .x sanss........... ...:3..hhmnwflwwmmm.””5535....r l 5:25.22. 333......... l 11.5.... . .... 7...: . .................._......sis.........rit. .- ..... cannon-”.033... ....a-u'o-aa t . : . . . . . ......wmmmgzzzzzet:.... .- .. ...-......uoeue-QQQQQKT rib...“g»............................................. .. -. . nococoa:-o-ooooooosaeoaaa: 559$“? . . e. ... .. . .. .Qitg\§glonan‘o—o—W}=Q Is.‘ ‘3: a _ a....o..o..-...u..:.o: on: essence-e“: $5.53: . . .s:use:sssss~=-§=s=-~uu . . no 00:0coco-aucuunreou..en..coa “....-r n ..c‘ QIOOQCO-oonoovn I 5.5:... . $5.55.......:so..:$:....:! ts..==:~= chunks. Revert... ... 3.2.... ...: .....z::“uumwflmfiaflnnnuoooox ....u...........u.£.c§.a ...... . use.es...sesssss--~s~nununmmuwWwwunUw .ooeoeoooauuuuuuv v we . sEieu...-I-.WW o ..usaflvsancaaoooooooso- . : ............................:.............................Hfl.1..:..::.. is... . Sci . . ...-...:.....sto...-::&drkovu % u ............3............. 1.. ...-age.-:..:::.:. ...»...efiems ooeoeoeen: a}: E'- are? .0000 eat-55.3.»; ......1................. .......... .............z:=.:.==.:..=..==:.... 6‘ up . wwww"«_.._.............. inn""woe-...ceoeo-oo:-:.o......eu . {...-.2 Ca 1.3.3.3....329‘33 on! is; ...... . or... -- ll... of»: Ooooooovfeeeo.ulcuoucase's-... faafl' 0300'..- ............................................as........ - 1.333353%“: ........... . tuna". z .fifi... . - ..‘bnskufishfis ....u e a e e . e . ......o...-~.....e.ea....o_ ago... ....5...Ieeoo.-.oooeue.-- sees-...... :- o. . as D '1 I l 2‘ E Lil-tr Ck... . ......rJELMII m map—«mahwnmcwh 150 100 30000 Time step (4 a.u. per step) 10000 45 Figure 3.3 cont. (c) 400K vs 100K I .. :5... $2.... . ...... ...-l: I.I.III I..|I.DE- High Temp. Si Low Temp. 81 ~ ...... oo. .33... .mmnnmfi . so. 035...”. .355: 3 or... ...... ....... .. ......e.....§u... aaseeeeeeeeuuuunuuuunimr:sgphu-nm~mh-hfl£eesees. .. : =5: 55%.... . .... .2 ..qwmwxmgga... ........__.::......tst. I33... fie.....==ss.:.. . ......z.3:333.2.23.533333 ...... I}; ...Sw . 22223.3: ii...émo:e #3:. . so... ...... ...suuuanunnuuntetsst teeter! 32...... ..Beflgkuuiiol ~nnnao... .eeeeegmfiwwa ”mate”? a colon ...'-.:.ao-. - so"; utensil-i .......:-:=. ....... .5 {...-i ......E.’ : uvvvoaooo‘ .... noses-:-.. ......uaeueufaOa‘I . ......e...... v a"? ........ it. as. ..e......................._.. la ...... 33...: .. a.u.-NM”; lI . . . ..... . #qsreee: .uu"..’ .. l: . ._ x. E: .55... ...... o a. c ‘- rrbtch finafitirzzr 3 gig ......nflhflflflflflflhaawvioooolxn. use: sssss tees: ==2nnu ........... 3.00000: ovoauavsresv. . ... . . . e. .u ...:3 .- 8.3.3: ..."..mmmaaw.~ ..h..............u..::.::........... 9.2.5.... ._LbPr_ ul_L.. me we a m m 3.322.th [IL 10000 20000 00 Time step (4 a.u. per step) This layer Another phenomenon worth discussion is that of the temperature oscillation. The magnitude of the oscillation should be inversely proportional to the square root of the number of atoms included in the calculation of the mean temperature of the ensemble. This can be expressed as AT = ixT osc N mean (3.2) This expression explains why the temperature oscillations in the high temperature layers are larger than in the low temperature layers. 47 Figure 3.4 Atomic vibration power spectra (Si-Si) 48 Nomallzed power spectrum (a) 100K vs 50K 0.4 b _ ............................ High Temp, Si _ Low Temp. Si 0.3 '- L 0.2 - L i- 0.1 - 0 ‘ l l I l l J I I l l l l J 0 2 4 6 8 10 Frequency (arbitrary uni) 49 Figure 3.4 cont. (b) 300K vs 100K 0.4 .0 on Normallzed power spectrum 0 o '_. io ............................ High Temp. Si Low Temp. Si l l l 1 1 l I 4 6 8 10 Frequency (arbitrary unit) 50 Figure 3.4 cont. (c) 400K vs 100K 9 bl ITITIlTjIIIllj'lel‘lTllllllr .0 M Normalized pgver spectrum 0 '_. -.. m ............................ High Temp. Si Low Temp. Si 6 Freqency (arbitrary unit) 51 3.3.2 Gil-atom, Gennanium—Germanium This system consists of pure Germanium atoms as seen in Figure 3.1(c). Two contacting layers are set to 300K and 100K separately at the beginning, and the temperature difference is maintained by thermostats for 16000 steps. Then the thermostats are released to allow energy transport (Figu'e 3.5). The vibration power spectra plot is presented in Figure 3.6. r i i i- : g 400 - .- L. . 2’1 . - :3 1% :, ............ High Temp.Ge i‘ :3 : 5s =: s ': : ~?- =5 i5 F s s is; a . : s; -2 -— Low Temp.Ge 350 i. = -E “55 3 = i .5: 5':= - 5 E: i; 2252 s E .3 as : ‘ i , 3. (525:3: : : 35 :31: : 5 '- 3 2 2 z 2;: 2 ’12: :' f: 3: ::': : ;: '. - :; :3:: ::- - .- :- -:. . -. -: : = :- :-:.-. :5 5 :: :5 . :::§ 5 E: y-E :5 55 :53 f‘S: : E: .=: -'. 351: z x: 3m .;§:~ :: 552:5: 5': “a: 5.2.95 : 2r is $553 .= 5555 3‘55 555?: ’gzisé sass £29 5 i=5 5;: fish: ‘ W -.--.i 5:: 5:“; :5 g =_: z: :2: i s ' ::s i : "~ ‘: ' : :: o -: :- :' - ::; -. z 0 3E isEE-sss fésf SE := §=.=§ E 5 ;‘ = 2.5 35 : 355 .5! E_= ;-'. 3 =53. id '§ §§§ 2%! E - :'E:: .5; -. s .. : z: .‘__ 2'. (B .= '-..= i is; s - 32‘ ' . h : 3 E ‘ ;' E : :z; z: i :3 ' ' '-: £5 '2' 0 : $;-}. P f 1 s E: L I l I , l l eée 0 _ I I L I l I I I l l I I I I I I 0 10000 20000 30000 Time step (4 a.u. per step) Figure 3.5 Temperature evolutions in contacting thin layers (Ge—Ge). Dot lines refers to high temperature Ge and solid lines represent low temperature Ge 52 .0 A I _ ............................ High Temp, Ge — —— Low Temp. Ge .0 w Normalized power spectrum o to l l T T F‘I l 1 I l .0 .3 L L L I I I I I J I I I 4 6 8 10 Frequency (arbitrary unit) Figure 3.6 Atomic vibration power spectra (Ge-Ge) From Figure 3.5, we see that the equilibration time required is very short (around 1000 stepz 100 fir ). This time is even faster than the Si—Si case which has a large speed of sound. We believe the faster equilibrium is due to the larger anharrnonic phonon contribution to the heat transfer. As seen in the pure silicon cases, the vibration coupling is very strong according to the power spectra plot. The very strong vibration coupling is determined by the identical properties, like bond strength and atomic weight, of the two layers. 53 3.3.3 “atom, Silicon-Germanium This system consists of 32 Silicon atoms and 32 Germanium atoms, see Figure 3.1(b). The Silicon layer and the Germanium layer were first set to (a) 100K vs 50K, (D) 50K vs 100K, (c) 700K vs 450K and (d) 1000K vs 500K in the four simulations. The time steps, during which the temperatures are controlled in these three cases, are 10000, 8000, 8000 and 10000, respectively. The thermostats were then removed and the thermal energy equilibration process began. Figure 3.7 shows the computed temperature profiles. The vibration power spectra were calculated and are shown in Figure 3.8. 54 Figure 3.7 Temperature evolutions in the two contacting thin layers (Si-Ge) 55 (a) 100K vs 50K ..-....H“2332.223mm.-. . . . .45.. «aaaununoio. .....mwfia. .. as n .. . .......:::t . 15‘. ... .... vamu‘fiu-W: s3... sag E... .1 ...... ......3:=:==. ......ixankafi... 3......Sk\ I...:::~Q=~\\K&K W Human? can”... 33:... ...... ..........:...... ......zfififii:Essf... :3... . M 3.0.0:. I. . rff! Qs§§c sewvwvvvvwvvgg “mm... ................. ................ concoooooooogooo'o . - 3009f. “#aaaaa’a’Ooo: ......gfl. “trawflwuafiaa. item-Warm» 9 a a 3...“... ...... anq‘ssss~:- .- u o u . . .- .....nofiw.’ .o-o fizyoo.‘ 3:32:223. ....u.......§. 3:... 3.33: . thaEE§§ . _ . 3 .. ..s. 30;...» ...L. .33.... ....5... 2 .2. .8. 2» .22.... t J a .- Ikufld. .3. o... P. w 1 p sis? .... ......uli... .:~.\..::coc3oooo‘ Q. :ssssoceoeuote\::. . p m — p L . w 1 £55..§.................. xx..........¢.....r:s:r £3.83“. :2 3.. P mmw uwcauufloafioh m II 50000 III 10000 20000 30000 40000 Time step (4 a.u. per step) 56 Figure 3.7 cont. (b) 50K vs 100K 180 Ge 1500(1) 200000 1CDOOO TIME 50000 57 Figure 3.7 cont. (c) 700K vs 450K 1000 Tempegratures § Ge I 1 I I I I I I I I I l I l 50000 100000 150000 Tlme step (4 a.u. per step) 58 Figure 3.7 cont. (d) 1000K vs 500K Ge 16w 1 1. mogspmacqaoe 4000060000 Time step (4 a.u. per step) 20000 59 Figure 3.8 Atomic vibration power spectra (Si-Ge) (a) 100K vs 50K 0.2 ..s 01 .0 .3 power spectral“ Ngmallzed 8 l ‘l I j T T I 2 h. ' 3 4 5 Frequency (a 61 6 rbltrary u 7 nl) 8 10 Figure 3.8 cont. (b) 50K vs 100K .0 a .0 co Normalized power spectrum 0 .0 '0 u '0 .I r I T I ' I I I T l 1 I I T I l I T I I T i 62 Si Figure 3.8 cont. (c) 700K vs 450K ' -- Si $0.3 — Ge $0.2 - 'o 0 s E o 0.1 2 0 2 4 6 8 10 Frequency (arbirary uni) 63 Figure 3.8 cont. (d) 1000K vs 500K 0.3 ............................ Si 8 .0 M Normalized gower SPGCUUB _I. U! lrlr11fl_7—[IlITTrllr'lj'Tfifljlj a. ...... ------------------- II 'I 0- I.. II. 0__.- 2 4 6 8 I 1 110 Frequency (arbitrary unit) In Figure 3.7 (a), the temperatures do not converge to a uniform temperature as was observed in pure materials. This suggests that the two neighboring layers have little communication from a thermal energy transport point of view, although they have different temperatures. From the vibration point of view, due to the mass difference and the low temperatures, as we can see from Figure 3.8 (a), there is little spectra overlap between the two contacting layers. Also, from the phonon point of view, at such low temperatures only the low-lying phonons were excited, which also suggested the small overlap of the power spech‘a. In these ultra-thin structures, only phonons with wavelength shorter than the structue dimensions are captured. These phonons should be acoustic phonons which can have smaller wavelengths. in case (b), we switched the set temperature for the two thin layers. After a comparatively very long time (more than 20125) equilibration, the temperature came to a uniform value. We examined the vibration spectrum and found that the overlap was greater than case (a). From the above two cases, it is believed that in such ultra-thin layer junctions, thermal energy transport across the interface is affected by how the temperature difference is applied. Such behavior is similar to a P-N junction in electron transport. In Figure 3.7 (c) and (d), which correspond to 700K vs 450K and 1000K vs 500K respectively, the temperatures in the two neighboring layers, after long time of equilibration, reached a uniform value. In case (0), the time required to reach 65 the equilibrium state is about 140,000 steps (~14ps), while in case (d), me time required is 50,000 steps (~ 5173 ). It is believed that the larger anharmonicity due to the himer temperature facilitated faster thermal communication. From the power spectra plots, Figures 3.8 (c) and (d), we see that the overiap in (d) is greater than that in (c), which suggests a more significant vibration coupling in (d) than in (c). The spectra overlaps of these two mses are greater than in case (a). However, the vibration couplings in case (c) and (d) are much less than in the pure material cases which explains why a long equilibration time is required. Also, because the materials are different, there is an acoustic mismatch. Thus when acoustic phonons travel from one side to another side, there will be scattering at the interface, and only a small part of the carrier’s energy is delivered to the other side of the interface. Apparently, the strength of the interface scattering depends on temperature. 3.4 Summary and Conclusion This portion of the project simulated the manual energy transport from high temperature nanometer thin layers to low temperature nanometer thin layers. The transport behaviors in the different cases studied differ from each other. The temperature range effect and material combination effect are observed. When there is no material mismatch, which is seen in the cases in case 3.3.1 and 3.3.2, the vibrations in the two layers coupled to each other and the thermal energy was transported with little impedance. However, in the cases where the two layers are made of different types of materials, the energy transport becomes quite difficult due to the weak coupling of the atomic vibrations in the two layers. In the very low temperature case 3.3.3 (a), due to the rare overlap of vibration power-spectra of the neighboring layers, the energy transport across the interface is barely observed. However, such insulation behavior also depends on how the temperature differences are applied. In case 3.3.3 (b), where the temperatures differences were applied in a way that the phonon spectra have more overlap, thermal communication was stronger, although it was a very slow process. At very high temperatures, cases 3.3.3 (0) and (d), the vibration resonance becomes greater than that in the low temperature situation, and the increased anharmonicity contributes more to the thermal transport The thermal energy is successfully transported from the high temperature side to the low temperature side, although such processes were shown to be very slow compared to the pure material case. From the phononic point of view, the acoustic mismatch which is a result of the material mismatch at the material interface, scatters the phonon wave packet and thus leads to large transport impedance. The scattering rate becomes larger at higher temperatures and thus more energy is penetrated through the interface. But at lower temperature, the scattering is frozen at a very low rate and thus the energy transport becomes much slower. This is another explanation for the slow or even zero energy transport. From this portion of the project, we can see that the DFT-MD simulation 67 reflected physically reasonable thermal transport phenomena in nano-scaie structures. It also demonstrated its ability of “on-theofly' calculation, that is, no empirical potential was needed. It is promising that the DFT-MD can handle more complicated thermal transport problems than the classical MD is able to handle due to the empirical potential function limitation. Also, with the application of DFT-MD, it is possible to include the electrons in the simulation of thermal transport which will help solve the long-standing problem of classical MD in simulating metallic materials. 68 CHAPTER 4 Equilibrium Molecular Dynamics Calculation of Thermal Conductivity/conductance of Au-SAM-Au Junctions 4.1 Introduction of Molecular Electronics Organic electronics [3] is steadily emerging as a promising candidate for low-cost information processing, storage, display, and communication devices. Electrically conducting molecules, oligomers, and polymers can be spin coated or self-assembled onto metallic surfaces, thus eliminating high-temperature processes required for Si-based devices. Transistors, photovoltaic cells, light-emitting diodes, and even diode lasers have been synthesized using such materials [44-46]. Integrated circuits can be printed on to flexible substrates, making them particularly advantageous over Si and Ill-V based electronics and photonics. One of the major challenges in high-density Si-based electronics is power dissipation. Thermal management is slowly emerging as a key bottleneck for pushing transistor gate lengths below 50 nm. This is despite the fact that single crystal Si has a room-temperature thermal conductivity of about 140 WIm-K, which is much higher than other relevant semiconductors (SiGe - 10 WIm-lc GaAs - 50 Wlm-K; lnP - 70 Wlm-K; lnGaAs — 10 WIm-K). This situation is aggravated by the fact that low-k dielectric materials, that have then'nal 69 oondw isolate of0.1- and name came van c Thes- mm of sh end. bone inter. deSll SUci mole the Uniq that conductivities on the order of 1 WIm-K, are now being used to electromagnetically isolate neighboring devices. Now consider the fact that polymers have thermal conductivity in the range of 0.1-0.2 WIm-K, which is much lower than other most inorganic semiconductors and metals. This could potentially create major challenges in thermal management of organic or molecular based electronics and photonics. Polymers consist of chains that have covalently bonded backbone of carbon-carbon bonds. However, chain-to-chain interactions are through weak van der Waals (vdW) bonds, which make them mechanically soft and compliant These interchain weak bonds also led to poor vibrational coupling, and thereby produce low thermal conductivity [47]. To overcome this, one could use oligrners or short chains of molecules directly bound to metals or semiconductors at either end. Since metal-molecule binding (generally through Au—thiol or Si-cabon bonds) is through covalent bonds, one can eliminate weak van der Waals interactions and thereby improve vibrational coupling. This also hapmns to be the desired configuration from the electronics view point, since charge mobility in such molecules is geneally higher that that in polymers made of the same molecules, which can be critical in device performance. In addition, by modifying the chemistry, the electronic states of these molecules can be manipulated in unique ways. This provides a huge palette of potential barrier shapes and sizes that are difficult to create through bandgap engineering of inorganic 70 semiconductors. While there are several fundamental issues about charge transport across molecules and metal-molecule contacts that remain unresolved [48,49] and controversial, the opportunities for novel electronic and optoelectronic characteristics are immense, which has now made molecular electronics into an area of intense interdisciplinary research [50-52]. While most of the focus so far has been on charge transport in such molecular devices, charge-vibration coupling and thennal transport has been received little attention. Although covalent coupling between conducting molecules and the metal electrodes could improve thermal coupling, it is unlikely to eliminate the power dissipation problem. in fact, as such devices are driven for electrical or optical performance, the power dissipation problem is bound to appear. Since many of the molecular configurations and isomers have low free energy barriers, the presence of local hot spots could lead to stmctural transitions that could seriously affect electronic/photonic performance. There have been intensive works focused on the electronic and structural properties of SAM-solid junctions [53-56]. However, the studies of thermal properties of such junctions are limited, and knowledge of thermal transport in these junctions is very important to the growing fields of molecular electronics and small molecule organic thin film transistors. Ge et. al [57] measured the transport of thermally excited vibrational energy across planar interfaces between water and solids that have been chemically functionalized with SAM using the 71 time-domain thennoreflectance. Wang et. al [58] studied heat transport through SAM of long-chain hydrocarbon molecules anchored to a gold substrate by ultrafast heating of the gold. Patel et. al [59] studied interfacial thenrlal resistance of water-surfactmt-hexane systems by non-equilibrium molecular dynamics. There are two typical kinds of SAM-solid junctions, one of them is the SAM-metal junction, which has been studied a lot [60,61]. The other kind is the SAM-semiconductor junction, which is relatively new, but has become very popular [62]. In this work, thermal transport in Au-SAM-Au junctions with alkanedithiols (-S -(CH2),. “'5 -) being the SAM molecules is studied using molecular dynamics (MD). The reason such junctions are chosen is that the structural properties, including the absorption site, tilt angle, coverage and etc, have been studied thoroughly [60,61 ,63,64], and a set of reliable classical potentials for MD simulation is available [65]. In such metal-SAM-metal junctions, the SAM-metal interfaces play important roles in thermal energy transport across the junctions, especially when the system sizes are in the nanoscele [66]. In solid materials, there exist two kinds of heat carriers in thermal energy transport: phonons and electrons. Phonons are the quanta of the lattice vibrational field [18]. Phonons dominate the thermal energy transport in semiconductors and insulators while electrons play important roles in energy transport in metals. In this work, where the metal-SAM junction exists, it is difficult for electrons in the metal to tunnel through the SAM molecules, and thus electron 72 ianspt iiSVN: iniiar injur poten' Howe tansl equm to he [72,7 Firs Chet 4t2' QTOI Ve1 transportislargelydepressedbythenonmetalSAM layer.Asaresult,thegoalof this work is to calculate the lattice thermal conductivities/conductance in both in-plane directions (x- and y-directions) and out-of-plane direction (z-direction) of the junctions. For lattice thermal conductivity calculations, classical MD with appropriate potential functions has been demonstrated to be a powerful method [39,67-73]. However, to our knowledge, no work has been done to investigate the thenrlal transport properties of metal-SAM—metal junctions by MD simulation. The equilibrium MD with the Green-Kubo [37] method has been applied successfully to heterogeneous systems, such as superlattices [71] and solid-liquid interfaces [72,73]. As a result, in the present work, equilibrium MD is chosen to simulate Au-SAM-Au junctions, and thermal transport properties are calculated using the Green-Kubo method. 4.2 Equilibrium Molecular Dynamics and Green-Kubo Method Classical MD is a computational method that simulates the behavior of a group of atoms by simultaneously solving Newton’s second law of motion (eq.(4.1» for the atoms with a given set of potentials. dzr —V¢=F=MIET (4.1) where ¢ is the potential energy, I? is force, M , is atomic mass, f is position vector and t is time. By processing the trajectory obtained from an MD 73 simulal difiusio relies l fluxes. curren“ [37}. calculi therm; Simule Z-dire the ]L V3085 4.3 s Estat hYUn inCor simulation using different statistical techniques, transport properties such as diffusion coefficient and thermal conductivity can be calculated. in the equilibrium method, a system is simulated in an equilibrium state. it relies on small statistical temperature fluctuations to drive instantaneous heat fluxes. The thermal conductivity of the system can be calculated from the heat current autoconelation (HCAC) function according to the Green-Kubo formula [37]. According to Green-Kubo relation, the thermal conductivity 1: can be calculated according to eq. (2.24). Since eq. (2.26) produces heat current vectors, thermal conductivities in x-, y-, z-directions can then be calculated by one single simulation. It needs to be noted that in this work, the calculated thermal conductivity in z-direction is a combination of thermal conductivities of the materials making up the junction. it should not be considered as a material property and its value varies when the system configuration changes. 4.3 System, Simulation and Potential In any MD simulation, potential functions are critical. In this work, the well established Hautman-Klein model [65] is employed. in this model, the light hydrogen atoms of the hydrocarbon molecules are not simulated explicitly but incorporated into the carbon backbone. Their masses are added to the carbon atoms which they bond to, and forming pseudo carbon atoms. The reason for 74 amht atoms mover demOl mswn poten Wpei gdd Lenn [79]. dista are ii Spec 120 r‘cid'll i0 3 09h such treatment is that the high frequency vibrational motions of light hydrogen atoms in the hydrocarbon groups are less important than the lower frequency movements of the carbon backbone [65]. Such a treatment has been demonstrated to be a valid approach to simplify simulations and to give good results [74,75]. Bond stretching, bond bending and Ryckaert—Belbmans torsion potentials are used in alkanedithiol molecules for the bonded interactions. Morse type interactions are used to simulate the interaction between sulfur atoms and gold substrates [76,77] and the interaction among gold atoms [78]. The Lennard-Jones potentials together with the Lorentz—Berthelot mixing rule [79], gab = gagb , Gab = .2-(0'0 +0b) , are used to simulate long distance and intermolecular interactions. The potential function and parameters are listed in Table 1. All the parameters are from ref. [65] and ['76], or as otherwise specified in the table. in our simulations, a neighborlist with cutoff of less than 12.00 ft is used to speed up the calculation (In our code, the neighborlist cutoii radius is compared with dimensions of the system, and it is adjusted automatically to avoid double counting). This neighborlist is not updated alter the structure is optimized due to the fact that there should not be large atomic displacements other than vibrations about the equilibrium positions in the solid phase system. For every simulation, 5 separate runs with different random initial conditions were performed. The resulting values are averaged over the 5 runs. 75 Table 4.1 Potential Functions with Parameters Used in the Simulation 76 Potential Function forms and parameters Bond Stretching S - C C—C Us =-;—ks(r—ro)2 (a) where 2 k, =14.00224eV/;\ r0 = 1.523 furor C -—C); 1.315 finer S - C) Bond Bending U, =%k,(9—6,)2 (b) where k9 = 5.388eV/rad2 6’0 =109.5(for C — C — C); 114.4(for S — C — C) Ryckaert- Bellemans Torsion S -— C —- C — C C-C—C—C Q=ZeWflw (o where (o is the dihedral angle a0 = 0.09617eV, a1 = 0.125988eV, a2 = -0.13598eV,a3 = —0.0317eV, a4 = 0.27196e V, a5 = —O.32642eV Table 4.1 cor [innald-lones ‘ l I with i Mrenlz-Beflhelo ,tmixing rule 1 Table 4.1 cont. Lennard-Jones vL-fitilz—i-i-f—ii with Lorentzaeflhelo where a and a are determined by mixing rule [79] . mixing rule for C: a = 0.00513eV,0' = 3.9143. for S that interact with other atoms other than S : a = 0.01086eV,0' = 3.5503. for Au: 5 = 0.001691e V,c— = 2.9343. [80] for S — S interaction : a = 0.01724e V,0' = 4.2503 Morse —a(r-r ) 2 A UM =De[(1—e ”) —1.0] (e) u-Au Au-S where De = 0.4753 V(for Au —— Au)[78]; 0.380e V(for Au — S )[77] 0 -l a = 1.583 A (for Au — Au); 1 1.470 .3. (for Au — S) rm 2 3 .0242 finer Au — Au); 2.650 furor Au — S) 78 Figure —S—(( sectior ydirec molecu contain follow.1 4.2(a)). with the thus fon Wee~fol absomti atoms E Optimize. equilibn'u was repc "0'“ 50k iElmored : ”“0”th '0 M1 The structures of Au—SAM-Au junctions studied in this work are shown in Figure 4.1. The junctions consist of two gold substrates with SAM in between. —S -(CH2), —S — is used as the SAM molecule for all simulations except those in section 4.4.7. Figure 4.1 shows systems with different cell sizes in x— and y-directions. In this figure, (a) is a system with 4 molecules; (b) is a system with 16 molecules and (c) contains 36 molecules. In all these systems, each substrate contains 12 layers of gold atoms. The procedures to prepae the junctions are as follow first, one Au(111) substrate is optimized using the Morse potential (Figure 4.2(a)). Then, a number of SAM molecules are implanted on the gold substrate with the sulfur heads placed in the three-fold hollow sites of the Au(111) surface, thus forming a lattice with lattice constant of 4.99 A [60,64,811 (Figue 4.3). The three-fold site was reported to be the most stable site for SAM molecules absorption [60,81]. The molecules are then relaxed at 100K until all molecule atoms attained their equilibrium positions (Figure 4.2(b)). Finally, the other optimized gold substrate is imposed on top of the SAM (Figure 42(0)). The equilibrium cell dimension in the z-direction was found to be about 68.10 31. It was reported that the tilt angle does not change much at a temperature range from 50K to 300K [81]. In this work, temperature effect on the tilt angle was ignored so that the simulation cell size in z-direction remains unchanged throughout the temperature range (50K-350K). In MD simulations, there is a limitation on the number of atoms that can be 79 handled due to the capability of the computing hardware. Thus the thickness of the gold substrate must be limited to tens of angstroms. However, in reality, the gold substrates are often thick (several hundred microns) and their properties are close to bulk solids. The periodic boundary condition (PBC) in z-direction was used to compensate the thickness limit in our simulations. The validation of this treatment is discussed in section 4.4.2. Besides PBC, free boundary condition is also used in some simulations for comparison. PBCs are used in x- and y- directions with no exception. A simulation system which contains 4 SAM molecules and 12 gold atoms in each substate layer is defined as a unit cell (see Figure 4.1(a)). It has a dimension of 9.99 3 in x-direction and 8.652 3 in y-direction. The unit cell is expanded in the x- and y- directions to obtain the simulation supercells. Supercells with sizes of 1x1, 2x2 and 3x3 unit cells, which correspond to systems with 328, 1312 and 2952 atoms are presented in Figure 4.1. 80 Figure 4.1 Simulated Au-SAM—Au systems of different sizes: (a). a 1x1 system of 328 atoms (4 alkanedithiol molecules and 288 gold substrate atoms), (b). a 2x2 system of 1312 atoms (16 alkanedithiol molecules and 1152 gold substrate atoms), (c). a 3x3 system of 2952 atoms (36 alkanedithiol molecules and 2592 gold substrate atoms). 81 (C) (b) (a) Figure 4.2 Procedures of preparing the Au-SAM-Au simulation system: (a). One gold substrate is optimized by Morse potential; (b). Alkanedithiol molecules are implanted on the substrate and the whole system is relaxed using the potentials specified in table 4.1; (c). The other optimized substrate is imposed on top of the alkanedithiol molecules 83 Figure 4.2 cont. (b) 85 Figure 4.2 cont. (0) Flgun circle TBDIE Simul inte Vail 4.4 Figure 4.3 Absorption sites and tilt directions of SAM molecules on Au(111) surface [81] (Open circles are gold atoms on the surface of substrate. Filled circles are sulfur heads which are absorbed on the substrate surface. Arrows represent tilt directions of SAM molecules. Dashed lines form the boundary of the simulation system and dotted lines represent SAM lattice constants.) The simulation procedure is as follow: (1) All atoms started moving from their equilibrium positions with random initial velocities. (2) Nose-Hoover thermostats were applied to the system for a long time (> 25ps) to make sure that the system reached the target temperature. (3) Thermostats were then released and an equilibration period of 150ps was performed. (4) A production run was performed in which the HCAC function was calculated. (5) HCAC function was integrated and the resulted thermal conductivity was plotted to find a convergence value. 4.4 Results and Discussion 4.4.1 Defining the Value of 'l'hermal Conductivity A typical normalized heat current auto-conelation (HCAC) function profile is 87 preset overal fluctuz gold 11 me P dispel difieri out-01 IUncti SCatii furthi direc smo. plot 00m SOrr Cor rUri presented in Figure 4.4. It can be seen that fast fluctuations are imposed on the overall profile, which make the integration in eq. (2.24) non-trivial. The fast fluctuations are believed to be from the optical phonons [39, 74]., Since the bulk gold will have only acoustic phonons [82], the optical phonons are the results of the presence of SAM molecules which has optical branches in the phonon dispersion [40]. The calculated thermal conductivity profiles are different in different simulations (see Figure 4.5). For thermal energy transport in the out-of-plane direction with free boundary conditions, phonons traveling across the junction will be scattered at two Au-SAM interfaces, and they will also be scattered at the free Au surfaces at the ends of the simulation cell. Then this will further reduce the out-of-plane thermal conductivity. However, in the in—plane directions, no such scattering exists, and the thermal conductivity exhibits a smoother profile in these directions. To get the thermal conductivity from different profiles, different methods were used to find the convergence areas. The thermal conductivity profile in Figure 4.5(a), where a flat area exists, are encountered in some integrations for z-direction thermal conductivities. Similar thermal conductivity profiles were also found by McGaughey et. al.[39], who used a running average to define the convergence region. In this work, since the flat area is obvious, we define the convergent thermal conductivity values by directly averaging the values over the flat area. For profiles like Figure 4.5(b), which appeared in many z-direction integrations, a convergence area is chosen where 88 (over; high 00.0w 00nd syste short in F1 aver the calm data the overall profile fluctuates around a mean value, and the overall fluctuation (overall fluctuation refers to the fluctuation with lower frequency compared to the high frequency fluctuations) magnitude becomes a minimum. Such a convergence area was similar to the “neck region” found in ref. [39]. Thermal conductivities are found by averaging values over the convergence area For systems with free boundary conditions in z-direction, the integration profiles are shown to be like Figure 4.5(c), where the value converges to 0. For profiles seen in Figure 4.5(d), which are found in most integrations in x—, y-directions, the averaged value around the first overall peak is used. Such a method is similar to the ‘first-dip’ method [70] which was found valid for the manual conductivity calculation for crystalline p—SiC. In all the different profiles except 4.5(c), the data is averaged over no less than 5000 steps, which equals to 2.5 ps. 89 Figure 4.4 A Typical Normalized Heat Current Auto—Correlation (HCAC) Function normalized HCAC In z-directlon _L TIII'IUTTTII‘III'IIVTIIIIT o .o O) (D Ijllfillrlli .0 a. .0 to llLlllllLlllulLl .6 N normallzed HCAC In z-dlrectlon o I J 1 -0.4 d h- -l as — 3 i- -r _. -l -0.8 :- '3 -1 '— l l l l J l l L l I J 1 l l l l I l L I l l l l l 0 50000 100000 150000 200000 250000 time step (0.5fs per step) 91 Figure 4.5 Different thermal conductivity profiles obtained from the calculations: (a). profile that has a flat area; (b). profile that has an overall fluctuation minimized area; (0). profile whose value fluctuates about 0; (d). profile with a first overall peak. 92 'l‘hannol Pl‘l“l'lll"i\/I‘"V k (W/IIIK) (a) Q95 A assuage aasfi 200000 100000 150000 time step (0.5fs per step) 50000 93 AY— E\g\ v— \a:~a—+l\-1ful\l\ ...rllllltlut Figure 4.5 cont. (b) I I T T l I I I I l I 1 r I I I 1 T I k .. 4 l— ..l S? _ 4 E E g ovelal flictuation minimized area i 2 . g 3 l 6 0 I E 5 _2 1 1 1 1 I 1 1 1 1 l 1 1 I l I r 1 1 l 0 50000 1 1 200000 time step (0.5fs per step) Figure 4.5 cont. (0) 4 I ‘r r 1 I r I 1 r r 4 r l g _ _ E 2- kvalue fluctuates aroundO - E x g 0 'D S 0 E E -2 0 5 — . l- -l .4 1 L l L 1 1 1 I l time step (0.5fs per step) 95 Figure 4.5 cont. (d) 40 I I I I 1 T I I I I I I I I l I I l — r S? - _ g 30 t- - g L _ 32' - first overall peak - L _ a _ / l _>_ 13' 20 - _ g _ o ~ l '5 ~ 1 E 10 P -l 0 s i ‘ r‘ I 0 ! 1 1 l 1 I 1 J 1 1 J 1 l L 1 I 1 1 l 0 50000 100000 150000 200000 time steps (0.5fs per step) 4.4.2 real rest jun: GXC 4.4.2 Boundary Conditions As seen in Figure 4.5(c), the simulation with the free boundary condition gave thermal conductivities which were approximately 0. There are two possible reasons for such results: (1) the free boundaries add very large boundary resistance to the entire Au-SAM-Au junction which makes the conductivity of the junction very small; or (2) the ultra thin (2-5nm) substrates depressed the excitation of many phonon modes so there were only very limited modes available for thermal energy transport. However, in practice, the substrates are not as thin as several tens of nanometers. To make our simulation reflect greater reality, PBC is used in z-direction. Wlth P80 in z-direction and an appropriate neighborlist cutofl‘ (less than 12 3in this work), the gold atoms in one substrate can interact with the image gold atoms of the other substrate while they do rid interact with the image SAM molecule atoms if the interaction cutoffs are smaller than the substrate thickness. In this sense, the gold substrates are not isolated thin layers with free surfaces at the junction ends, but rather, they work as thick chunks of gold. 4.4.3 Finite Size Effect ln MD simulations, the sizes of the simulation cells usually affect the results. To investigate the finite size effects in x- and y-directions, thermal conductivity calculations were performed at 100K on 1x1, 2x2 and 3x3 systems with 12 layers of gold atoms on each substrate. Both P80 and free boundary condition were 97 [85 fro the used in the z-direction. PBCs are used in the x- and y-directions for all cases. The results are listed in table 4.2. The results from 2x2 and 3x3 systems do not differ from each other significantly when the errors are taken into account This means that, in x- and y-directions, the 2x2 system is sufficient for our simulations. In the rest of the work, all simulations use the 2x2 system. System Size 1x1 2x2 3x3 2x2 3x3 Boundary poo poo 9130 Free Free condition in z-dircction # of atoms 328 1312 2952 1312 2952 kx (W/m.K) 4911.0 17.6130 22517.4 12512.7 13.8127 ky(W/,,,.K) 4611.0 19413.4 l4.6:l:3.3 10.8116 9.5117 kz(W,,,,,K) 1210.2 1810.3 1.8104 Table 4.2 Finite Size Effect on Thermal Conductivities 4.4.4 Substrate Thickness Effect In the simulations, the z-direction is a special direction since it is the cross direction of the Au-SAM-Au junction. Energy transport has to cross two Au-SAM interfaces in this direction. In order to investigate thermal transport across the interfaces, the influence from the limited thicknesses of the substrates should be minimized. The substrate thickness effect on thermal conductivities/conductance 98 for + $0. sin Filml . Igrlr for a 2x2 system was studied by changing the number of gold layers in the substrates. All the following cases have PBCs in all three directions and all simulations are carried out at 100K To investigate the effect of SAM molecules on the substrate in—plane thermal conductivities, the system without SAM molecules is studied. The results are presented in Table 4.3. # of gold layers 6 12 18 24 36 12(no SAM) thickness L of the 39.42 68.10 96.73 125.67 182.90 68.10 0 junction (A ) k (W/mx) 9111.5 17613.0 23915.6 27.9143 36.3183 18611.6 x ky(W/,,,.K) 7.811.] 19413.4 235145198162 25317.4 19814.7 k(W/,,,.K) 1310.2 1810.3 2510.6 3811.2 4910.5 0 Z G 327158 261146 258166 302195 268130 (MW / m2K ) Table 4.3 Substrate Thickness Effect on Thermal Conductivities and Thermal Conductance To compare the z-direction results, thermal conductance should be the property to be compared, as it is mentioned in section 4.2 that the thermal 99 th llr BI Cc conductivity kz is structural dependent rather than an intrinsic property of the junctions. Thermal conductance, G, is defined by q=GoAT where q is the heat flux normal to the junction interfaces, and AT is the temperature difference across the junction. G is related to k2 by G = k2 / L (4.1) where L is the thickness of the junction. In Table 3, the thermal conductance does not emibit monotonic decrease when the substrates are thickened from 6 layers to 36 layers of gold as one would intuitively expect It is believed that the discrepancies of thermal conductance among systems with different Au thicknesses are not from the thickness size efiect but from other factors such as the difficulty of defining the convergence thermal conductivity value and the limited number of runs used to get the mean values. The effects of these factors are reflected in the error bars. The total thermal resistance of the junction can be written as a serial combination of resistances of different parts that make up the junction: Rtotal = 2 x Rsubstrate + 2 X Rinterface + RsAM (4.2) Consider the relation between resistmce and conductance G =1/R (4.3) Eq. (4.2) becomes 100 Du be 3] it 1 1 1 1 = 2x +2x +——— Gtotal , Gsubstrate Ginterface GSAM Due to the large conductance of gold substrates, the first term, “Gm... (4.4) becomes negligibly small. (In reality, bulk gold has a thermal conductivity of 318W/m-K , which leads to a conductance G of about 1.2x105MW/m2-K for a substrate with thickness of 26.2 3. The corresponding resistance R is 8.3x10'6m2-K / MW and it is negligible compared to other terms. As stated in section 1, the electron trmsporl contribution to the thermal transport is ignored since it is hard for electrons to tunnel through the SAM molecules. Due to this and the substrate thickness effects, the Au lattice thermal conductivity should be smaller than 318 W/m-K. However, even if the calculated thermal conductivities (~20 W/m-K ) in the x- and y-directions are used, the resistance comes out to be ~1.3x10‘4 mzoK / MW , which contributes only 3% of the total resistance). Wang et al [6] found that the energy transport along the SAM molecule chains was ultrafast (0.95nmlps), which suggests that the thermal resistance inside the molecule itself is very small. Henry et. al. [83] also found that a single polyethylene molecule has a thermal conductivity is higher than 100 W/m-K along the chain using an equilibrium molecular dynamics simulation with Green-Kubo formula. Therefore, the junction thermal resistance is dominated by the SAM-Au interfaces. Since there are no strong inter-molecular interactions 101 at among the discrete SAM molecules, the energy transport from one molecule to another should be weak Comparing the results of the systems with and without SAM molecules (the 6‘h case in Table 4.3), one can see that the existence of SAM molecules does not affect the in-plane thennal transport. It is also found that for the system without SAM molecules connecting the two Au substrates, the out-of-plane thermal conductance becomes 0 as expected. We then believe that the in-plane (x-, y—direction) thermal transport mostly happens in the crystalline gold substrates and SAM molecules present channels for out-of-plane thennal conduction. From Table 4.3, it can be seen that in the x-, y-directions, the supercell with thickness of 68.10 3 still suffers from the finite size effect for the in-plane thermal transport. However, if the dimension in z-direction is too large, the calculation becomes too demanding to handle with the current code. For the junctions in our work, what is really important is the out-of—plane (z-direction) thermal transport, and the 12 layer substrates have shown to be thick enough to ignore the thickness effect on the thermal conductance in z-direction (see Table 4.3). As a result, substrates with 12 gold layers are used in all the rest simulations. 4.4.5 Temperature Effect One of the topics investigated in this work is the temperature dependence of thermal conductivity/conductance. Figure 4.6 shows the in—plane results of a system with free boundary conditions. The solid lines are power fits to the discrete data. Data for simulations that use PBC in the z-direction are shown in Figure 4.7. 102 Fr 11‘. or From the fitted line, the temperature dependence of the calculated lattice thermal conductivities in the in—plane directions (x-, y-directions) shows the same trend as that found in simulations with free-bounday conditions (Figure 4.6), which declines with increased temperature, and the same trend is found in thermal conductivities of crystalline solids [70,34]. For thermal conductance in the out-of-plane direction (Figure 4.7(b)), a comparison with Wang et al [84] experimental data on Au-SAM-GaAs junctions was done (see Figure 4.10). In Wang’s work, Au-alkanedithiols—GaAs junctions were studied which are different from the Au-alkanedithiols-Au junctions studied here. In practical experiments, the alkanedithiols lattice on the GaAs substrate is not as perfect as that on the Au substrate in our work, and the molecules are not always upstraight while some of them lay down. As a result, the effective molecule-substrate contacts are not as good as the ones studied in this paper which are perfect contacts [85]. So, it is not a surprise that our data are much larger than Wang’s experimental data. Although the absolute values are not comparable, the trend of the discrete data is similar. The mean thermal conductance increases at low temperatures as the temperatures raises, and becomes almost unchanged at higher temperatures (>150k). We are not aware of any experimental thermal conductance data on exact the same junction that is studied here at a temperature range from 50K to 350K However, for metal-nonmetal interfaces, the experimental thermal 103 conductance is reported to be 8 < G < 700 MW/(mzK) [86-88]. In our work, as can be seen in Figure 4.7(b), the thermal conductance ranges from 200—300 MW/m’K which falls in the above range. Ge et al [57] reported Au-(hydrophobic SAM)«water interface to have thermal conductance of 5015 MW/(mzK) and a conductance of 100120 MW/(m’K) for Au—(hydrophilic SAM)-water interface. Wang et. al. [4] reported a thermal conductance of 2201100 MW/(m’K)foraAu-SAM junction. ltcanbeseenthatourdataisonthe same order of available experimental data. aol I r1 ITI'I I IT‘ I I I ITI TI IrI I III IfiITI I IT.-+ 25 :— .1 g _ —->—— k_x a E + k_y ' ‘_ a Ego .1 ‘” ' 1 o .. a _ a '5 1 - 15 - o ' \ j E 10 _ \_ _ o : “EH—...” .3 a 5 _ _. l. r i o I I l I I I I I I I I I I I I LII I I II I I I I I I II II I I I I 50 100 150 200 250 300 350 400 “WMWOIKI Figure 4.6 Temperature Dependence of ln-plane Thermal Conductivities of Systems with Free Boundary Condition in z-direction 104 Figure 4.7 Temperature Dependence of Thermal Conductivities and Conductance of Systems with PBC in z-direction 105 (a) Thermal Conductivities in x-,y-,z-directions 30 II'IIII'IWTIjTIl‘TITTIITTWIIjTTITTII 1; k_x 4 IL)! 01 ,/ 8 IIIl‘rIWIleIIlIIII'IjI—IIIFI \ 1 \ ~~~in . IIIILIIIIIIIIIIIJLJ ..L O I thermal conductivities (WImK) 61 UI LIIIII, 111111111I111111111I1144I11L111111I1111 50 100 150 200 250 300 350 400 temperature(K) O O 106 Figure 4.7 cont. (b) Thermal Conductance of the Au-SAM-Au Junction 888888 8 lIfIIIIIIIIIITIjII'IIIlIIIIIIlr I conductance G (MWI(m2K)) 119200 I IITIIIITI1IrTIIIIIITIIII'IIII'FIIIIIITT- ‘ fi .1 11- d -T- _‘ fi - 111111111I1111Il1111111111111I1111I111 50 100 ll 1LL1_LLLI ILLJILLIILLIIIJIIIII § 150 200 250 300 350 temperature (K) 107 Sll a: 4.4.6 Simulated Normal Pressure Effect The external normal pressure effect on the thermal conductivity/conductance was simulated by decreasing the dimension of the simulafion cell in z-direction at 100K because there is no way to apply forces on the free surfaces of the Au substrates when PBC is used in the z-direction. In this way, the junction would feel 'pressure" as the junction is compressed by the decreased cell dimension in z-direction. The z-direction dimension is decreased by 1 3. with an interval of 0.2 3. The 0.2 )1 decrease of dimension corresponds to pressure increases ranging from 10 to 350 MPa. The calculated results are shown in Table 4.4 and Figure 4.8. In Figure 4.8, where the z-direction thermal conductance presented, no pressure dependence was observed when error bars are considered. The reason for such a result could be that the SAM-Au interface resistances are the dominate factors which impeded the thermal energy transport across the junction (detailed discussions were presented in section 4.4.4). Although the junction is under pressure, due to the flexibility of the chain-like alkanedithiol molecules, the structure will adapt itself to the small dimension change in z-direction and thus the local dynamics around the SAM-Au interface does not change much. As a result, the interface resistance is not affected much. The thermal conductivities in x- and y- directions do not show any pressure dependence (see Table 4.4). 108 ‘. ihickne l67.30 i67.70 [67.10 Flgl Thickness k k k G (MW /(m2K)) x (W/m-K) y (W/m-K) z (W/m-K) (3) 68.10 22213.4 16.4130 1.8103 265144 67.90 14.5148 21.81102 2010.6 295188 67.70 15913.2 21.7165 21104 310159 67.50 14.7166 18.1163 2.2103 326144 67.30 18114.6 17314.0 2310.2 342130 67.10 18318.0 18117.2 2210.2 328130 (ligwflngl) 5 8 8 NT—rrl[rrrTrrlr—rlr—rr Thermal Conductance 8 150 Table 4.4 Simulated Pressure Dependence of Thermal Conductivities/conductance in x-, y- and z-directions TTI I'IIIITI 1 I I I I I I I L I l J I I I L I L I l I I 67.2 67. 4 67.6 67.8 68 Slmulatlon Cell Thickness In z-drection (angstrom)2 Figure 4.8 Out-of-plane Thermal Conductance vs. Simulation Cell Thicknesses 109 4.1 4.4.7 Chain Length Effect Systems with different alkanedithiol molecule chain lengths were also studied. The thermal conductivities of junctions with —S—(CH2)8 --S -, “S‘(CH2)9 ‘S" . “S’(CH2)10 “S“ were calculated at temperatures ranging from 50K to 350K The in-plane thermal conductivities at 100K are tabulated in Table 4.5, and they do not exhibit any chain length dependence. The out-of-plane thermal conductance is plotted in Figure 4.9. Solid lines are power fits of discrete data. Chemical S2 (CH )8 S2 (CH )9 S2 (CH )10 Formula kx 17.6130 14914.1 22.1182 (W/ m-K ) ky 19413.4 19313.4 16311.7 (W/ m-K ) Table 4.5 Chain Length Effect on Thermal Conductivities in x-, y-directions at 1 00K 110 mt, II'IIIIrI‘III‘IIII'IIII'IjérrfilIfrI'IITI-1 t : 2400; +3210"! _‘ fi C +$ZICH19 : .5. i n sacrum , 23507 1 O - ”300 i. 3250 “o I: r §zoo £150 .- _ 1mI-IIIIIIIIIlIIIIIIIIIAIIIIIIIIJIIIIJIIILII-1 O 50 100 150 200 250 300 350 400 temperature(K) Figure 4.9 Temperature Dependence of Au-SAM-Au Junction Thermal Conductance with Different Alkanedithiol Chain Lengths 30 18183531131111 E5201 i¥_“:::-*“' a 1 §§15 - 1" “E3. 5 0 60 100 150 200 250 300 350 Temperature (K) Figure 4.10 Experimental Data of Temperature Dependmce of Au-SAM—GaAs Junction Thermal Conductance with Different Alkanedithiol Chain Lengths [84] 111 Considering the errors, the thermal conductance values of junctions with different chain lengths are not very different from each other. However, the mean values revealed weak chain-length dependence: as the chain become longer, the out-of-plane thermal conductance decreased slightly. This makes sense if one considers the extreme cases of chain lengths equal to 0 and infinity. The thermal conductance should decrease from the maximum value for junctions with no SAM to 0 when the substrates are separated by infinitely long chains. Such a trend coincides with the experimental measurements of Wang et. al [84] on Au—SAM-GaAs junctions (see Figure 4.10). The reason that there is no obvious thermal conductance change should still be that the SAM-Au interface resistance dominates the energy transport ability in z-direction while the chain—like molecules themselves have very small resistances. It can be concluded that the limited chain length Change (from-S—(CHQ, "S- to ‘S"(CH2)10 ’S“) does "Ot have significant effect on thermal transport ability of the junctions. 4.4.8 SAM Molecule Coverage Effect To change the molecule coverage on the Au substrate, the number of molecules attached to the substrate surface was changed. With PBCs in x- and y-directions, every alkanedithiol molecule is equivalent. Molecules were deleted symmetrically so as to keep symmetries. Thermal conductivities/conductance of systems with 16, 14, 12, 10, 8 alkanedithiol molecules were calculated and presented in Table 6. 112 Coverage kx ky [(2 G (# of thiolsl (W/m-K) (W/m-K) (W/m-K) (MW/(mzK)) simulation cell) 16 17.6130 19413.4 1.8103 261146 14 13.1116 10.6122 1.8103 257148 12 20914.0 13.7148 1.6102 233126 10 14514.7 15814.] 1.3102 185126 8 17316.5 17.8142 1.010.] 148119 Table 4.6 Molecule Coverage (refers to number of alkanedithiols in a 2x2 simulation cell) Dependence of Thermal Conductivities and Thermal Conductance From Table 4.6, no trend of in-plane thermal conductivities is found. The mean value of out-cf-plane then'nal conductivities/conductance decreases with decrease of coverage as expected. In Au—SAM-Au junctions, thermal energy is transported from one substrate to another through the discrete alkanedithiol molecules. These molecules are like channels through which energy passes. When the number of molecules decreases, energy transport channels are decreased. As a result, thermal energy transport woomes more diflicult and conductance becomes smaller. 4.4.9 Vibration Coupling Analysis To investigate how efficiently the thermal energy is transported from the 113 substrate to the discrete molecules, the vibrational power spectra of Au atoms in the substrate, sulfur atoms and carbon atoms in the alkanedithiol molecules were calculated (see Figure 4.11). The calculations were done by performing Fourier transforms of the velocity autocorrelation functions. Two layers of Au atoms and 16 sulfur heads which form an interface were chosen as the samples from which the velocity autocorrelation functions were calculated. Trajectories of 16 carbon atoms bonded to the sulfur heads were sampled to calculate the vibration power spectra of the carbon atoms. It can be seen that the vibration power spectrum of the Au substrate atoms has broader peaks than the sulfur, while the spectrum of the sulfur atoms has many discrete spikes. The overlap of these two spectra is limited to frequencies lower than 10THz. It is believed that the relatively large coupling is due to the strong Morse interaction between the Au and S atoms. Although the vibration coupling between Au and S appears to be strong, the discrete SAM molecules limit the number of channels available for thermal energy transport, and thus they significantly influence the heat transfer efficiency. As a result, the Au-SAM interface presents a large resistance to the junction. It also can be seen that the overlap of S and C spectra spans over frequencies up to about 22 THz, which facilitate the thermal transport inside the molecule chains. The calculated thermal conductance values in this work are close to those of the interface conductance values of water-surfactant heads (300140 MW /(m2K) ), 114 hexane-surfactant tails (370140 MW /(m2K)) and benzene-sufactant tails (200130 MW /(m2K)) which also have lags vibration coupling [59]. The CB4-orgnanic solvent interface thermal conductance found by Huxtable et. al. [89] (10— 20 MW /(m2K)) has almost no vibration spectra overlap, and is much lower thanourdata. 0.008 '- 3|. =3..- Marl-— 53m.— 1r lllll ":3. _ :3?“ p 8 or iiillllllllillrllllr llzed Vibration Power 3 g IIIITfIIrII ‘ ~ ' h 02:9 2%.: . .r.%u.-.I.I.I.o.r.r‘r'r.n.b:.-.-.u.-'-.n.r...-.....T:...r.l.I.-:...c.o.o.l'o.l.'."0.0.'.','.'.'. — I ~ g — o — o H IliOllleallllIIOIllllh ouruuvflt r'l'r'c'a'f“ owI‘lt_ . h . .'.°.'.~'-'- If. I Norma P 8 M I w... up rt or l Ell—LIJLIIILII_I_LIILII1LLIIJIIJ_II “;:n l‘l' llllll""' lllrljr: I I I . z? '5... 2 Lil-...“-.. II::':‘ -;.I 1L»..- 20 25 30 Frequency (TH 2) Figure 4.11 Normalized Vibrational Power Spectra ofAu Substrate, Sulfur Heads and Carbon atoms 4.5 Summary and Concluslon L This portion of project calculated thermal conductivity/conductance of Au-SAM (alkanedithioI)—Au junctions using equilibrium classical MD with Green-Kubo method. PBCs were used in x- and y-directions. Both free boundary 115 condition and PBC were used in the z-direction. The effect of simulation finite size was investigated. Vibration corpling was analyzed to explore the mechanism of thermal energy transport between gold substrate and SAM. Due to the limited thermal transport channels presented by the discrete SAM molecules, thermal energy transport across the interface is not efficient even though the Aues vibration coupling appears to be strong, and thus the interface resistance is large. From the phononic point of view, the Au-SAM interfaces present scattering sites which scatter and reflect the phonon wave packets. Only a part of the phonon energy is transmitted througl the interfaces and this leads to interface resistance. Temperature dependence, molecule chain length dependence, substrate thickness dependence, pressure dependence and alkanedithiol molecule coverage dependence of manual conductivities/conductance were studied. The results show that the thermal conductance in the out-of-plane direction (z-direction) increases with temperature increase at temperatures below 150K and remains almost unchanged at temperatures above 150K The in-plane thermal conductivity displayed a bulk crystal lattice thermal conductivity behavior which decreases with the increase of temperature. It was also observed that the junction thermal conductance does not have obvious molecule chain length dependence (chain length from-S-(CH2)3 -S- to "S"(CH2)10 -S-). Substrate thickness does not seem to affect the thermal conductance across the junction. Pressure dependence is also not obvious in Au-SAM—Au junctions. 116 These three observations demonstrate that it is the Au-SAM interface that dominates the thermal transport across the junction. Alkanedithiol molecule coverage has an effect on the out-of-plane direction thermal transport The thermal conductance decreases obviously with coverage decrease due to the reduced number of energy transport channels. All the calculated thermal conductance values are between 200 and 300MW/m‘2-K which is inside the experimentally measured range of metal-nonmetal interfaces. 117 CHAPTER 5 NON-EQUILIBRIUM MOLECULAR DYNAMICS SIMULATION OF THERMAL ENERGY TRANSPORT IN AU—SAM-AU JUNCTIONS 5.1 Non-Equilibrium Molecular Dynamics and Simulation System In a NEMD simulation for thermal conductance calculation, either a constant temperature difference [90-93] or a heat flux [59,93-97] is imposed by altering the atomic dynamics in localized heat sink and source regime. In the constant temperature difference method, the temperatures of the sink and source regions are themostated at the target values by rescaling the atom velocities in these two regions. The energy differences before and alter the velocity rescaling processes in both the sink and source regions generate heat fluxes. Although the velocity rescaling is an abrupt disturbance to the simulation system, it was proved to be a valid algorithm which does not significantly modify the local thennal equilibrium of the sink and source region [90,92,98]. In the method which imposes an energy flux to the simulation system in this work, the velocities of the atoms in the sink and source regions are rescaled so that the same amount of energy is added to the source and taken out from the sink. The amount of energy difference before and after the rescaling is determined by the energy difference between the fastest atom in the sink region and the slowest atom in the source region. In this portion of project, the second algorithm was mainly used while the first one was also used 118 for comparison purpose. In a NEMD, after the simulated system reaches a steady state, the thermal conductance can be calculated by Foun'er’s law (eq. 5.1). J = GoAT (5.1) where AT is the temperature difference of the heat sink and the heat source, J is the heat flux which is defined as the amount of heat energy transfer rate per unit area perpendicular to the direction of the heat flux and G is the thermal conductance in the heat flux direction. A typical set-up of the simulation system in this work is shown in Figure 5.1. The simulation system consists of three Au substrates with SAM connecting them. The SAM is made of alkaneditiols (“S " (CH 2 )n "’ S "" ). Periodic bounday conditions (PBC) are used in all three spatial (x-, y-, z-) directions. Wrth PBC in z-direction, the left substrate and the right substrate make up a united substrate which is exactly the same as the substrate in the center. Four Au layers at the center are chosen as the sink region and two Au layers at each ends are chosen as the source region. Atoms in these two regions are subject to velocity rescaling while all the rest atoms move freely. For the system shown in Figure 5.1, eq. (5.1) can be expressed as dE J = ZA-At =GAT (5.2) where Ar is the time interval between two consecutive velocity rescalings, A is 119 the area of the system perpendicular to the heat flux and dE is the energy taken out from the heat sink and put in to the heat source. Calculation errors are obtaimd by analysis the error propagation relation through eq. 5.2: e6 = "lo'dEIZ +|0'M|2 xG (5.3) where eg istheabsolute errorof G, 0'“! and 0M are relative errorsof (E and AT respectively. «...—- O'uhiw“ Jllblbss‘. I'I".‘ , _. ‘\.l’J-" .00.. E Figure 5.1 ATypical Set-up of the Simulation System In this part, the same unit cell is used as defined in chapter 4 section 4.3. The simulation supercells are obtained by expanding the unit cell in x and y directions. The system in Figure 5.1 is a 2x2 system. The potential functions and parameters are the same as those in chapter 4 section 4.3. The simulation procedure is as follow (1 ) All atoms started moving from their equilibrium positions with random initial velocities. (2) Thermostats were applied to the whole system to make sure that the system reached the target mean 120 temperature. (3) Thermostats were then released and an equilibration period was performed. (4) Thermostats for heat sink and source were switched on. (5) After the system reached a steady state, a production period during which the heat current and temperatures at different 2 coordinates are recorded. (6) Thermal conductance was calculated according to eq. (5.3). For all the simulations, neighborlists with radii less than 12 A were used to speed up the calculation and the time step was set to 0.5 fs. 5.2 Results and Discussion 5.2.1 Finite Size Effect To investigate the finite size effect in all three directions, interfacial thermal conductance values of systems with different cross section sizes and different substrate thicknesses are calculated. Twle 5.1 shows the results from calculations using the constant temperature difference method, and Table 5.2 shows results from calculations with the heat flux method. 121 Table 5.1 Thermal Conductance of Systems with the Constant Temperatures Method 122 Source Sink Cross Cross section Junction G temp. temp. section arcadxfi.) lengthUok) (MW/(1)1210) (K) (K) size 130 70 2x2 19.98 x 17.30 136.20 340 130 70 2x2 19.98 x1730 251.34 352 130 70 3x3 29.97 x 25.95 136.20 344 140 60 2x2 19.98 x1730 136.20 330 150 50 2x2 19.98 x1730 136.20 342 123 Table 5.2 Thermal Conductance Systems with Heat Fluxes Imposed 124 Cross Cross Junction G section section area length (3,) (MW /(m2K)) Size (Ax A) 2x2 19.98 x17.3i 136.20 348 :1: 80 2x2 19.93x17.3i 251.34 360i85 3x3 29.97 x 25.95 136.20 352 :l: 80 125 difiere in cm same ydire show thick: subs excit in th is al only Sysl lEm hea su§ COD. the A13 In Table 5.1, the interfacial thermal conductance values from systems with different sizes are close to each other. For case 1 and 3, which are only different in cross section area, the interfacial thennal condumnce values are almost the same. This suggests that the finite size effect is not important in the x- and y-directions for the 2x2 systems. The steady state temperature profile of case 4 is shown in Figure 5.2 (a). For case 1 and 2, which are only different in substrate thicknesses, the interfacial thermal conductance of the system with thicker substrate has a slightly higher value. Since thicker substrates allow phonons with larger mean free path to be excited in the simulation, more phonons are involved in thermal energy transport in the system with thicker substrates than that with thinner substrates. However, it is also found that the thermal conductance difference between case 1 and 2 is only 3.5%. As a result, we can imore the substrate thickness effect when the systems have lengths of 136.23.. In Figure 5.2 (b), which shows the steady state temperature profile of case 2, the substrate temperature profiles away from the heat sink and source or the SAM-Au interfaces are very flat. Such an observation suggests that the thermal conductance in the Au substrates is very large compared to the thermal conductance of the Au-SAM interfaces where large temperature gaps exist. The temperature profiles are nonlinear near the sink and the source regions due the artificially altered dynamics by the velocity scalings. Also, the temperature profiles become nonlinear near the interfaces which 126 suggesttl'latthedynamicsoftheAuatoms neartheinterfacesaredifferentthan those inside the substrate due to the surface effect. It should be noted that the temperature nonlinearities in Au substrates do not affect the validity of applying eq. (5.2) on the Au-SAM interface since only the temperature difference between two points are involved to calculate the interfacial thermal conductance. Simulations with different heat sink and source temperatures were also perfomd (see case 1, 4 and 5), and the calculated thennal conductance values are not very different fiom one another. This suggests an almost linear relation between the heat flux and the interfacial temperature difference, which demonstrates the validity of the Foun'er’s law on the AuSAM interface. The simulations with the heat flux method yield results (Table 5.2) that are similar to those from constant temperature difference method. Although these two nonequilibrium algorithms yield similar results, the large fluctuations of heat fluxes in the constant temperature difference method (Figure 5.3 (a)) lead to relative errors larger than 1. Moreover, it is very diflicult to make sure that the simulation system reaches a state in which the mean energy taken out from the heat sink equalsthatputintotheheatsource. Theheatfluxmethoddoesnotsufferfromthe above two problems (Figure 5.3(b)), and thus the rest of the work uses the heat flux method. As discussed above, a 2x2 system with thickness of 136.20 A is enough to ignore the finite size effect, and such systems are used throughout the rest of the simulations. 127 Figure 5.2 Steady state temperature profiles: (a). Temperature profile of a 2x2 0 system with a junction thickness of 136.20 A. (b). Temperature profile of a 2x2 0 system with a junction thickness of 251.34 A . 128 (a) Temperature (K) Source Sink Source 180 -' I 1 I 1' T ‘ I I l I I I ' r 1' T l I f I I _ : ... 140 :4 . - k. - I i‘ i A“ 100 T A 44‘“ J. 7 _. . y A j _, “ “‘ ' — so ; ‘.‘ ‘ - ; _l 1 I l 1 J 1 J. L i l J L l I 1 I l l 1 l l 1 — o 20 40 60 so 100 120 129 Figure 5.2 cont. (b) Temperature (K) 40 Source Sink Z (angstrom) 130 Figure 5.3 Heat fluxes: (a) heat flux of heat source and heat sink in a simulation with constant sink and source temperatures, (b) heat flux of a simulation with the heat flux method. 131 (a) 3E+11 2E+11 lll'rl'TIj— 1E+11 -1E+11 Heat Flux J (JImAZIs) O -2E+1 1 IIIITIIIjIV -——--- J-Il ummmmmml m .II ii: lIih-JII- W'IIIII r 351-1 1 2E+1 1 151-11 I‘lllllllll -1E+11 -2E+1 1 [Ill—iji -+11 E 0 132 I 800535+11 Figure 5.3 cont. (b) 2.4E+1o .- I f I 1 I I I 1 I 1 1 l I I I C 2 2.2E+10 - 4 A " _ w t - \ . f" 2E+10 F -_ a i - \ - .. C 1.8E+10 - j >< I F .3. 1.0910 ": 4 u. ‘3 MEI-10% W h 0) :1: . N 1.2E+10 :- 1E+10 :- H C l l l L I I l I l L A I l l I 83090 20000 40000 60000 8000 time step 133 5.2.2 Temperature Effect For NEMD simulations in which there are temperature difierences, an equilibrium temperature of the system is not defined. Systems with different mean temperatures are studied. The thermal conductance values are plotted against the mean temperatures, as shown in Figure 5.4. 5w—I1TII—erI1TFTITTIIIll1IIIIIIIIIIll'llllllIll- . 1— .4 2500 - “ "‘ j N _ -_ E. : 3 £450 :. .. E ; r l o " .. g t __ I F i - g4“) l" i ...1 '0 - -1 S350 :. l- .1 0 — T - 1 E300 P ‘“ J“ “ i 9 .1: t —— - l- J j 2w 1 I I I-‘l—I I I I l I I I I l I I I I l I I I I I I I I I l I I I I l I I I I l I I I I 0 50 1 00 150 200 250 300 350 400 450 Temperature (K) Figure 5.4 Au-SAM interfacial thermal conductances at different mean temperatures. The interfacial thermal conductance increases with temperature increase when the mean temperatures are below 250K The thermal conductance values do not have significant changes at temperatures above 250K Such a trend of 134 thermal conductance is similar to that found in the experiments by Wang et. al. [84] on GaAs-SAM—Au junctions and from our previous EMD study of the same Au-SAM-Au junctions [99]. The plateau in Figure 5.4 appears at a higher temperature than those in the aforementioned two references. it should be noted that in the nonequilibrium steady state, the two inner Au-SAM interfaces are at temperatures lower than the system mean temperature, while the outer two interfaces are at temperatures hidier than the system mean temperature (see Figure 5.2). As a result, we estimate that the thermal conductance stops increasing with temperature at a temperature around 200 K, which is close to that found in the EMD study [99]. it should also be noted that the GaAs-SAM junctions and the Au-SAM junctions have different degrees of anharmonicity. Since the anharmonicity is a factor that will counter afiect thermal transport across the junction (will be discussed in the next paragraph), it is easy to understand why the plateau appears at different temperatures in different junctions. In the experiment on Au.SAM junctions, which was concluded by Wang et. al [4], an estimated interfacial thermal conductance of ZZOiIOOMW/(mzK) at 800 °C was reported. Considering the approximations made in estimating the thermal conductance in ref. [27], our high temperature data (400 :t 120MW /(m2K) ) are in reasonable agreement with their estimated value. If we assume that the thermal conductance of the Au-SAM-Au junction is dominated by the Au-SAM interfaces, the junction thermal conductance in this work is approximately 135 a200i60MW /(m2K) at high temperatures. In ref. [100] by Segal et. al., where a quantum mechanical model was used to predict the thermal conductance of an alkane chain attached to two electrodes, a thermal conductance of about 3.5xlO'“W/K per alkane chain at 300K was reported. Considering the area per chain to be 2.l6> (5.5) The VDOS’s of surface Au (the top layer of Au substrate at the Au-SAM surface), 3, c1 (the 1" carbon atom from the sulfur head), c2 (the 2"cl carbon atom from the sulfur head) and C4 (the 4‘“ carbon atom from the sulfur head) atoms at temperatures of 50K, 150K, 250K and 350K The VAF’s are obtained from 136 equilibrium runs at the aforementioned temperatures. The results are presented in Figure 5.5. The VDOS’s are not only proportional to the population of vibration modes but also proportional to temperature since they are proportional to the square of the velocity [101]. To make the VDOS’s of different elements comparable, VDOS’s of different elements are weighted by their respective mass. 137 Figure 5.5 vooss of surface Au (1" column), s (2“I column), one" column), 02(4th column) and 04(5‘h column) atoms at temperatures of 350K (top row), 250K (2"d row), 150K (3" row) and 50K (bottom row). C1 refers to the 1" carbon atom from the sulfur head, C2 refers to the 2"L1 carbon atom from the sulfur head and C4 refers to the 4‘" carbon atom from the sulfur heads 138 35 383on 139 It can be seen that the populations of vibration modes at all frequencies of all elements grow with temperature increase. As a result, more phonons are involved in thermal energy transport, which leads to the increase of thermal conductance with temperature. However, this does not explain why the thermal conductance stops increasing at high temperatures. To facilitate analysis, we divide the molecular (including 8 and C atoms) vibration modes into three regions: low-frequency (LF) modes ranging from 0 to 15THz, intermediate-frequency (IF) modes ranging from 15 to 30THz aid high-frequency (HF) modes with frequencies higher than 30THz It is found that the HF modes are not excited at temperature up to 400K, so only the LF and IF modes are visualized. In Figure 5.5, both the substrate Au atoms and the SAM molecule atoms have largely populated LF modes, leadingtoaresonancetypeofthennaltlmsportbetweentheAu substrates and SAM molecules. It also can be seen that the LF modes extend overall all the molecular elements, suggesting a highly delocalized feature of thesemodes. ThelFmodescouldalsotransportthermal energy, however, the direct coupling between the IF modes to the Au substrate is not possible since there are no vibration modes in the Au substrate with frequencies ranging from 15 to 3OTHz. However, the anharmonicity makes the LF and IF phonon-phonon interactions possible, and thus felicitates energy communication between LF modes and the IF modes. One can picture the following themlal energy transport mechanism: some 140 BI‘IE' pho mol Iunl H01 cor tun thrl usi int be m. tr; 1h lh energyoftheLFmodesistransferredtolFmodesatoneinterfacedueto phonon-phonon interactions. Then IF modes carry the energy through the SAM molewlesandreleasetheenergytoLFmodesattheotherinterface. Sucha tunneling-like transport decreases exponentially with chain-length increase [100]. However, in chapter 4, EMD simulations on Au-SAM-Au found that the themial conductance did not depend on molecule chain-length, which suggests that the tunneling-like transport does not play an important role in the thermal transport through the Au-SAM-Au junctions with SAM molecular chain-length of 8 or more carbon atoms. To increase the thermal conchctance of the solid-SAM junctions, using substrate materials which can utilize the IF vibration modes for resonance type transport should be an effective strategy. In our model, the potentials used to describe the SAM intramolecular interactions are mainly harmonic interactions (see Table 4.1). Although the LJ and Ryckaert-Bellemans Torsion potential could result in anharmonicity, the effect is believed to be small because these interactions are very weak Segal et. al [100] compared the full alkane force field (with anharmonicity) and the pure harmonic model, and it was found that the limited anharmonicity does not affect the thermal transport along the alkane chain. So, we believe themlal energy transport inside the SAM molecule is dominantly harmonic (ballistic transport), which transports thermal energy efficiently. This explains why the temperature difference inside the SAM molecule is very small (see Figu‘e 5.2). However, the Morse bonds, which 141 are used to describe the Au-Au and Au—S interactions, are not harmonic, resulting in anharmonicity. The anharmonicity which scatters phonons counter affects the thermal energy transport from Au into SAM molecule [102]. It should also be noted that the anharmonicity increases with temperature increase, and this explains why the thermal conductance does not keep increasing at high temperatures. From Figure 5.2, we can see that the largest temperature drops exist at the Au—SAM interfaces. We can also see from Figure 5.5 that there are strong couplings between LF molecule modes and the Au substrate. It should be noted that the ratio of the number of the Au substrate atoms to the numbers of the SAM molecule atoms is 3.6. As a result, the SAM molecular spectra work as bottlenecks, which have much smaller vibration modal populations compared to that of the Au substrate (one could imagine multiplying the Au spectra by 3.6 and examining the overlap between Au spectra and SAM molecule spectra). As a result, the populations of the SAM molecular LF modes that can resonate with the substrate vibration modes are limited, and this leads to relatively large interfacial thermal resistance compared to the resistances of Au substrates and SAM molecules. 5.2.3 Simulated Pressure Effect The normal pressure effect on the thermal conductance is simulated in this study by decreasing the dimensions of the simulation supercells in the 142 z-directions. In this way, the junction would feel “pressure” as it is compressed by the decreased cell sizes. Simulations were performed at mean temperatures of 100K and 300K The calculated Au-SAM interfacial then'nal conductance values are plotted in Figure 5.6. It is found that the then'nal conductance does not depend on the simulated external pressue. From structural point of view, although the junction is under pressure, due to the flexibility of the chain-like SAM molecules, the structure can adapt itself to the small dimension changes in the z-direction so that the dynamics at the SAM-Au interface are not changed much. Thus, neither the vibration modes nor the strength of the anharmonicity should change much. As a result, the phonon transport in the junction remains the same, and this leads to the same thermal conductance. 143 Figure 5.6 Interfacial thermal conductance versus simulation supercell thicknesses. 144 Jherrgal Conductance (MW/(11112 K)‘; a 8 a e a a 8 § § Supercell Dimension in Z-direction (angstrom) 145 I 1 I I I I I Tfi r I T I I r I 1 r I I T f r I I I g + 100K + 300K ; F '1 :- 7— 11‘ 7’ T" € C : : _[ ._ __ _ C F —— l "' ‘1 ; 2 *- 1 I: I -_ 4 L :1 ’— I l l I I I I I I I I J I I I I I I I J I I I I 4 J. 134 134.5 135 135.5 136 135.5 5.2.4 Coverage Effect The SAM molecular coverage effect on thermal conductance is also studied by changing the number of alkamdithiol molecules on the substrates. With PBC in x- and y-directions, every alkanedithiol molecule is equivalent in position. Molecules are deleted symmetrically so as to keep symmetries. Au—SAM interfacial thermal conductance of the M systems with 16, 14, 12, 10, 8 alkanedithiol molecules in each SAM are calculated at a mean temperature of 100K and the results are presented in Figure 5.7. 5w _ f r T l l T l 1 T I I l I I T 1 f I l I l l T 2450- -§ : T 1 E400 E- .2 3 : " : £350 1: $ .: d! .. .. o — s c .. —_ _ g3“) _-_- f _ g E .- —— I l— — 8250: T : .. '_' ‘F I TL 1 E200 -_- j .3 : I __ - 1-150 _- ".1. b I I I l I I I I I I I I I I l J I I L I l L J 1006 8 10 12 14 16 18 Number of SAM Molecules Figure 5.7 Coverage dependence of Au-SAM interfacial thermal conductance at 1 00K 146 In Figure 5.7, it is seen that the thermal conductance increases when the number of the SAM molecules increases. In the Au-SAM-Au junctions, thermal energy is transported from one substrate to mother through the SAM molecules. These molecules work as “channels" through which energy passes. When the number of molecules decreases, the number of thennal transport channels decreases. Due to the weak intermolecular interactions (van de Waals (vdW) interactions), the intermolecular thermal energy transport is difficult. This makes the "channels' relatively isolated. Since each molecule has equal capacity of transporting themlal energy, the thermal conductance is almost a linear function of the number of molecules. The observations in Figure 5.7 can also be explained from the vibration point of view. As discussed in Section 5.2.2, the populations of LF molecular vibration modes work as bottlenecks which limit the thermal transport efficiency. As the number of molecules decreases, the total population of the LF molecular vibration modes to couple with Au substrate decreases, leading to the decrease of thermal conductance. 5.2.5 Au-SAM Bond Strength Effect As discussed in Section 5.2.4, each SAM molecule works as an isolated themlal transport 'channel”. As a result, the contacts of the molecules to the Au substrate are important parameters which could influence the energy transport efficiency of each “channel". To study the effect of the Au-SAM bond strength on thethennalenergytransportacrosstheinterfaces, thebindingenergy(De)ofthe 147 Morse potential, which describes theAu—S bond, was increased by 10%, 20% and 30%. The calculated interfacial thermal conductance data are plotted in Figure 5.8. It can be seen that the interfacial thermal conductance increases as the Au-SAM binding energy increases. The VDOS’s are calculated for the surface Au atoms and the S atoms from the MD run with the original bond strength (l-De) and that with 1.3-De (Figure 5.9). It can be seen, that the vibration modes in the Au spectra with frequencies lower than 5THz shift to the left and those with frequencies higher than 5THz shift to the right when the bond strength is increased. The 8 spectra also exhibit similar trends. It is easy to understand the spectra shifts to the higher frequencies due to the increased bond strength, but it is not intuitive to believe that some modes should shift to lower frequencies. The spectra shifts can potentially change the LF modes resonance between Au and S atoms, and thus change the interfacial thermal conductance. However, because the shifts are so tiny, the thermal conductance change due to the spectra shift is expectedtobesmall. Ontheotherhand, duetotheincreasedbondstrengths, the motions of the surface Au atoms and the S atoms are more confined around the potential minimum (equilibrium position), and thus the anharmonicity at the interface is decreased. This leads to smaller phonon scattering at the interface due to the anharmonicity, hence higher them'lal conductance. As a result, we believe that the thermal conductance change due to bond strength change is mainly because of the change of the anharmonicity from the Au-S bond. Some 148 modesshiflingtolowerfrequenciescouldalsobearesultoftl'lechangeofthe anharrnonic scattering at the interface. 550 PT I l i I I j j I T T I 1 I r— I r l I T t : 9500 f "T j N : % s. . - ~~450 - —- ‘ - 3 3 F j g _ -- _ o 2. ‘ g400 : fi 1 $350:- I 1' 1' ~: C b '1 O +- 4 2300 3- .1 a _ ..- E : _I_ J‘ .1. £250 7 1 2m -I I I I J I 4 L J J I l I I I I J l I I1 1 1.1 1.2 1.3 Binding Energy (De) Figure 5.8 Binding Energy Dependence of Au-SAM Interfacial Thermal Conductance at 100K De is the binding energy of the original Morse potential. 149 Figure 5.9 VDOS’s of the surface Au atoms and the S atoms for the original binding energy and 130% binding energy. 150 6 4 Frequency (TH z) 3.5 eases moo> 4.1—.4qq_«__l_.jfi1a___._nn"a r l T ,1 r... ......... m m r e I D I m*— I .I \I . Wu 2 ... 5H I _ _ 1T - s s ( r . ....... m. . e r u I I 0 q s 1 r ...: m l . ........ .. ......usatssfim." u . . . .1.qu F I nnnnnnn bane-nun"..nuummw.u.. l 5 _r_rP+___—h___—pL__Ps___P__. 0 m m m m m m o m 0 0 0 0 w a m m m was rages 89 151 5.2.6 Heat pulse To study the thermal energy dissipation in the Au-SAM-Au junctions, a heat pulse is introduced to the source region after the system reached equilibrium state at 100K The heat pulse is imposed by sealing the temperature of the source region to 800K To visualize the temperature profiles at different instaices while minimize the influence of temperature fluctuations, time steps are divided into blocks with each block containing 1000 time steps (=500fs), temperatures a‘e then averaged over the 1000 steps. Snapshots of temperature profiles are presented in Figure 10, and the time block in which the heat pulse was imposed was set to t=0. 152 Figure 5.10 Snapshots of Temperature Profiles 153 800 A =5. . . a t ‘ 3 . . ‘ g 400 .. '1 P. F O. A . ‘e‘e. ..‘Aa‘ e a h. M *- ‘ggggmlb‘AAldblbfllbtfifi ”......‘Ab‘fi .‘Q.....“ 0 I l I A I 1 I l I L I A l l 1 l I Z (angstrom 800 I I 1' I I I Y T I I I I I I I I I I g _ t=6.5PS t=1 1.5ps t=19.0ps ‘ 9 3 g 400 ~ -_ « ~ « m r O. b A. ‘1 E ..“..A ‘ .‘ .A “ten... .“.‘1i“ebb u‘..‘. .2 ”a....“..‘... “.......g‘..~ “geese....ss‘ 0 4 1 1 1 r r L I r 1 1 I I 1 1 I 1 I Z (angstrom 800 I I 1 r I I T m I I I I I I I I I T 3 “27-595 .. t=36.5ps t=74.0ps 2 g 400 - . I 4 ~ . h 0 Q. g L....‘¢“AAAQA“.¢.“Q A‘.“fl‘6e‘sesbsee“b“‘.‘ :A“‘w“.5l“hs‘~““‘ l- 0 I 1_ I I J L I 1 I J I 1 1 1 J 1 I J Z (angstrom) 154 From Figure 5.10, it can be seen thatthetime requiredforthe heatto spread over the Au substrate is less than 1.5 ps. The time required for energy to cross the interfaces is much longer, which is betwmn 36.5 ps and 74.0 ps. Temperature gradients inside the SAM molecules are always small in these snapshots, suggesting fast heat dissipation inside the SAM molecules. These observations prove that thermal resistances of the Au substrates and the SAM molecules are much smaller than the Au-SAM interfaces resistance. 5.2.7 Further VDOS Analysis VDOS’s are also calculated for the Au atoms at the Au—SAM interface and those locate inside the substrate (Figure 11). To identify the locations of the Au atoms, the 12-layer Au substrate is divided into 6 slices with the first slice being the source region. The VDOS’s of the Au atoms of the third (Au_3) and the fourth (Au_4) slices are calculated, representing the inside atoms. The VDOS’s of the single layer Au atoms at the interface are calculated as the surface atoms. The VDOS’s are calculated from the trajectories from an equilibrium MD run at 150K By comparing the spectra, it can be seen that the vibrational features of the inside Au atoms are almost the same (Au_3 and Au_4), while the dynamics of the Au atoms at the interface are obviously different. Compared to the VDOS’s of the inside Au atoms, the VDOS of the surface Au atoms has two obvious peaks around 2THz, and the modal populations from 3 to 6THz decreased. Overall, the VDOS moved to the low frequency side when the Au atoms are at the interface. 155 Theeffectivetotal bondstrengthsoftl'tesurfaceAuatomsare reducadduetothe redumd number of bonds compared to the atoms inside the substrate. This leads to the reduced intensity of higher freqmncy vibrations and the emphasis of the lower frequency vibrations. Such change of spectra leads to additional boundary thermal resistance, and this is reflected in the nonlinear substrate temperature profiles near the interfaces in Figure 5.2. 1m I I I 1 I I I I —I T I l I l I I I I I 1’ Au_3 d ------- Au_4 “ Au_su'faco ‘ 80000 i- q :3 " : i J 5 * fl :1 l :1. ‘ am - u '1 :: - s * a = , e ' ' ' . '1 it: ‘ g r .1‘ : i 1 ID 40000 :- 1‘. , it .. 8 i > ’ f‘ _ 1 20000 '- r _ 1 fl ”’r I r j‘1 0 I. 1‘ l'l l l 0 2 4 6 8 10 Frequency (TH 2) Figure 5.11 VDOS’s of the Au atoms at different locations of the substrate 156 Since the areas under the VDOS lines are proportional to the temperatures, wecan comparetl'tetemperaturesoftheLFmodesandthelFmodesofthe SAM molecules separately. Such a frequency dependent temperature of a certain species a is defined in the same way as described in ref. [101]: w+Am I Dame (a) )da) Ta ,eq 0+Aa) 5.6 Da 3,4 (a) )da) ( ) T (0))= to Da ,ne (0)) and Da ,eq ((0) are the VDOS’s from non-equilibrium and equilibrium runs, 121(0) is the frequency dependent temperature of the non-equilibrium system, and Ta,eq is the temperature in an equilibrium run. The temperature of the LF modes is obtained from eq. (5.6) with integrals taken from O to 15THz, and the temperature of the IF modes is the integral from 15THz to 30THz. The non-equilibrium VDOS’s are calculated from a run with a mean temperature of 100K The temperatures are calculated for the C atoms in the SAM molecules, and the results are shown in Figure 5.12. The x-axes represent the numbering of the C atoms which is defined in the way as described in Section 5.2.2. It is very interesting to find that the temperatures of the LF modes have a symmetric profile with a convex at the center. The higher temperatures of the 157 centralCatomssuggestanextraheatsourceattl'tesepositions.Thisshouldbe the result of the delocalized LF modes which is reported to delocalize over four to five mrbon segments [4,100]. As a result, we believe some energy is transported to the central C atoms of the SAM molecules directly from the Au substrate without passing through other atoms in the molecule. Wang et. al. [4] also concluded that Au layer did not transfer its heat to an individual atom at the base of the SAM. However, we cannot confirm such a finding given the information from Figure 5.12. In the IF modal temperatures, the atoms that are closer to the hot substrate (solid line) have slightly higher temperatures than the atoms close to the cool substrate (dashed line). It seems that the temperatue of the IF modes are weakly influenced by the substrate temperatures. Since there is no direct coupling of these modes to the substrate, the aforementioned finding is believed to be related to the different strengths of the anharmonicities of the two substrates. Overall, both temperatures of the LF modes and the IF modes are largely decoupled with the substrate temperatures. This is because that the ballistic energy transport inside a SAM molecule has a much faster rate than the interface energy transport. One could picture the thermal transport across the Au-SAM-Au junction as follow: thermal energy from the high temperature substrate is transported to the low temperature substrate through the SAM molecules (not necessarily along the molecule chains), meanwhile, energy redistributes inside the SAM at a much 158 faster rate. 120 L i low-W T j ——O—— intenmdde-freqmncyT ] @110 b d 1’ I I 3 _ g1 _ j 00 o P t.- .. 90 - 4 : ————a 80 C1 CZ C3 C4 C5 C6 C7 C8 Label of C atoms Figure 5.12 Frequency dependent temperatures of the LF and IF modes of the C atoms in the SAM molecules. The solid line is the surface temperature of the high temperature substrate and the dashed line is the surface temperature of the low temperature substrate. 5.3 Summary and Concluslon The present work studied interfacial thermal conductance of Au-SAM(octanedithiol)-Au junctions using NEMD simulations. Temperature dependence, simulated pressure dependence, SAM molecule coverage dependence and bond strength dependence of the Au-SAM interfacial thermal conductance were studied. It was found that the thermal conductance increases 159 with temperature increase at temperatures below 250K, and it did not change much at temperatures from 250K to 400K The increase of thermal conductance was due to the increased populations of phonons which were involved in the thermal energy transport While at high temperatures, the anharmonicity depressed the further increase of the thermal conductance. Pressure dependence was not found for the interfacial thermal conductance. Changing SAM coverage was found to be an effective way of changing the efficiency of mermal energy transport across the Au-SAM interfaces. The thennal conductance decreased obviously with coverage decrease due to the reduced number of SAM molecules, which acted as energy transport “channels”. Bond strength also had an impact on the interfacial thermal transport. Stronger Au-S bonds lead to smaller anharmonicity at the interface, and thus facilitate thermal transport. VDOS’s were analyzed to explore the mechanism of the Au-SAM thenhal energy transport and energy transport inside the SAM molecules. The low VDOS’s of the SAM molecules worked as bottlenecks, and thermal energy transport from Au to SAM molecules was not efficient. Thus the interfacial resistance was large. The VDOS’s also showed that the vibration modes were largely delocalized in the SAM molecule, and the harmonic (ballistic) transport redistributed energy efficiently inside the molecules. The junction response to a heat pulse demonstrated that the heat dissipation inside the SAM molecules and inside the Au substrates were very fast, but the energy transport across the Au-SAM 160 interface took a much longer time. These observations indicated that it was the Au-SAM interface that dominated the thermal transport across the junction. Our calculated thermal conductance agreed reasonably well with available experimental data and other theoretical predictions. All the calculated thennal conductance values were between 150 and 450MW/(m2 -K), which was inside the range of the experimentally reported then'nal conductance of metal-nonmetal interfaces 8 < G < 700 MW /(m2K) [86-88]. 161 CHAPTER 6 DENSITY FUNCTIONAL THEORY STUDY OF THIOL ADSORPTION ON GAAS(001) SURFACE 6.1 Introduction In the previous two chapters, thermal energy transport in AuSAM (dithiol)-Au junctions was studied using classical MD. The prerequisite of such classical MD simulations is a set of reliable empirical potential functions. Besides SAM on Au substrate, another popular molecule-solid junction is SAM on GaAs substrate. However, the MD simulation of thermal energy transport across SAM-GaAs interfaces is cumbered by the lack of a reliable empirical potential to describe the SAM-GaAs bond interactions. As discussed in chapter 2, the ab-initio method employing DFT is able to simulate any system without the need of empirical potentials, because the forces acting on the nuclei are calculated by solving the structures of the electronic systems. Ab-initio method can calculate the bond interactions with great accuracy. It has also been pointed out that the drawback of such a method is the demanding computational load. Such a disadvantage makes the ab-initio MD unsuitable to simulate thermal energy transport in large systems (thousands of atoms). Although there have been many experimental studies focusing on the thiol 162 8C lh adsorption on the GaAs (001) surfaces [103,104], theoretical simulation on the thiol-GaAs interaction is very limited. Only available theoretical simulation is by Voznyy et. al. [105], who used the DFT to study the structue, bonding nature aid binding energy of thiols on As-rich GaAs(001) surfaces. In classical MD simulations, the potential functions are usually chosen empirically, and the parameters are often fitted to either experimental data or results from higher level calculations, like drinitio calculations. However, the bonding data for thiol-GaAs from experiments are not readily available. Although the ab-initio study by Voznyy at al. [105] is thorough, their data are not enough for empirical potential parameter fittings. As a result, the purpose of this chapter is to use the ab—initio method with DFT to characterize the bonds between thiols and the GaAs (001) surface, and to construct an energy hypersurface so that a set of empirical potential parameters can be established. 6.2 System and Simulation In this work, the DFT method is employed to study the GaAs (001) surface reconstruction and thiol-GaAs bonds using the QUANTUM-ESPRESSO package [106]. The Vanderbilt ultrasoft pseudopotentials [107] and planewave basis are used, which are common methods used to reduce the computational load of ab—initio calculations. Both local density approximation (LDA) [108] and generalized gradient approximation (GGA) with either Perdew-Burke-Ernzerhof (PBE) [109] or 163 PW91 [110] exchange-conelation functionals are tested. It has been reported that both PW91 and PBE generally overestimate lattice constants relative to experimental values [111], while LDA underestimates them [112]. It is also reported that PBE is better in description of polar surfaces, which is the case for the GaAs (001) surface. Another advantage of GGA exchange-correlation functionals over the LDA is that they usually predict the binding energy with higher accuracy. A planewave cutoff energy of 60Ry is used for all geometry optimization, and it is increased to 100Ry in the binding energy calculation to converge the binding energy better than 0.01 eV. 2x2x2 Monkhorst-Pack [113] K-points grids are used to converge the energy calculation better than 0.001 eV. Spin polarization is used in all binding energy calculations. In the simulations, an 8—layer GaAs substrate is first optimized (see Figure 6.1). PBCs are used in all three directions. The simulation box has a dimension of 35 .3. in the z-direction, leaving a vacuum gap of about 15 A between the topmost and the bottommost layers to avoid any interaction between the system and its image in the z-direction. The dimensions in the x- and the y-directions are chosen to be 2 times the lattice constant. The 8-layer substrate is chosen to make sure that the two surfaces in the z-directions are decoupled with each other. In the optimization, the bottom two layers are fixed, and the rest of the atoms are relaxed. 164 Figure 6.1 An 8-layer GaAs system (green-Ga; white-As) After the optimization, the bottom 4-layer atoms are removed, leaving the top 4-layer atoms. The bottommost Ga atoms are then saturated with H atoms. The hydrogen saturation is necessary because the polar nature of the bottom surface will influence the behavior of the top surface. Then the whole system is relaxed again (see Figure 6.2) except that the bottommost Ga atoms are fixed. W .1. Figure 6.2 GaAs substrate with H terminated on the bottom (green-Ga; large white-As; small white-H). After attaining the optimized GaAs substrate, thiol molecules are placed on top of the substrate. The whole system is then relaxed until its geometry is 165 optimized (see Figure 6.3). Figure 6.3 Thiol on GaAs (001) substrate (green-Ga; large white-As; blue-C; yellow-C; small white-H). 6.3 Results and Discussion 6.3.1 GaAs Bulk Property Reproduction To test different exchange-correlation functionals, bulk properties, including lattice constant, bulk modulus, and specific heat are calculated. These results are presented in Table 6.1. By comparing the results calculated using PBE, PW91 and LDA, it can be seen that both PBE (5.740 A) and PW91 (5.736 A) overestimate the lattice constant, while LDA (5.548 A) underesfimates the lattice constant compared to the experimental value of 5.653 A. These findings coincide with the findings from other references [111, 112] mentioned in the previous section. The bulk modulus predicted by PBE and PW91 is smaller than the experimental value while LDA predicted a bulk modulus very close to the experimental datum. All three functionals produce almost the same specific heat values which are very close to the experimental value. It is seen that all the functionals have their own limitations. Since PW91 is a more recent version of GGA, and it produces almost the same results as PBE, PBE is dropped and PW91 is uwd to represent the GGA exchange-correlation functional. Since the GGA is better in simulating polar surfaces and polymer structures [105], PW91 is used in all binding energy related calculations. PBE PW91 LDA Reference Woe Com (3) 5.740 5.736 5.548 5.653 [114] Bulk Modulus (GPe) 59 61 78 74.8 [115] Specific Heat (J /(g-K)) 0.317 0.316 0.320 0.33 [116] Table 6.1 Calculated bulk properties using different exchange-correlation functionals. 6.3.2 GaAs(001) Surface Reconfiguration To study the thiol adsorption on the GaAs(001) surface, it is important that the simulation produces the correct surface configuration. The surface could be terminated by either As atoms or Ga atoms, or the combination of these two types of atoms. The surface reconfiguration also depends on the initial states of the surface in the simulations. For example, the As-terminated surface reconfiguration can be either the case in Figure 6.4 (a) or (b). For convenience, the configuration in Figure 6.4 (a) is designated as the type I surface, and that in (b) is designated as the type II suface. The surface mometry optimizations are 167 carried out on the 8—layer systems (see Figure 6.1). (a) (b) Figure 6.4 As-terminated GaAs(001) surfaces: (a) surface type I; (b) surface type II (white-As; green-Ga. Note: only the atoms in the top two layers are visualized for clearer view) Both As-terminated and Ga-terminated GaAs(001) surfaces are studied with initial configurations close to both type I and type II surfaces. The final total energies of the optimized surfaces are listed in Table 6.2. It is found that the As-terminated surfaces are much more stable than the Ga—ten'ninated surfaces (see Table 6.3). It has been shown that standard etching procedures used for thiol deposition on GaAs result in an excess of As on surface, and XPS also shows that sulfur binds to As rather than to Ga [117,118]. Our calculated energetics reflect the physics of these experimental findings. It is found alter the geomefi'y optimization that the Ga-ten'ninated has half of the Ga surface atoms move inwards to the surface while the other half pop out. The As-terminated surface has uniform As surfaces. It is also seen from Table 6.2 that the two types of As-terminated surfaces also has a relatively large energy gap (0.84 eV), making the transition from the type I to the type II surfaces difficult. Since the 168 Ga-ten'ninated surface is much less possible in realty, the As-terminated GaAs(001) surface is used in the rest of the study. Stu-face type Configuration energy (eV) As-I 0 As-II 0.84 Ga-I 2.65 68-11 1.76 Table 6.2 Energetics of Reconfigued GaAs Surfaces (note: the energy of the most stable surface type is set to OeV) The final geometries of the As-terminated GaAs(001) surfaces are presented in Figure 6.4 with only the topmost two layers shown for better viewing. The calculated characteristic dimensions of the surfaces are listed in Table 6.3. It is seen that PW91 and LDA produces comparable values with data from PW91 uniformly larger than LDA data. This is believed that PW91 tends to overestimate lattice dimensions while LDA tends to underestimate them, as discussed in Section 6.3.1. For surface type I, the distances between As atoms are, (123 = 4.0911, d24 = 8.10A from pw91 and d23 = 3.93A, at24 = 7.83131 0 from LDA. These values coincide well with the values (6123 = 4A, 6124 = 8A) 0 from ref. [105]. The values of 6124 z 3A are close to the SAM lattice constant (3 .3.) found in ref. [62], which demonstrates that the As-terminated type I surfaces 169 are achieved in the experiments. As a result, we use such a surface configuration for the study of thiol adsorption on GaAs(001) surfaces. sin-race Structure 1 dlz d 23 5124 PW91 (.31) 2-35 4.09 8.10 LDA ( R) 2.48 3.93 7.83 Tcrsofi'(A) 2-42 3.99 7.97 Surface Structure 11 d” d23 d24 PW91 (3.) 2-57 4.35 8.12 1.1511(3) 2.57 4.18 7.32 Tel'sofl'(;\) 2.45 4.25 8.08 Table 6.3 Characteristic Dimensions ofAs-terminated GaAs(001) Surfaces 6.3.3 Thiol Adsorption on GaAs(001) Surface After attaining the optimized GaAs(001) surface, thiols with different chain lengths are placed on top of the bare surfaces, and the whole system is relaxed to find the most stable configurations. The binding energy between the thiol molecules and the GaAs surface is calculated according to eq. (6.1) E binding = GaAs + Ethiol — EthioI—GaAs (6.1) Where EGaAS is the total energy of the relaxed GaAs(001) substrate, Erma] 170 is the total energy of the relaxed thiol molecules, and Erma—GaAs is the total energy of the relaxed thiol-on-GaAs system. For the optimized S-CH3-on-GaAs (see Figure 6.5), the S atom locates right on top of the As to which it bonds. The molecule backbone lies in the plane which is determined by the As dimmers and the 8 atom (call it "dimmer plane” for convenience). The As-S bond length is 2.2543, and the calculated binding energy is 2.28 eV. When S-CH2—CH3 is optimized on the GaAs(001) surface, the S atom is also on top of the bonded As atom, and the molecules backbone still lies in the dimmer plane. The As—S bond length is 2.251131, which is only 0.13% shorter than the As-S bond of the S—CH3—on—GaAs system. The binding energy 2.32 eV is higher compared to the S-CH3-on-GaAs system. Such findings are similar to the findings by Cometto et. al. [119], who found that S—CH2-CH3—on—metal systems had slightly higher binding energies and slightly shorter S-metal bonds than those of the S-CH3-on-metal systems. 6 O. ‘ .- n". y (J {/C/ Figure 6.5 S-CH3 adsorption on As-terminated GaAs(001) surface type I 171 Beyond the S-CH2-CH3 thiol molecule, molecules with longer chains are optimized on the GaAs(001) surfaces. Different from the location of S atoms in the previous two systems, S atoms of the S-(CH2)2-CH3, S-(CH2)3-CH3, S-(CH2)4-CH3 molecules are not located on top of the bonded As atoms (see Figure 6.6 (b)). The backbones also do not lie in the dimmer plane. Franzen [120] also found that the short-chain thiols (with 1 or 2 carbon segments) had different equilibrium geometries and had different binding energies from longer chain thiols. (a) (b) Figure 6.6 S-(CH2)2-CH3 adsorption on As-terminated GaAs(001) surface type I. (a) side view and (b) top view For better comparison, the S-As bond lengths and binding energies are listed in Table 6.5. It is seen that both the bond lengths and binding energies have very little changes when the thiol molecules have 3 or more carbon segments. The S-As bond length are shorter than the S-Au (2.531) [119,121] and S-Cu (2.353.) [119,122]. The shorter bonds also suggest stronger bond strengths and this is proven by the binding energies calculated in this study (S-As: 2.34eV; S-Au: 1.64-2.3eV [119-121,123-125]; S-Cu: 2.19eV [119,122]). It is also shown in Table 6.4 that thiols with longer chains will have shorter S-As bonds and higher binding energies. Thiol molecule As—S bond length ( 12‘) Binding enagy (eV) S-CHB 2.254 2.28 S-CHZ-CHB 2.251 2.32 S—(CH2)2-CH3 2.238 2.34 S—(CH2)3-CH3 2.239 2.35 S-(CH2)4—CH3 2.238 2.33 Table 6.4 bond lengths and binding energies of thiol-on-GaAs To identify which atoms are involved in the thiol-GaAs bond, the electron isodensity surface corresponding to 0.02 a.u. is visualized in Figure 6.7. It is seen that at such a low density of electrons, only the S atom and the bonded As atom are sharing electrons. This phenomenon suggests that the thiol-GaAs bonding is formed mainly by the 8 head and the bonded As atoms, while other atoms do not contribute much to the bond formation. This argument is proven by the binding energy data in Table 6.4 which stays almost constant when the chain length increases. It is possible that the first carbon segment may affect the thiol-GaAs bond because S-CH3 has a slight lower binding energy than longer chains (note that the first carbon segment of S-CH3 is CH3 while the first carbon segments of longer chains are CH2). I73 However, such a difference is so tiny that we can consider the S-As bond as the thiol-GaAs bond. Figure 6.7 Electron Isodensity Surface corresponding to 0.02 a.u. of a S-(CH2)7-CH3 bonded on the GaAs(001) surface type I. 6.3.4 Empirical Potential Fitting As discussed in the previous section, the thiol-GaAs bond interaction is governed by the S-As bond, and characterization of the thiol-GaAs binding energy is equivalent to the characterization of the binding energy of the S-As bond. The thiol—GaAs binding energy as a function of the S—As bond length is calculated using the S-(CH2)3-CH3 thiol molecule. The change of S-As bond length is achieved by moving the thiol molecule in the + z-direction in a 0.2121 interval. The calculated binding energies are plotted in Figure 6.8 as dots. 174 UM (eV) 7/ d (an gstrom) Figure 6.8 S-As bond energy as a function of bond length (dots: data from DFT calculation; solid line: fitted Morse functional) In classical MD simulations, the Morse potential [126] is often used to describe covalent bonds. Since the thiol molecule and GaAs surface are connected by the S-As covalent bond, the Morse potential is chosen to describe the S—As bond. The functional form of the Morse potential is listed in Table 6.6 as eq. (a). The equilibrium binding energy De and the equilibrium bond length rm is prefixed to 2.35eV and 2.2393 as calculated from DFT (see Table 6.5). The parameter a is fitted using the least-square method so that the Morse functional reproduces the one-dimensional energy hyperspace (see Figure 6.8). The final fitted parameters are listed in Table 6.5. It should be noted that the Morse potential with the current set of parameter does not produce very good binding energy when the bond length deviates away the equilibrium value greater than 1 3. However, in reality, the chance of the As-S bond to become 1 3 longer than the equilibrium value is very small. As a result, the current Morse potential should be good enough for simulations at moderate temperatures. The quality of this 175 potential is verified in the simulation of SAM-GaAs system in next chapter. Potential Function forms and parameters Mm UM = De [(1 — eff—i») )2 — 1 .0] (a) As—S where 0 —1 De = 2.356V,a =1.62724A , rm = 2.239 3 Table 6.5 Functional and parameters for As-S Morse Potential 6.4 Summary In this chapter, the ab-initio method employing DFT is used to calculate the energetics of thiol molecules on the GaAs(001) surfaces. Different configurations of GaAs(001) surfaces, including the As-terminated and the Ga-terminated type I and II surfaces, are optimized. The calculated energies are compared, and the As-terminated type I surface is found to be the most stable configuration. Such a configuration is used for all the thiol adsorption simulations. Thiol molecules with different chain lengths are adsorbed on the GaAs surface. It is found that thiols with 1 and 2 carbon segments have similar optimized geometries, but thiols with longer chains have different geometries. Except thiols with 1 and 2 carbon segments, all adsorbed thiols have almost the same S-As bond lengths. Thiols with 2 or more carbon segments have essentially 176 firesamebindingenergywhenbmdedtoflreGaAfiOOflswfaceThidwifllufly 1 carbon has a slightly lower binding energy. A0.2 a.u. electron isodensity suface is visualized, and it is found that the thiol-GaAs bond is downated by the S—As bond. To simulate the S-As bond using an empirical potential, an energy hypersurface is constructed for the S-(CH2)3—CH3-on—GaAs system using ab-initio calculations. A Morse potential is then fitted to the energy hypersurface, and a set of parameters is attained. 177 CHAPTER 7 NEMD STUDY OF THERMAL ENERGY TRANSPORT IN GAAS-SAM-GAAS JUNCTIONS 7.1 Introduction Besides the SAM-on-metal, SAM-on-GaAs is another type of emerging molecular electronics. Different from the thiol-Au interface, there are no empirical potentials directly available for thiol-GaAs bonds. However, an accurate empirical potential for the thiol-GaAs interaction is necessary for an MD to be carried out on the system so as to study the thermal energy transport across the interface. Such a potential has been successfully determined in the previous chapter by fitting a Morse potential functional to an energy hypersurface created in DFT calculations. In this chapter, this fitted Morse potential together with other potentials is employed in the NEMD to study the thermal energy transport in the GaAs-SAM-GaAs junctions. 7.2 Results and Discussions 7.2.1 Potential for GaAs Substrate The Stillinger-Weber (SW) potential [131] was developed to simulate bulk crystalline semiconductor silicon. The SW potential parameters have been expanded to simulate all the Ill-V compound semiconductors [132]. Moreover, SW only consists of two-body and three—body potentials, which makes its 178 implementation easy, and makes the simulation fast. The SW potential is first tested to calculate the bulk properties of GaAs crystal. The calculated bulk properties are listed in Table 7.1. It can be seen that the SW potential reproduce GaAs bulk properties with high accuracy. sw potential Tersoff Reference Lame Comm, (3,) 5.582 5.636 5.653 [114] Bulk Modulus (GPa) 77.3 75.5 74.8 [115] Shear Modulus (GPa) 40.7 48.3 32.9 [133] Young’s Modulus (GPa) 85.6 86.6 85.9 [133] Poisson Ratio 0.32 0.31 0.31 [133] Table 7.1 Bulk properties predicted using empirical potentials However, it is found that the SW potential is not able to reproduce the dimmer As surface as shown in Figure 6.4(a). Such a drawback of SW potential is due to its nature of not considering the bonding environment at the material surfaces. Alternatively, a more sophisticated Tersoff potential [134,135] taking account of the bonding environment by including the bond orders is used to simulate the GaAs(001) substrate. The Tersoff potential functionals can be found in ref. [136] (please note that the ,6 in eq. 6 of ref. [136] should be 7) and parameters are listed as follow: 179 , c = 1.226302 , d = 0.790396 h = —0.518489 n=6.317410, 0:1.560903, 7:0.357192, D, =2.18, r,=2.34 , S=1.641366 , R=3.5 , D=0.1 .16, = 0.2982 , 2,, = 0.3458 The parameter mixing rule of Tersoff [134] was used. All the parameters are taken , from ref. [136], except that ,l. was relaxed so that the dimmerized As surface is more stable at high temperatures as suggested by Kuronen et. al. [155]. The calculated bulk properties using the Tersoff potential is also presented in Table 7.1. It is seen that the Tersoff potential produces bulk properties with equivalent accuracy to SW potentials. The Tersoff potential is also used to optimize the As-terminated GaAs(001) surface type I and II. The optimized dimensions are listed in Table 6.4, and it is seen that the Tersoff potential simulate the GaAs(001) with accuracy comparable to DFT calculations. An optimized As-terminated type I surface is presented in Figure 7.1. Figure 7.1 Optimized As-terminated type I GaAs(001) surface (white-As; green-Ga) 7.2.2 Thiol Adsorption on GaAs(001) Surface Up to this point, all empirical potentials necessary for the thiol-on-GaAs system are obtained. The Tersoff potential are used to simulate the GaAs(001) substrate. The Hautman-Klein potentials [65] (see Table 4.1) are used to model the interaction for thiol molecules. The thiol-GaAs interaction is simulated by a Morse potential as obtained in chapter 6 (see Table 6.6). The vdW interaction is simulated by the LJ potential. The LJ parameters are taken from ref. [65] and [80], and they are listed in Table 7.2. The Lorentz-Barthelot mixing rule [79] is used to combine LJ parameters for unlike atoms. Lennard—Jones UL_J = 45119.)” — (2)6] (a) Wlth where 6' and 0' are determined by mixing rule [79] Lorentz-Berthelot , , , for intramolecular mteractlons C2: mixing rule a = 0.00513eV,a = 3.9143 for intermolecular interactions C2 : a = 0.0046e V,0' = 3.43093 [80] for C3: 5 = 0.0076eV,0' = 3.9143 for Ga: 3 = 0.0180eV,0' = 3.9053 [80] forAs: a = 0.0134eV,0' = 3.7693 [80] for S that interact with other atoms other than S : a =0.01086eV,c=3.5503 Table 7.2 LJ parameters for different atoms 181 To test the above set of potentials, SAM consisting of S-(CH2)8-CH3 molecules are placed on the GaAs(001) surface. The potentials are successful in predicting the geometry of thiols on GaAs(001) surface (see Figue 7.2). The tilt angels to the surface normal of the thiols are calculated by taking the running averages of an MD run, and the calculated value is 60°;t4°. This value is comparable to the experimental findings of 57° :t3° [104]. It should be noted that the core size of C2 segment is reduced (as reflected in Table 7.2) to stabilize the thiol array for the intermolecular interactions. Figure 7.2 Thiols adsorbed on a GaAs(001) surface 7.2.3 NEMD Study of SAM-GaAs Interfacial Thermal Transport 7.2.3.1 Program Test Since the Tersoff potential is complicated and need relatively large computational resources, it is not practical to implement it to our own code which is a serial program. To solve this problem, the NEMD method is implemented to the GULP [130] code. To test the program, a crystalline Argon system is simulated using the modified GULP program. The thermal conductivity of the solid Argon is calculated according to the Fourier’s law, eq. (7.1) J k=— 7.1 . VT ( ) Since it is very difficult to guarantee a perfectly linear temperature profile, a linear fit is used to obtain the temperature gradient (see the black line in Figure 7.3). The calculated thermal conductivity is 0.43W/m-K'. The experimental thermal conductivity for crystalline argon ranges from 0.14 to 2.2 W/m~K [137]. Considering the large temperature gradient, our calculated datum is comparable to the experimental results, which shows the validity of the implementation of NEMD in GULP program. 183 Figure 7.3 System and temperature profile of NEMD on crystalline Argon 184 Temperature (K) 140 120 100 80 20 ceeecececeeemecaeseecmeceeoeceeeeeseecceesec «30‘s- Mac's "03";- ‘wli- it“; V. to ‘f‘d ‘5 ~ Kr"; ‘c‘r'r‘v ‘— {I'- il‘ ‘.-"v 3'1’” Flt-“7"". Okabbe-tao “Beyer-guts, u-a-I'huvaum hfla‘av'u owl—wid- c'u'uer'v'or'd‘dv var tarnish-v3 ‘4 v9 .Uh‘a'abhwh ““801“th 54"! v‘rU‘o‘v‘iU ‘IV‘I‘ VN'U'Z- U‘vl‘ifikd'v‘b’vuk'fi-MSpV'l’ e‘g 9!,7QVV'KfM 9",: I"... infigbca'iv 'atxa.c"ac(n'lv'cc \“OPR1’5"\! 95 V“) Q‘,LI 5"“ r...“ b‘cw‘r'fi {ale 0 .huhuuswuuNVszc995w»;vvhovahULuMCsbseeacauyee9 .‘v H's. ”fir 3v"e‘n 15.04‘ L,-‘ -L;\—r"~c‘ic‘-‘.:cc’t.r'u ink l'u‘a ‘a‘d‘k‘biurvs 9‘: ifi‘ehQ-Uh‘v 3:593ra‘iE‘J” ‘Hc rah 'a’a‘ur‘Jk- ‘c'fpr'or'iufir 'vhr Back 30".) z‘ ‘v"-~‘F'\-'v"\v lv‘- ‘0‘— M‘e rcwis'bSH’V id‘s!" 3‘": '5'” ‘vutl’yffiyvyyhv‘ubuhvw095» use 9% L-‘v'fi UV‘rUi'ic-‘vls'mhs/‘lnh v 'vi-I‘d‘i-bk'o GfiuwuaubeawuSeavyvuevebwwenuuNCUuwaouygvcube“he“ .‘c'B'vV‘v‘r I’M vu‘w Internet 'efih‘cfii '- 1.: ".r‘u L-‘e “u Iii-“9".) 'o-‘e' v “1‘- ‘2"!“51-5'4"; It“: "uh-K10 ..UOIOOflOIQICODOUVUVUOVUOGOOOVQOBQGOGQQQOUfiOEEU. 1 z-coordinates 185 7.2.3.2 NEMD Calculation of SAN-GaAs Interfacial Thermal Conductance To calculate the thermal conductance at the SAM—GaAs interfaces, NEMD simulation is performed on a GaAs-SAM—GaAs junction with SAM consisting of octanedithiols (—S -(CH2 )3 —S - ). PBCs are used in the x- and y-directions, and the free boundary condition is used in the z-direction. Both the heat flux method and the constant temperature method discussed in Chapter 5 are used to create the nonequilibrium state. The left four layers of atoms are chosen as the heat source region, and the right four layers of atoms are chosen as the heat sink region (see the upper portion of Figure 7.4). The velocity scaling interval is 150fs, and the steady state temperature profile is averaged over 10ps. 186 Figure 7.4 System and temperature profile of NEMD on a GaAs-SAM-GaAs junction in (a) heat flux method (b) constant temperature method 187 a) ThmmnfimeCK) 500 4m 300 200 .. .. ..L “3;" F3. 33:3“ 3.3.3.3-. ‘ . r-3-i-3.3:l 1-3.3-3 a .3: .W: . 3_L3-3.3.k.1_0_3.k.3-6.. LILLLJ -anfi- LL3-L1 teas—6-6-0.:LJ- 5.3 ..3 3-3.43-3-3- 0_6«0._ hi—l—a-L—}—t- e.— p 1 *«i—I-t—t—t— ‘ 1 t_iLt_l-l- ..1-l_;Ej"j- ~i-l—e—l-Li»-el» . --.LthathL CFO.— l1111l11 fTI'TYll'TrlllIIII'IIII'. ’___,.. pf, 11.111111111111111: IItl'. zaomdhmws Figure 7.4 cont. (b) Temperature (K) 500 300 200 .aO-I ’ < Source Sink IIII IIII IIII IIII IIrTIIlI IT‘I -I I j I I T I { 1AlllJllILILAJJIILilllJilIl|lAIJ III'IIII'IIIl'llll'TYlI'IIIl'lflI'IIl .1111 hlllllllllIlllJJlllllllllLllllllllLl z-coordinates 189 The thermal conductance is calculated according to eq. (5.2). The calculated value of the left SAM-GaAs interface from the heat flux method is 4691:150MW/m2°K (460K) and 290i95MW/m2-K (340K) for the right interface. The values from the consmnt temperature method is 51 lMW / m2 °K (390K) and 218W /m2 -K (280K) respectively. Following the thermal conductance values are the mean temperatures of the interfaces in brackets. It should also be noted that the perwntage errors are larger than 1 in the constant temperature method as discussed in Section 5.2.1. It isseenthatthethermalconductancesoftheGaAs-SAMinterfaceshave values in the same order as that of the Au-SAM interfaces. The possible reason is that SAM-GaAs and SAM-Au bonds have very similar strengths. Another observation from these data is that the thermal conductance, within the same method, has larger values when the interface is at a higher temperature. This is consistent with the trend of thermal conductance as a function of temperature studied in Section 5.4. The reason should be the same as discussed in Section 5.4 which says that at the higher the temperatures, more phonons are involved in the thermal transport and thus higher thermal conductance as observed. 7.2.3.3 VDOS and Analysis To analyze the then'nal energy transport across the GaAs-SAM interfaces from a vibration coupling point of view, VDOS are calculated according to eq. (5.4) for atoms at different z-coordinates from an equilibrium run at 300K The VDOS of 190 selected groups of atoms are presented in Figure 7.5. The ldreling scheme of atoms is the same as discussed in Figure 5.5. ltisseenflwatduetothebondenvimnmentchangeJhevibrationspectraof the surface GaAs atoms are shifted to lower frequencies. This is because the surface makes the effective bond strengths of the surface GaAs atoms weaker, which leads to lower natural frequencies. It is also seen that C2 and C7, C4 and CS have almost identical VDOS profiles. This is because, due to geometrical symmetry, C2 and C7, C4 and CS have exactly the same bond environment. For more complete discussion of the VDOS properties, the vibrational modes are again divided into three categories (LF, IF and HF) as did in Section 5.2.2. It is found that there is no HF mode as those found in ref. [100]. The HF modes in ref. [100] are the results of the fast vibrating hydrogen atoms. The absence of the HF modes is because the united atom model used in the present work does not consider the high frequency vibration of hyd'ogen atoms. It is also found in ref. [100] that the HF modes are highly localized, which does not contribute much to the thermal energy transport. As a result, the united atom model is an appropriate model for the thermal energy transport simulation in this study. It is seen that the LF modes are directly coupled to the GaAs phonons, and these modes are believed to be the main carrier for thermal energy transport across me GaAs-SAM interfaces. GaAs substrate have vibration mode with frequencies up to about 12THz. Compared to Au stbstrates, which have frequency cutoff of about 8THz 191 (see Figure 5.5), the GaAs substrate provides more modes for coupling to the SAM vibration mode. This leads to higher thermal conductance of GaAs-SAM than Au-SAM interfaces. 192 Figure 7.5 VDOS of atoms at different z-coordinates (GaAs-in: GaAs inside the substrate; GaAs-sun GaAs at the surface; S-Sulfur, C-Carbon with the same labeling scheme in Figure 5.5). 193 AN 1.5 3:258...— _. 4. wvvvv'fvvvv’Vvvvvwvvvv—‘fi— ILA-AAIAAAIIIIAA‘Ilfilfl 8 N0 .Fbk'b’b’lhl’lbhtb Pppbht ’P...’..p I.’lp'b|.lbl’ "‘I‘ ‘ 1“ 1“ llAlAl-AAIIIIIAIALAL-l‘ wwvvvaivvvvvvvvvv'v . 3 m aria 'bP-I’I..’.Ib’Ipb.llr£l’b..h.pl.prlr. h.F..L.l."bhl’b.h>.’.lb.h”...bhl’br}b hblpp”’Ibl’lbb-P’Ph’....P..Db'l”.'hl ’l’ ’I ’tl’b’tb A-AAIA-IA-AA-IIJLIMAAAA 194 I‘ll-IIIAIAIIIAALILAL-A The thermal energy transport inside the SAM molecule should be the same as that discussed in Section 5.2.2 for the Au-SAM-Au junctions: both LF modes and IF modes are delocalized, and the thennal transport is almost harmonic which is very efficient. To better quantify the degrees of localization, the participation N ratio [159] part is calculated, which describes the number of atoms participated in a certain vibration normal mode with frequency a) . As in ref. [159], the energy moment is defined as Mp ((0) = Z [8i(w)]p 7.2 atom i ( ) where 3,.(w) is the kinetic energy of atom 1' associated with frequency a). As discussed in Section 5.2.7, the mass weighted VDOS is directly proportional to the kinetic energy. As a result, the magnitude of the weighted VDOS is used as 51(0) . The participation ratio is calculated as Np(w)=M12 /M2 (7.3) By analyzing eq. (7.3), it can be easily shown that if the mode with frequency a) is a traveling wave (fully delocalized), N p (4)) equals the total number of atoms. If the mode is purely localized, N p (4)) equals 1. To investigate the level of modal delocalization among different groups of atoms, three groups of atoms are sampled to calculate the participation ratio. In Figure 7.6 (a), the participation ratio of a total of 14 atoms from a SAM molecule 195 and the first two layers at the GaAs(001) surface are calculated (two Ga atoms, two As atoms, two S atoms and eight C atoms). In Figure 7.6 (b), a total of 10 atoms from a SAM molecule are sampled for the participation ratio calculation (two S atoms and eight C atoms). In Figure 7.6 (c), the eight C atoms in a SAM molecule are sampled to calculate the participation ratio. To better compare the effect of different atoms on the participation ratios, the ratio is normalized against the total number of atoms sampled. As a result, the fully delocalized mode would have a participation ratio of 1. It is seen in Figure 7.6(c) that in pure C atoms, the modes with frequency from O to 8THz and from 23 to 27THz have participation ratios near 1, suggesting highly delocalized features. This makes both LF and IF modes have a highly delocalized portion. These delocalized modes can transport thermal energy very efficiently among these C atoms. As a result, the temperature gradient in the C segment is very small (see Figure 7.4). Due to presence of S atoms, which has only a few lF modes, the participation ratio of the SAM molecule atoms is reduced for the IF modes (see Figure 7.6(b)). Although the participation ratio of modes with frequency from 6 to 1OTHz decreased due to the spectra mismatch between S and C atoms, the LF modes are still greatly delocalized. As a result, the temperature difference between S and C atoms are still small (see Figure 7.4), although not as small as that among C atoms. 196 Figure 7.6 Participation ratio of (a) GaAs-SAM-GaAs junction (14 atoms), (b) SAM molecule (10 atoms) and (c) pure carbon segments (8 atoms). 197 =1...— Secs—.2... or o o Vf O l n Q C) 6.: 1 orieu uomfprima neareuuou 1- - 198 However, when the surface atoms of the substrates are considered, the participation ratio of the IF modes are greatly mduced due to the lack of IF modes in the substrate (see Figure 7.6(a)). One can image that if a large number of substrate atoms are sampled, the IF modes would eventually have participation ratiosofO.TheparticipationratiosofLFmodesarealsoreduced,butthese modes still have enough delocalized feature to facilitate efficient thermal transport between the substrate and the SAM molecules. 7.3 Summary In this chapter, a set of ernpiricel potentials is established to simulate the GaAs-SAM-GaAs junctions. This set of potentials is tested by optimizing thiol geometry on a As-terminated GaAs(001) surface, and is demonstrated to produce good results compared to experimental measurements. The GULP program is modified to implement the NEMD method, and the program is tested to predict thermal conductivity values of solid Argon which are comparable to experimental data. GaAs-SAM interfacial thermal conductances are calculated at moderate temperatures using both the heat flux method and the constant temperature method. It is found that these values are uniformly larger than the Au—SAM interfaces. VDOS’s of atoms from different components of the junction are calculated. It is found that only the LF modes of the SAM molecules can directly couple to the GaAs phonons which only have modes with frequencies lower than 12 THz. Compared to Au VDOS which has a cutoff frequency of about 8THz, 199 GaAshasmaepoptflatedLFmodeswhidlresultsinbettervibrationcmrpling between the SAM molecules and the substrates. By calculating the modal participation ratios among different groups of atoms, the level of delocalization of differentvibrationmodesaremeasured. ltisfoundtl'lattheboththeLFand IF modes are highly delocalized among the C atoms, which lead to efficient thermal energy distribution among the C segments. The LF modes are also delocalized to S atoms, which facilitate thennal energy transport along the whole SAM molecule chain. Due to the spectra mismatch between the substrate and the SAM molecules, the delocalization feature of the modes is geatly mduoed. The thermal transport between the substrate and the SAM is governed by the direct LF modes coupling to the low frequency substrate phonons. 200 CHAPTER 8 NEMD STUDY OF THERMAL ENERGY TRANPORT ACROSS GRAPHENE-POLYMER INTERFACES 8.1 Introduction Polymer nanocomposite materials have been manufactured to enhance the mechanical, thermal and electrical properties of pure polymers [138]. These materials have potential applications in aerospace [139], automotive manufacturing [140], and electronics [141]. Carbon Nanotube (CNT) polymer nanocomposites have already received significant attention from researchers, because CNTs have unusual mechanical and electronic properties [142]. Recently, Stankovich et. al. [143] presented a general approach for the preparation of graphene-polymer composites via complete exfoliation of graphite. The graphene-based composite materials have become more and more popular because single graphene sheets have extraordinary thermal conductivity (~ 3000W/ m-K ), mechanical stiffness (lococpa) [144], electronic transport properties [145-149], and moreover, inexpensive cost (graphite is dollars per pound) compared to CNTs (hundreds of dollars per pound). Since bulk polymer has very low thennal conductivities ~ 0.2W / m . K , graphene sheets intercalated polymers are expected to improve the thermal transport properties. 201 Studies on the graphene-polymer composite have been focusing on the structural [150-152] and electronic [153,154] properties, leaving the thermal transport a relatively unexplored area. Since thermal conductance in bulk polymer and pure graphene sheets are well studied, knowledge of graphene-polymer interfacial thermal transport is critical in understanding the thermal transport in the graphene intercalated polymer composite materials. In this chapter, different graphene-polymer interfaces 3e constructed, and are used in MD simulations. NEMD is applied to these interfaces to study the thermal energy transport The pupose of this chapter is to provide some preliminary results on the themial transport across graphene-polymer interfaces to lay a foundation for future research. 8.2 Simulation and System The Brenner’s REBO potential [129] is chosen to simulate the bonded interactions in the carbon system. Since the Brenner potential does not consider the vdW interactions, LJ potentials are added to the Brenner’s model to Simulate the vdW interactions. LJ potentials are applied to intermolecular interactions, and the intramolecular interactions for atoms apart more than 2 bonds. The parameters for atoms in the polymer molecules are from ref. [80]. The parameters for graphene carbon atoms are from ref. [156]. The Lorentz-BertheIot mixing rule is used to combine the LJ parameters. The LJ parameters together with the functional form are listed in Table 8.1. 202 The polymer Simulated in this study is paeffin with pure C3,,H62 molecules. The predicted density using the chosen potentials at 298K and 1atm is 0.86 glcc, which agrees well with experimental data of 0.87-0.91 glcc [114]. The potentials are also used to optimize the graphene interlayer distance, and the calculated value is 3.428513. which is very close to the experimental measurement of 3.354}. [157]. Lennard-Jones with Lorentz-Berthelot mixing rule UL—J=4E[l%lu-l%ll where 6' and 0' are determined by mixing rule [79] for H : s = 0.0019096eV,a = 2.571113. [80] for C in polymer molecules: s = 0.0046eV,a = 34309.4( [80] for C in graphene: s=0.00239eV,c =3.4121.3. [156] Table 8.1 LJ parameters for polymer and graphene simulations. The generation of the amorphous polymer is critical to the Simulations. The bulk amorphous polymers are constructed using the modified Markov process [26] with bond conformational probabilities chosen to account for both intramolecular and intermolecular non—bonded interactions. 203 8.3 Results and Discussions 8.3-1 Thermal Conductivity of Pure Polymer To further validate the potential model, thermal conductivity is calculated for a pure amorphous polymer system using the NEMD with the heat flux method. The simulation supercell consists of 30 (3301162 molecules with periodic boundary conditions applied in all three spatial directions (see Figure 8.1). The heat source is set at the ends of the supercell, while the heat sink is located at the center. Thermal energy flows in to the ends and flows out from the center. The steady state temperature profile (after 80ps) is visualized in the lower portion of Figure 8.1. A linear fit is used to determine the temperature gradients VT, and the thermal conductivity is calculated according to eq. (7.1). 204 Figure 8.1 Pure amorphous polymer system and the NEMD temperature profile. 205 Source Sink Source 5m _ I l I l l I I l l I I l I 1 '7 l V I I I I r l I T l l l '4 450 LA 3*- : '1 Al“ - _ g, - .- .‘3 A .1 A400 L' ‘11“ ‘1 '2 l- ‘ u 5 - A .x . 0 _ AA A A “350 " An " 3 ~ ‘A A : 2 I t A‘ 2 g“ T ‘AA it: 1 o ' ‘A .‘ ‘ r- : A A A 3 25° ' ‘2 x“ 1 I A A . : A‘s-1 A : 2m "' fi gait "' ' iii : : {_ . 1w l 1 1 1 1| 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 l 1 1 l 1 - z-coordinate 206 The thermal conductivity is thus calculated to be 0.40:0.09W/m-K. The experimental thermal conductivity of paraffin wax at room temperature is around 0.25W lm-K [160,161]. However, considering the large temperatures imposed on the simulation system and the temperature dependence of polymer thermal conductivities [162], our calculated result is comparable to the experimental data. To better control the temperature range of the system, the constant temperature method is used with source region and sink region being fixed at 350K and 250K respectively. The calculated thermal conductivity is 0.27W/m-K which agrees very well with the experimental value. As a result the present potential model is capable of simulating thermal energy transport in the polymer systems. 8.3.2 Thermal Transport across Polymer-Graphene Interfaces There are different possible topologies of polymer-graphene interfaces since the graphene layers can take any orientation. In the present work, the interfaces formed by paraffin polymers and graphene basal planes are studied. As discussed in Section 8.1, the method has already been developed for complete exfoliation of graphite [143]. This means that single graphene sheets can exist in the polymer-graphite composite materials. As a result, the interfaces formed by a single graphene sheet and amorphous polymer are investigated. The Simulation supercell consists of two graphene sheets and two blocks of amorphous paraffin polymer (see Figure 8.2) with PBCs used in all three spatial directions. In this system, each graphene sheet forms two interfaces with the polymer blocks 207 surrounding it. The heat sink and source are located at the centers of the polymer blocks. 5m _1 I I I 1 I I I I I I I I I I l I I T r '—4 450 L + -3 3 A“ “A ; 400 E- l *i '3 15350:- f l '5 g : l 3 a» — ‘ ~ i300 7‘ 1 E l : 5250 L" l,‘ '1 l- - . - 200 "- AAA“ BR“ —‘ : A ‘A‘ 3 150 E- l ' € l- -l Pl 1 1 l I I l l l I I l l L l I l l l J 1 10° 0 10 2o 30 40 z-coordinate Figure 8.2 Single graphene sheet-polymer interface system and the steady state temperature profile. From the steady state temperature profile, after 80ps (see Figure 8.2), it is seen that the largest temperature changes occur at the graphene—polymer interfaces. The temperature discontinuity at the interfaces suggests large interfacial thermal resistances compared to the polymer blocks. Using eq. (5.2), 208 the calculated mean thermal conductance at the interfaces is 269 i“ 70MW / m2 ° K . Although the only interaction that can transfer thermal energy between graphene sheets and polymers is the weak vdW interaction, the calculated thermal conductance is of the same order as the GaAs-SAM and Au-SAM interfacial thermal conductances. The reason is to be explored in future studies. Although complete exfoliation of graphite is possible, the possibility of incomplete exfoliation still exists [163], leaving stacks of graphene sheets in the polymer composite. As a result, thermal energy transport across the graphite (made of stacks of graphene sheets)-polymer interface is studied. The simulation settp is presented in the upper portion of Figure 8.3, with the source and sink regions locate in the grmhite segments. However, after 50ps of production, the temperatures of different parts still remain closely to the initial mean temperature of 100K which are apart from the sink and souce temperatures. This suggests that the thermal energy transport is impeded by the graphene layers. It demonstrates that the grqihite interlayer interaction (vdW) is too weak to transport thermal energy efficiently. Theoretically, it has been predicted by Sun. et. al [158] that the thermal conductivity along the C-axis is 4 orders of magnitude smaller than in the graphite basal plane. 209 g! n '. 'i'rFr, .7 rFI ... y. . ,3: “$35! U Irv Source Sink Source 300 I. I T r v I r r I r r T I I I 7 j 250 T A T, i l l. “ -l g 200 l" u g . g. 150 - '. GE, : . an.“ A AM ‘ 1 E-l 100 ,_ 5““ Alla Am‘ ‘A‘ 4 : ‘ ‘A I l- A AA fl .1 50 L .. “ .... j l J l aJ l L J l 1 0 20 40 60 z-coordinates Figure 8.3 Graphite-polymer interfaces and temperature profile From these two case studied, it has been shown that the interfaces formed by the graphite basal planes and polymers will mduce the effective thermal conductivity due to the interfacial thermal resistance. In the two cases studied, the interfacial interactions are only weak vdW interactions, which are not good at thermal energy transport. Another possible configuration of the graphite-polymer interface is formed by the edge of the graphene and polymers. In such a topology, the graphite-polymer bond can be formed. The bonding interactions are usually more effective in thermal energy transport as studied in previous chapters. As a result, we believe that the improved thermal conductivity of graphite-polymer 210 composite material is due to the presence of such bonding interfaces. Thermal energy transport across such a graphite edge—polymer system is to be investigated in future studies. 8.4 Summary In this chapter, a set of potentials to simulate the graphite-polymer system is tested. This set of potentials predicts a polymer density (0.869/cc) which agrees well with experimental results (0.87-0.91glcc). It also mproduces good graphite interlayer distance and polymer thermal conductivity. NEMD are performed on two different interfaces formed by graphene basal planes and polymers. It is found that the graphene-polymer interfacial thermal conductance has a moderate magnitude. It is also found that the graphite interlayer flierrnal resistance is very large and greatly impedes the thermal transport in the C-axis. More work needs to be done to explore the mechanism of thermal conductivity improvement of graphite—polymer composite materials. 211 CHAPTER 9 CONCLUSIONS In the present study, thermal energy transport across different material interfaces are studied using different approaches, including ab-initio MD, equilibrium classical MD and non-equilibrium MD. The material interfaces studied in this work include semiconductor-semiconductor interfaces, Au-SAM interfaces, GaAs-SAM interfaces and graphene-polymer interfaces. From the different simulation results, the folan conclusions are reached: (1). Thermal energy transport is more efficient in pure materials than thermal energy transport across material interfaces because the interfaces scatter energy can'ying phonons. (2). lnfonnation contained in the atomic vibration power spectra is useful in explaining the thermal transport phenomena. If the contacting two materials making up the interfaces have large spectra overlap, which indicates strong vibration coupling, thermal energy transport beMen these two materials will be efficient. If there is little or no overlap, the themial communication between the materials will be difficult. (3). CPMD simulation of thermal energy transport gives physically reasonable results showing that it is a promising method of studying nano-scale thermal transport phenomena and calculating thermal transport properties. (4). Thermal conductance calculation of Au-SAM-Au junctions showed that the Au-SAM interface has a much larger thermal resistance than the Au substrates and the SAM molecules. The As-SAM interface is the main barrier for 212 the thermal energy transport across the junctions, which was demonstrated by the junction’s response to a heat pulse. From the vibration coupling point of view, the low VDOS’s of the SAM molecules worked as bottlenecks, and thermal energy transport from Au to SAM molecules influenced by the VDOS mismatch of these two materials. (5). Thermal transport along the SAM molecule chains is nearly ballistic, and thus very efficient. The LP modes are delocalized which transport energy rapidly. The IF modes are also delocalized but it does not have direct coupling to Au substrate phonons, and thus does not contribute much in energy transport across Au-SAM interfaces. (6). The Au-SAM interfacial thermal conductance increases with temperature increase at low temperatures, and it becomes almost constant at high temperatures. The increase of thermal conductance at low temperatures is due to the increased population of thermal transport phonons. As more phonons are excited to participate in thermal energy transport, the transport becomes more efficient. On the other hand, anharrnonic phonon scattering impede thermal transport efficiency across the interfaces. As temperature increases, the anharrnonic scattering becomes stronger, and then further increase of interfacial thermal conductance is suppressed. (7). The Au-SAM interfacial thermal conductance is not affected by the normal pressure exerted on the substrate. (8). Chain length (from 8 carbon segments to 10 carbon segments) of SAM molecules does not have an obvious effect on the thermal energy transport in the Au-SAM-Au junctions. This is due to the fact that thermal energy transport is dominated by ballistic transport along the SAM molecule which is thus not sensitive to the limited length change. 213 (9). SAM molecular coverage on the substrate greatly influences thermal transport across the Au-SAM interfaces. It is found that every molecule works as separate energy transport channels. The Au-SAM interfacial thermal conductance is proportional to the number of available channels. (10). Au-SAM bond strength has an impact on the interfacial thermal transport. Stronger Au-S bonds lead to smaller anharmonicity at the interface, which means less phonon scattering, and thus smaller interfacial resistance. (11). All the calculated Au-SAM interfacial thermal conductance values are between 150 and 450MW/(m2-K), which is inside the range of the experimentally reported thermal conductmce of metal-nonmetal interfaces 8 < G < 700 MW/(mzK) [es-es]. (12). DFT calculation of GaAs(001) surface reconfiguration shows that the As-terminated type I surface is the most stable surface structure. Such a result is consistent with experimental findings [117,118]. (13). DFT calculation of thiol binding to GaAs(001) surface showed that thiols with 1 or 2 carbon segments have different geometries compared to longer chains. Except thiols with 1 and 2 carbon segments, all adsorbed thiols have almost the same S-As bond lengths. Thiols with more than 2 carbon segments have essentially the same binding energy when adsorbed on the GaAs(001) surface. (14). Visualization of a 0.2 a.u. electron isodensity surface for the thiol-on—GaAs system shows that the thiol-GaAs interaction is dominated by the S-As bond. (15). A Morse potential is fitted to the energy hypersurface from DFT calculations, and a set of parameters is attained. This set of potentials is able to 214 reproduce the tilt angle of the thiols adsorbed on an As-terminated type I GaAs(001) surface. (16). Calculated GaAs-SAM interfacial thermal conductances by NEMD range from 218 to SllMW/(m2 'K), which are in the same order of Au-SAM interfacial thennal conductmces. The thermal transport mechanism in the GaAs-SAM-GaAs junctions is believed to be the same as that in the Au-SAM-Au junctions. (17). Preliminary results of NEMD study on thermal energy transport in polymer-graphene systems suggests that special bonding between graphene sheets and polymer molecules are necessary to improve the thermal transport ability of the graphite-polymer composite materials. Besides all the results obtained from this study, a regime map of Tlen [164] is further completed to direct modeling of thermal energy transport (see Figure 9.1). 215 Figure 9.1 Regime Map for Thermal Energy Transport 216 've Ballistic- ' Fourier Heat l Conduction Hyperbolic Heat " conduction as a 1% Classical MD *i’i Ab—initio MD iii a .r .t .. x . x \ .. . s, . .. ° .. . . . .. a . . . . a s. . _.,. _ . . . x x. .. s x x . . .. . a .... 2‘ . . .. ..1 r, .. . . r. I . x. p \t. . .s .1 .. ..\ .\. . . \. .. K . .... I . I . . 3.. x , . . .. \ . . .. . s .. .. n I] .V J!» ‘1'! 14 . . . . X // ./ / ,/ ... r. // / . z . . ., . ... x x. \ . a In. .. I. 4r. 1 .. x .. \ t \ I i L: II. {Ill 1! .... . . 1.0, O O. , 1 . .1 sneer-a, seeeal 217 At short length scales the Boltzmann equation describes energy transport. At short length scales for longer time frames, the Ballistic Diffusive Equation is reasonable, while for short times and large length scales the hyperbolic models work reasonable well. For long times and large length scales, the classic Fourier heat conduction models the diffusive processes very well. At a time scale up to several hundreds of picoseconds and a length scale up to several hundreds of nanometers, classical MD is capable of simulating thermal phenomena. Given a set of reliable potentials, classical MD can give results with betteraccuacythanmeBoltzmannequafionsbecauseitdoesnotrelyonflie empirical parameters, such as mean free path and characteristic timescale, as needed in the Boltzmann equations. At a time scale rm to several picoseconds and a length scale up to several nanometers, the ab-initio MD can be used to simulate thermal energy transport phenomena with even greater accuracy than classical MD. 218 CHAPTER 10 FUTURE DIRECTIONS 10.1 QMMM in Thermal Energy Transport across Material Interfaces In this project, CPMD is proven to be a promising tool to study thermal energy transport in nanoscele systems. However, more work needs to be done to improve the method so that it can be used to study larger systems. A possible approach is to employ the so called QMMM [127] (Quantum mechanics/Molecular mechanics) method to take advantage of the accuracy of quantum mechanics and the speed of the classical molecular dynamics. In simulations of thermal energy transport across material interfaces, the most difficult part to simulate is the structues and bonds at the interfaces. It is often very difficult to find a set of potentials which can simulate the interface behavior accurately. It is possible to utilize the QMMM towards the simulation of interface problems. For the bulk materials that are far from the interfaces, the atoms can be Simulated by classical MD to save computational time compared to ab-initio MD. For the atoms around the interfaces, ab-initio quantum MD can be used to attain higher accuracy thm classical MD. 10.2 Electron-Phonon (e-ph) Interaction in Thermal Energy Transport MD is sometimes combined with Boltzmann equations to study phonon transport [128]. To use the combined approach to study interfacial thermal transport, phonon mean free path and interfacial reflectivity are usually obtained from MD calculations, and they are input to Boltzmann equations so that thermal 219 transport properties can be calculated. The phonon mean free path and interfacial reflectivity are modal dependent, i.e. phonons with different frequencies have different mean free paths and different reflectivity. When studying the phonon transport in an electron active system, the e—ph interaction affects the phonon transport Due to the phonon scattering by the e-ph interaction, phonon mean free paths are changed compared to the cases where no active electrons are present The e-ph interaction strength an be calculated from the ab-initio method through the QUANTUM-ESPRESSO package [106]. If a method can be developed to relate the strength of e-ph interaction to the phonon mean free paths, the e-ph interaction effect on phonon thermal energy transport can be investigated. 10.3 Electron Contribution to Thermal Energy Transport Electron current are account for the electron transport, however, at the mean time, it can transport thermal energy, especially in metals. The electrical thermal conductivity can be calculated according to the Wledemann-Franz Law (eq.(10.1)) [18] kc = N LaeT (10.1) where ke is the electrical thermal conductivity, N1, is the Lorenz number, 0e is electrical conductivity, and T is temperature. The electrical conductivity 0' e of a junction can be calculated using the ‘pwcond’ implemented in the QUANTUM-ESPRESSO package [106]. 220 10.4 A More Sophisticated Potential Model for Au Substrate In chapter 4 and 5, the empirical potential used to simulate the Au substrate is a pair-wise Morse potential. Such a potential is used because it is simple and easy to implement Moreover, it reproduces Au bulk properties with satisfactory accuracy. However, a more sophisticated empirical potential, the embedded atom method (EAM), can be used to simulate the Au substrate. Phonon dispersion can be calculated using both Morse potential and EAM, and the one can compare which potential provides more accurate prediction. A potential with better accuracy in predicting the phonon dispersion relation is expected to produce better thermal conductivity values. 10.5 A Full Anharmonic Potential for SAM Molecules In this project, the potentials used to simulate the SAM molecules are mainly from ref. [68]. This set of potentials is near-hennonic (see Table 4.1), and the hydrogen atoms are integrated into the carbon backbones. Segal et. al. [100] proved that even if a full aiharmonic potential model is used for molecular wires, the ballistic phonon transport inside the molecule chains are not affected compared to the harmonic models. However, using a full anharmonic potential, such as the Brenner’s potential [129], to study the thermal transport across SAM is still worth doing, and it is expected to give more insight into the phonon scattering in SAM molecules. Moreover, in the Brenner’s potential, hydrogen atoms are explicitly modeled, which present a potential scattering source for phonon transport. 221 10. BIBLIOGRAPHY Moore’s Law. 2007. http://www.intelcom lntemational Technology Roadmap for Semiconductors. 2005. hmzllgublicjtrsnet A J. Heeger, 2000, 'Semiconducting and metallic polymers: The fourth generation of polymeric materials," Nobel Lecture. h :lew. l. chemi Ilaureat -Iecture.html Z. Wang, J.ACarter,A Lagutchev,Y.KKoh, N. H. Seong, D.GCahill, and D. D. Dlott 2007. Ultrafast Flash Thermal Conductance of Molecular Chains. Sg'ence, Vol. 317. no. 5839, pp. 787 — 790 A Henry and G. Chen, 2008, High Thermal Conductivity of Single Polyethylene Chains Using Molecular Dynamics Simulations, Phys. Rev. Lett. 101, 235502 P. L. Kapitza, Zh. Eksp. Ther. Fiz. 11, 1 (1941) [English Transl.: J. Phys. USSR. 4, 181(1941)] I. L. Bekarevich and I. M. Khalatnkov, 1961, in Proceedings of the 7‘" International Conference on Low Temperature Physics, G M. Graham and A C. Hollis Hallett, Eds. (University of Toronto Press, Toronto), p. 480 E. T. Swartz, R. O. Pohl, 1989, "Thermal boundary resistance,” Rev. Mod. Phys. 61, 605. D. G Cahill, W. K Ford, K E. Goodson, G D. Mahan, 'A. Majumdar, H. J. Maris, R. Merlin, S. R. Phillpot, 2003, “Nanoscale Thermal Transport,” J. Applied Physics, Vol. 93, pp. 793-818. 0. M. Wilson, X. Hu, D. G Cahill, P. V. Braun, 2002, "COIIOidaI metal particles as probes of nanoscele thermal transport in fluids," Phys. Rev. B 66, 224301. 222 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. PK Schelling, S.R. Phillpot, P. Keblinski, 2002, ‘Phonon Wave-Packet Dynamics at Semiconductor interfaces by Molecular-Dynamics Simulation,” Applied Physics Letters, Vol.80, No.14, 24842486. D. Poulikakos, S. Arcidiacono, S. Maruyarna, ”Molecular Dynamics Simulation in Nanoscale Heat Transfer. A Review,” Microscale Thermophysical Engineering, 7:181-206 (2003). S. Volz, J-B. Saulnier, M. Lallemand, ”Transient Fourier-Law Devices by Molecular Dynamics in Solid Argon,” Phys. Rev. B, Vol. 54, No. 1 (1966). JR Lukes, D.Y. Li, XG Liang, C.L Tlen, ”Molecular Dynamics Study of Solid Thin-Film Thermal Conductivity,” J. Heat Transfer, Vol. 122 (2000). Y. Chen, D. Li, J. Yang, Y. Wu, J.R. Lukes, A Majumdar, ”Molecular Dynamics Study of the Lattice Thermal Conductivity of Kr/Ar Superlattice Nanowires,” Physica B 349, 270-280 (2004). AR. Abramson, CL. Chen, A Majumdar, ”Interface and Strain Effects on the Thermal Conductivity of Heterostructures: A Molecular Dynamics Study,” J. Heat Transfer, Vol. 124, 963-970 (2002). D. Marx and J. Hutter. 2000. Ab initio Molecular Dynamics: Theory and Implementation. Modern Methods and Algorithms of Quantum Chemistry, NIC series, Vol. 1, pp. 301 -449 J. M. Ziman. 1960. Electrons and Phonons. QxLorg University Press, New York R. Car, and M. Parrinello. 1985. Unified Approach for Molecular Dynamics and Density-Functional Theory. Physical Review Letter. 55, no. 22—25: 2471 -2474 G Pastore, E. Smargiassi, and F. Buda. 1991. Theory of ab-r'nitio molecular-dynamics calculations. Phys. Rev. A 44, 6334 S. Nose. 1984. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81, 511 W. G Hoover. 1985. Canonical dynamics: Equilibrium phase-space 223 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. distributions. Phys. Rev. A 31, 1695 W. Kohn and L. J. Sham.1965. Self-Consistent Equations Including Exchange and Correlation Effects. Phfi. Rev. 140, A1133 D. C. Langreth and J. P. Perdew. 1980. Theory of nonuniform electronic systems. I. Analysis of the gradient approximation and a generalization that works. Phys. Rev. B. 21, 5469 J. M. Haile. 1997. Molecular Dynamics Simulation, Wiley, New York M. P. Allen and D. J. Tildesley. 1987. Computer Simulation of Liquids, Oxford University Press, New York D. A Broido and A Ward. 2005. Lattice Thermal Conductivity of Silicon from Empirical lnteratomic Potentials. Phys. Rev. B, 72, 014308 C. J. Gomes, M. Madrid, J. V. Goicochea, and C. H. Amon. 2006. ln-Plane and Out-Of-Plane Thermal Conductivity of Silicon Thin Films Predicted by Molecular Dynamics. Transactions of th_e A§M§, vol. 128, 1114 R. J. Stevens, P. M. Norris, and L. V. Zhigilei. 2004. Molecular-Dynamics Study of Thermal Boundary Resistance: Evidence of Strong Inelastic Scattering Transport Channels, Proceediggs of IME§E04 P. K Schelling, S. R. Phillpot, and P. Keblinski. 2002. Phonon wave—packet dynamics at semiconductor interfaces by molecular-dynamics simulation. Applied Physics Letters, 80, 2484—2486 P. K Schelligg, S. R. Phillpot, and P. Keblinski. 2004. Kapitza conductance and phonon scattering at grain boundaries by simulation. J. Appl. Phys. 95, 6082 Y. Chen, G ng, D. Li, and J. R. Lukes. 2006. Thermal Expansion and Isotopic Composition Effects on Lattice Thermal Conductivities of Crystalline Silicon, Proceedings of IMEQE2006 K V. Tretiakov, S. Scendolo. 2004. Thermal Conductivity of Solid Argon from Molecular Dynamics Simulation, ,1. Chem. Phys. 120, 3765 224 37. 38. 39. 40. 41. 42. H. Kaburaki, J. Li, S. Yip. 1998. Thermal conductivity of solid argon by classical molecular dynamics. Materials Research may Smggium Proceeding, 538: 503 P. K Schelling, S. R. Phillpot, and P. Keblinski. 2002. Comparison of atomic-level simulation methods for computing then'nal conductivity. Physical Review B, 65, 144306 M. Matsumoto, H. Wakabayashi and T. Makino. 2005. Thermal Resistance of Crystal Interface: Molecular Dynamics Simulation. Heat Transfer-Asian Rfll’d‘l, 34, 3 Ryogo Kubo, Mario Yokota and Sadao Nakajima. 1957. Statistical-Mechanical Theory of Irreversible Processes. ll. Response to Thermal Disturbance, J. Phys. &. 4m. 12, 1203 S. G Volz, G Chen. 2000. Molecular-dynamics simulation of thermal conductivity of silicon crystals. Phyg'ggl Review B, 61, 2651 A J. H. McGaughey, M. Kaviany. 2004. Thermal conductivity decomposition and analysis using molecular dynamics simulations, Part II. Complex silica structures. IW' ,Mrnal of Heat and Mass Transfer. 47, 1799 CPMD, http://www.cpmd.orgl, Copyright IBM Corp 1990-2008, COpyright MPI fClr FestkOrperforschung Stuttgart 1997-2001. Humphrey, W., Dalke, A and Schulten, K, ”VMD — Visual Molecular Dynamics”. W 1996. vol. 14. pp. 33-3 N. Troullier and J. L. Martins. 1991. Martins-Trouiller. Phys. Rev. B, 43, 1993. JP. Perdew and A Zunger. 1981. Local Density Approximation. Ph . Rev. B 23, 5048 A. Dodabalapur, L. Torsi, and H. E. Katz, ”Organic Transistors: Two-Dimensional Transport and Improved Electrical Characteristics,” Science 268, 270 (1995). 225 45. 47. 48. 49. 50. 51. 52. 53. 55. F. Gamier, R. Hajlaoui, A Yassar, and P. Srivastava, ”All-Polymer F ield-Effect Transistors Realized by Printing Techniques,” Science 265, 1684 (1994). H. Sirringhaus, N. Tessler, and R. H. Friend, ”Integrated Optoelectronic Devices Based On Conjugated Polymers,” Science 280, 1741 (1998). K Kurabayashi, 2001, Anisotropic Then'nal Properties of Solid Polymers, lntemational Journal of Thermophysics, Vol. 22, No. 1, pp 277-288. R. F. Service, ”Molecular Electronics: Next Generation Technology Hits an Early Mid-life Crisis,” Science 302, 556 (2003). P. Lugli, A Pecchia, M. Gheorghe, L. Latessa, A DiCarlo, ”Electronic Transport in Molecular Devices: the Role of Coherent and Incoherent Electron-Phonon Scattering,” Semicond. Sci. Technol, 19, 8357-8361 (2004). MA Reed, C. Zhou, C.J. Muller, T.P. Burgin, J.M. Tour, ”Conductance of a Molecular Junction,” Science, Vol. 278, 10 October 2000, pp. 252-254. MA Reed and J.M. Tour, ”Computing with Molecules,” Scientific American, June 20, 2000 J. Heath and M. Rather, ”Molecular Electronics,” Physics Today, May 2003. Y. Loo, J. W. P. Hsu, R. L. Willett, K W. Baldwin, K W. West and J. A Rogers. 2002. High-resolution transfer printing on GaAs surfaces using alkane dithiol monolayers. J. Vac. Sci. Technol. B, 20, 2853 O. S. Nakagawa, S. Ashok, C. W. Sheen, J. Martensson and D. L. Allara. 1991. GaAs interface with Octadecyl Thiol Self-Assembled Monolayer: Structural and Electronical Properties. Jag. 4. App. Phy. Vol. 30, pp. 3759-3762 A Koike and M. Yoneya. 1996. Molecular dynamics simulations of sliding friction of Langmuir-Blodgett onolayers. J. Chem. Phys. 105, 6060 226 57. 58. 59. 60. 61. 62. 63. 65. C. L. McGuiness, A Shaporenko, C. K Mars, S. Uppili, M. Zhamikov and D. L. Allara. 2006. Molecula Self-Assembly at Bare Semiconductor Surfaces: Prparation and Characterization of Highly Organized Octadecanethiolate Monolayers on GaAs(001). J. Am. Chem. Soc. 128, 5231 -5243 Z. Ge, D. G Cahill, and P. V. Braum. 2006. Thermal Conductance of Hydrophilic and Hydrophobic Interfaces. Phys. Rev. Lat 96, 186101 2. Wang, J. A Carter, A Lagutchev, Y. K Koh, N. H. Seong, D. G Cahill, and D. D. Dlott. 2007. Ultrafast Flash Thermal Conductance of Molecular Chains. Elm Vol. 317. no. 5839, pp. 787 -— 790 H. A Patel, S. Garde, and P. Keblinski. 2005. Thermal Resistance of Nanoscopic Liquid-Liquid Interfaces: Dependence on Chemistry and Molecular Architecture. Nam; Lottie. 5, 2225 Y. Yourdshahyan, H. K Zhang, and A M. Rappe. 2001. n-alkyl thiol head-group interactions with the Au(111) surface. Phys. Rev. B. 63, 081405 H. Gronbeck, A Curioni, and W. Andreoni. 2000. Thiols and Disulfides on the Au(111) Surface: The Headgroup—Gold Interaction. J. Am. Chem. Soc. 122, 3839 C. W. Sheen, J. Shi, J. Martensson, A N. Parikh, and D. L. Allara. 1992. A New Class of Organized Self-Assembled Monolayers: Alkane Thiols on GaAs(100). ,1. M. Chem. Soc. 114, 1514 W. Andersoni, A Curioni, H. Gronbeck 2000. Density Functional Theory Approach to Thiols and Disulfides on Gold: Au(111) Surface and Clusters. lntemational Journal of Quantum Chemi§t_ry. 80, 598 Scott Reese, Marye Anne Fox. 1998. Self-Assembled Monolayers on Gold of Thiols Incorporating Conjugated Terminal Groups. 4. Phfi. Chem. B, 102, 9820-9824 J. Hautman and M. L. Klein. 1989. Simulation of a monolayer of alkyl thiol chains. ,1. Cm. Phys. 91, 4994 227 67. 68. 69. 70. 71. 72. 73. 74. 75. 76. D. G Cahill, W. K Ford, K E. Goodson. 2003. Nanoscale thermal transport. 409ml oprplied Physics. 93(2), 793. S. Volz, J. B. Saulnier, M. Lallemand, B. Perrin, P. Depondt and M. Mareschal. 1996. Transient Fourier-law deviation by molecular dynamics in solid argon. Phys. Rev. B, 54,340 R. H. H. PM and Liam. 1994. Interplay of disorder and anharmonicity in heat conduction: Molecular dynamics study. Phys. Rev. _B_, 50, 15757 J. P. Crocombette, G Dumazer, and N. O. Hoang. 2007. Molecular dynamics modeling of the thermal conductivity of inadiated SiC as a function of cascade overlap. 4. App. Phys. 101, 023527 J. Li, L. Porter, 8. Yip. 1998. Atomistic modeling of finite-temperature properties of crystallinep—Src, ll. Thermal conductivity and effects of point defects. aMural of Nuclg Materials. 255, 139 A.J.H. McGaughey and J. Li, 2006, Mobcular Dynamics Prediction of the Thermal Resistance of Solid-Solid Interfaces in Superlattices, proceedings of IMECE, 13590. J. L. Barret and F. Chiaruttini, 2003, Kapitza Resistance at the Liquid-Solid Interface, mlecular Physics, 101, 1605-1610 L. Puech, G. Borrfait and B. Castaing, 1986, Mobility of the 3He Solid-Liquid Interface: Experiment and Theory, J. Low Temp. Phys, 62, 315-327 B. L. Huang, A J. H. McGaughey, and M. Kaviany. 2006. Thermal conductivity of metal-organic framework 5 (MOP-5): Part I. Molecular dynamics simulations. Int. 4. HQ and Mass Transfer. 50, pp. 405-411 A R. Leach. 1996. Molecular Modeling Principles and Applications. Addison Wesley Longman Ltd. MA I. H. Sung, D. E. Kim. 2005. Molecular dynamics simulation study of the nano—wear characteristics of alkanethiol self-assembled monolayers. 228 77. 78. 79. 80. 81. 82. 83. 85. Appl. Phys. A 81, 109 R. Mahafl‘y, R. Bhatia, and B. J. Garrison. 1997. Diffusion of a Butanethiolate Molecule on a Au(111) Surface. ,1. Phys. Chem. B. 101, 771 R. C. Lincoln, K M. Koliwad, and P. B. Ghate. 1967. Morse-Potential Evaluation of Second- and Third-Order Elastic Constants of Some Cubic Metals. Phys. Rev. 157. 463 M. Rieth: Neno-Engineering in Science and Technology (World Scientific, Singapore 2003) A K Rappe. C. J. Casewit, K S. Colwell, W. A Goddard, and W. M. Skiff. 1992. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. fig. 114, 10024 (unit conversion involved) L. Zhang, W. A Goodard, S. Jiang. 2002. Molecular simulation study of the c(4x2) supperlattice structure of alkanethiol self-assembled monolayers onAu(111). J. Chem. Phys. 117, 7342 V. B. Gohel, C. K Acharya and A R. Jani, 1984, On Phonon Dispersion in Noble Metals, J. Phys. F: Met. Phys, 15, 279-285 A Henry and G. Chen, 2008, High Thermal Conductivity of Single Polyethylene Chains Using Molecular Dynamics Simulations, Phys. Rev. Lett. 101, 235502 R. Y. Wang, R. A Segalman, A Majumdar. 2006. Room temperature thermal conductance of alkanedithiol self-assembled monolayers. AppL Phys. Lett. 89, 173113 Private discussion with Robert Wang from Mechanical Engineering of UC Berkeley. R. J. Stoner and H. J. Maris. 1993. Kapitza conductance and heat flow between solids at temperatures from 50 to 300 K Phys. Rev. B. 48, 16373 229 87. 88. 89. 91. 92. 93. 95. 97. 98. 99. R. M. Costescu, M. A Wall, md D. G Cahill. 2003. Thermal conductance of epitaxial interfaces. Phypig Revigg B. 87, 054302 A. N. Smith, J. L. Hostetler, and P. M. Norris. 2000. Thermal boundary resistance measurements using a transient thennoreflectance technique. Microscale Thermpmysical Epgineering. 4, 51 S. T. Huxtable, D. G Cahill, S. Shenogin, P. Keblinski, 2005, Relaxation of vibrational energy in fullerene suspensions, Chem. Phfi. Lett, 407, 129-134 P. Chantrenne and J. L. Barret; Journal of Heat Transfer. 2004, 126, pp. 577-585 H. J. Castejon; J. Phys. Chem. B. 2003, 107, pp. 826-828 C. Olischlger and J. C. Schon; Physical Review B. 1999, 59, pp. 4125-4133 E. Lussetti, T. Terao, and F. M. Plathe; J. Phys. Chem. B. 2007, 111, pp. 11516-11523 S. Mahajan, G Subbarayan, and B. G Sammakia; The Tenth Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronics Systems. 2006. pp. 1269—1274 F. Zhang, D. J. Isbister, aid D. J. Evans; lntemational Joumal of Themophysics. 2001, 22, pp. 135-147 F. M. Plathe; The Journal of Chemical Physics. 1997, 106, pp. 6082-6085 M. Zhang, E. Lussetti, L. E. S. de Souza, and F. M. Plathe; J. Phys. Chem. B. 2005, 109, pp. 15060-15067 B. C. Daly and H. J. Maris, 2002, Calculation of the thermal conductivity of superlattices by molecular dynamics simulations. Physig B, 316, 247-249 M. D. Muhlbaier. 2005. Self—Assembled Monolayers for Nanofabrication 230 100. 101. 102. 103. 104. 105. 106. 107. 108. 109. 110. and Nana-Scale Electronics. D. Segal, A Nitzan and P Hanggi. 2003. Thermal conductance through molecular wires. ,1. Chem. Phys. 119, 6840-6855 N. Shenogina, P. Keblinski, and S. Garde. 2008. Strong frequency dependence of dynamical coupling between protein and water. M Phys. 129, 155105 B. Hu, B. Li, and H. Zhao. 1998. Heat conduction in one-dimensional chains. Phys. Rey. E. 57, 2992-2995 C. W Sheen, J. K Shi, J. Martensson, 1992, A New Class of Organized Self-Assembled Monolayers: Alkane Thiols on GaAs (100), J. Am. Chem. Soc, 114, 1514-1515 A. Abdelghani, C. Jacquin, 2000, Structural characterization of GaAs/thiollelectrolyte interface, Materials Letters 46, 320-326 0. Voznyy and J. J. Dubowski, 2006, ”Structure, Bonding Nature, and Binding Energy of Alkanethiolate on As-Rich GaAs (001) Surface: A Density Functional Theory Study’, The Journal of Physical Chemistry B, 110, 23619-23622 S. Baroni, A Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G Ballabio, S. ScandoIo, G Chiarotti, P. Focher, A Pasquarello, K Laasonen, A Trave, R. Car, N. Marzari, A Kokalj, ht_tp:llwww.Mcf.orgl D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B. 41, 7892 (1990). P. Hohenberg and W. Kohn. Inhomogeneous Electron Gas. Phys. Rev., 136, 3864, 1964. J. P. Perdew, K Burke, and M. Emzerhof, Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865 (1996) J. P. Perdew. Electronic Structure of Solids ’91. Akademie Verlag, Berlin, 1991 231 111. 112. 113. 114. 115. 116. 117. 118. 119. 120. 121. 122. Edward F. c. Byrd and Betsy M. Rice, J. Phys. Chem. C 2007, 111. 2787-2798 G Zollo, J Tarus, R M Nieminen, J. Phys: Condens. Matter 16(2004) 3923-3932 H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations. Phys. Rev. B 13, 5188 (1976). West R C 1987-1988 CRC Handbook of Chemistry and Physics 68'" edn (Boce Raton, FL: CRC press) P. K Lam, M. L. Cohen, and G. Martinez, ‘Analytic relation between bulk moduli and lattice constants,’ Physical Review B, 1987, 35, 9190-9194 V. S. Vavilov, ”Handbook on the physical properties of 68, Si, GaAs and InP by A Dargys and J. Kundrotas”, PHYS-USP, 1996, 39 (7), 757-757 K Adlkofer and M. Tanaka. Stable Surface Coating of Gallium Arsenide with Octadecylthiol Monolayers. Langmuir 2001, 17, 4267 Y. Jun, K Y. Zhu, and J. W. P. Hsu. Formation of Alkanethiol and Alkanedithiol Monolayers on GaAs(001). Langmuir 2006, 22, 3627 F. P. Cometto, P. Paredes—Olivera, V. A Macagno, and E. M. Patrito, Density Functional Theory Study of the Adsorption of Alkanethiols on Cu(111), Ag(111), and Au(111) in the Low and High Coverage Regimes. J. Phys. Chem. B 2005, 109, 21737-21748 S. Franzen. Density functional calculation of a potential energy surface for alkane thiols on Au(1 1 1) as function of alkane chain length. Volume 381, 2003, 315-321 Y. Yourdshahyan and A M. Rappe. Structure and energetics of alkanethiol adsorption on theAu(111) surface. J. Chem. Phys. 117, 825 (2002) A Ferral, P. Paredes-Olivera, V. A Macagno, E. M. Patrito. Chenisorption and physisorption of alkanethiols on Cu(1 1 1). A quantum mechanical investigation. Surface Science, Volume 525, Issues 1-3, 10 February 2003, Pages 85-99 232 1 23. 124. 125. 126. 127. 128. 1 29. 130. 131. 132. 1 33. J. Gotrschaldt and 3. Hammer. A density functional theory study of the adsorption of sulfur, mercapto, and methylthiolate on Au(111). J. Chem. Phys. 116, 784 (2002) Henrik Gronbeck, Alessandro Curioni, and Wanda Andreoni. Thiols and Disulfides on the Au(111) Surface: The Headgroup-Gold Interaction. J. Am. Chem. Soc, 2000, 122 (16), pp 3839—3842 Harrell Sellers, Abraham Ulman, Yltzhak Shnidman, James E. Eilers. Structure and binding of alkanethiolates on gold and silver surfaces: implications for self-assembled monolayers. J. Am. Chem. Soc, 1993, 115 (21), pp 9389-9401 P. M. Morse, Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys. Rev. 1929, 34, 57-64. Altoe P, Stenta M, Bottoni A, Garavelli M. (2007). ”A tunable QMIMM approach to chemical reactivity, structure and physico-Chemical properties prediction”. Theor. Chem. Acc. 118: 219—240. Yang, R., and Chen, G., 2004, ”Thermal Conductivity Modeling of Periodic Two-Dimensional Nanocomposites,” Phys. Rev. B, 69, pp. 195316. D. W Brenner, O. A Shenderova, J. A Harrison, S. J. Stuart, B. Ni and S. B. Sinnott. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J. Phys: Condens. Matter 14 (2002) 783-802 GULP- a computer program for the symmetry adapted simulation of solids, J.D. Gale, JCS Faraday Trans, 93, 629 (1997) F. H. Stillinger and T. A Weber, 1985, ‘Computer simulation of local order in condensed phases of silicon’, Phys. Rev. B, 31, 5262 M. lchimura. Stillinger-Weber potentials for Ill-V compound semiconductors and their application to the critical thickness calculation for lnAs/GaAs. Physica Status Solidi (a), 153 (1996) 431 -437 h J .ioffe.rulSV S micondl aAsl 233 134. 135. 136. 137. 138. 1 39. 140. 141. 142. 143. 144. 145. 146. J. Tersoff, 1989, ‘Modeling solid-state Chemistry: lnteratomic potentials for multicomponent systems’, Physical Review B, 39, 5566-5568 J. Tersoff, 1988, ‘New empirical approach for the structure and energy of covalent systems’, Physical Review B, 37, 6991 -7000 M. Sayed, J. H. Jefferson, A B. Walker, and A G Cullis, 1995, ‘Molecular dynamics simulations of implantation damage and recovery in semiconductors’, Nucl. lnsb'. Meth. Phys. Res. B, 102, 218-222. Y. Touloukian, Thermophysical Properties of Matter, Vol. 3, Plenum, New York, 1 (1980) Komameni S 1992 . Nanocomposites. J. Mater. Chem. 21219—30 Njuguna B and Pielichowski K 2003 Adv. Eng. Mater. 5 769-78 Goettler L A, Lee K Y and Thakkar H 2007 Polym. Rev. 47 291-317 Chung D D L 2004 J. Mater. Sci. 39 2645-61 Carbon Nanotubes: Synthesis, Structure, Properties, and Applications; Dresselhaus, M. 3., Dresselhaus, G., Avouris, P., Eds; Springer. New York, 2000. Stankovich, S. et al. Stable aqueous dispersions of graphitic nanoplatelets via the reduction of exfoliated graphite oxide in the presence of poly(sodium 4-styrenesulfonate). J. Mater. Chem. 16, 155—158 (2006). Sasha Stankovich, Dmitriy A Dikin, Geofl‘rey H. B. Dommett, Kevin M. Kohlhaas, Eric J. Zimney, Eric A. Stach, Richard D. Piner, SonBinh T. Nguyen & Rodney S. Ruoff. Graphene-based composite materials. Nature 442, 282-286 Novoselov, K S. Et al. Two-dimensional gas of massless Dirac fermions in graphene. Nature 438, 197—200 (2005). Novoselov, K S. Et al. Electric field effect in atomically thin carbon films. Science 306, 666—669 (2004). 234 147. 148. 149. 150. 151. 152. 1 53. 154. 1 55. 156. 157. Zhang, Y., Tan, Y.-W., Stonner, H. L. & Kim, P. Experimental observation of the quantum Hall effect and Beny's phase in graphene. Nature 438, 201—204 (2005). Zhang, Y., Small, J. P., Amori, M. E. S. & Kim, P. Electric field modulation of galvanomagnetic properties of mesoscopic graphite. Phys. Rev. Lett. 94, 176803 (2005). Zhang, Y., Small, J. P., Amori, M. E. S. 8 Kim, P. Electric field modulation of galvanomagnetic properties of mesoscopic graphite. Phys. Rev. Lett. 94, 176803 (2005). Mansfield K F and Theodorou D N 1991 Macromolecules 24 4295-309 Daoulas K C, Harmandaris V A and Mavrantzas V G . 2005. Detailed Atomistic Simulation of a Polymer Melt/Solid Interface: Structure, Density, and Conformation of 8 Thin Film of Polyethylene Melt Adsorbed on Graphite. Macromolecules 38 5780-95 Harmandaris V A, Daoulas K C and Mavrantzas V G 2005. Molecular Dynamics Simulation of a Polymer Melt/Solid Interface: Local Dynamics and Chain Mobility in a Thin Film of Polyethylene Melt Adsorbed on Graphite. Macromolecules 38 5796—809 Goki Eda and Manish Chhowalla. Graphene-based Composite Thin Films for Electronics. Nano Lett, 2009, 9 (2), pp 814—818 US Patent 4704231 — Low-density graphite-polymer electrical conductors KURONEN A., TARUS J., NORDLUND K Defect creation by low-energy ion bombardment on GaAs (001) and Ge (001) surfaces Nuclear instruments 8 methods in physics research. Section B, Beam interactions with materials and atoms. 1999, vol. 153, pp. 209-212 L. A Girifalco, M. Hodak, and R S. Lee. 2000. Carbon nanotubes, buckyballs, ropes and a universal graphite potential. Physical Review B, 62,13104-13110 N. N. Greenwood and A Eamshaw, Chemistry of the elements 235 158. 1 59. 160. 161. 162. 163. 164. (Pergamon, New York, 1984) K Sun, M. A Stroscio, M. Dutta. 2009. Graphite C-axis thermal conductivity. Superlattices and Microstructures, 45, 60-64 R. J. Bell and P. Dean. 1970. Atomic vibrations in vitreous silica. Discuss. Faraday Soc. 50, 55—61 S. Himran, A Suwono. 1994. Characterization of Alkanes and Paraffin Waxes for Application as Phase Change Energy Storage Medium. Energy Sources, 16, 117-128 hpp:llwww.mstdiluvian.orfl~masonlmaterialslthennal conductiylty' .html P. Dashora, N. S. Saxena, M. P. Saksena, K Bela, K Sachdev, P. R. Pradhan and G. D. Ladiwala. 1992. A Theoretical Study of the Temperature Dependence of the Thermal Conductivity of Polymers. Physica Scripta, 45, 399-401 J. Xiang. Private discussion. C. L. Tim, 1993, Microscale Thermal Phenomena in Contemporary Technology, Robert Henry Thurston Lecture, ASME Vifinter Annual Meeting. 236