

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


FIRSTPRINCIPLES STUDY OF ELECTRONIC PROPERTIES OF ONE DIMENSIONAL NANOSTRUCTURES By JUNWEN LI Bachelor of Science in Applied Physics Qingdao University Qingdao, Shandong, China 1999 Master of Science in Physics Institute of Physics, Chinese Academy of Sciences Beijing, China 2002 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial ful llment of the requirements for the Degree of DOCTOR OF PHILOSOPHY May, 2010 FIRSTPRINCIPLES STUDY OF ELECTRONIC PROPERTIES OF ONE DIMENSIONAL NANOSTRUCTURES Dissertation Approved: Dr. John W. Mintmire Dissertation Advisor Dr. Timothy M. Wilson Dr. Yin Guo Dr. Nicholas F. Materer Dr. A. Gordon Emslie Dean of the Graduate College ii ACKNOWLEDGMENTS I am pleased to acknowledge my appreciation for those who have helped this dissertation take shape. First and foremost, my sincere thanks go to my advisor, Dr. John Mintmire, for spending enormous time and e orts in guiding my research. Every week I obtain valuable feedback from our discussion at individual and group meetings. I will always remember the scenes when he derived formulas and elucidated my puzzles on paper or the board. This dissertation originates from these gradual accumulations. He provides strong support for my research projects with his insight and understanding of the research eld and grand knowledge. His rigorous scholarship is what I myself hope to achieve in my future research career. Dr. Mintmire also sets up an example for me to follow as an individual. I would like to thank my committee members: Dr. Timothy Wilson, Dr. Yin Guo, and Dr. Nicholas Materer for their time and patience in serving my committee and reviewing my research proposal and dissertation. I appreciate their valuable comments and advice on my work. I would like to thank Dr. Jacques Perk, Dr. Xincheng Xie, and Dr. Robert Hauen stein for their help and encouragement. I would also like to thank Dr. Carter White and Dr. Daniel Gunlycke at Naval Research Laboratory and Dr. Vincent Munier at Oak Ridge National Laboratory for our collaboration. My thanks must go to all the former and current members in the research group and fellow graduate students for their friendship and indispensable support, but I have to single out Dr. Shelly Elizondo, Dr. Thushari Jayasekera, Vasantha Jogireddy, iii and Pillalamarri Pavan for their help, friendship, and making a pleasant work envi ronment. Finally, I would like to thank my parents, my grandfather, my elder brother and his family for always being there for me, even though whole continents and oceans stand between us. Thanks to my wife, Jing Zhu, for her love, patience, and encouragement. She has always been an invaluable source of strength and support for me during the years we have been together. iv TABLE OF CONTENTS Chapter Page 1 INTRODUCTION 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . 4 2 THEORETICAL BACKGROUND AND APPROACHES 5 2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 HartreeFock Approximation . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 DFT Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.2 ThomasFermi Model . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.3 HohenbergKohn Theorems . . . . . . . . . . . . . . . . . . . 15 2.3.4 KohnSham Formulation . . . . . . . . . . . . . . . . . . . . . 17 2.4 Exchange Correlation Approximations . . . . . . . . . . . . . . . . . 21 2.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.2 Local Density Approximation . . . . . . . . . . . . . . . . . . 21 2.4.3 PBE Generalized Gradient Approximation . . . . . . . . . . . 23 2.5 Gaussian Basis Sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.6 Helical Band Structure Methods . . . . . . . . . . . . . . . . . . . . . 28 2.7 Landauer Approach for Quantum Conductance Calculation . . . . . . 29 3 GRAPHENE NANORIBBONS 34 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 v 3.2 Planar Zigzag Graphene Nanoribbons . . . . . . . . . . . . . . . . . . 36 3.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.2 Energetic Results . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.3 Nonspinpolarized LDF Band Structures . . . . . . . . . . . . 40 3.2.4 Ferromagnetic LDF Band Structures . . . . . . . . . . . . . . 45 3.2.5 Antiferromagnetic LDF Band Structures . . . . . . . . . . . . 49 3.2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3 Planar Armchair Graphene Nanoribbons . . . . . . . . . . . . . . . . 55 3.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.3.2 Model Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.3.3 Local Density Functional Results . . . . . . . . . . . . . . . . 58 3.3.4 Consideration of Longer Range Interactions . . . . . . . . . . 61 3.3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4 Twist E ect in Electronic Properties of Armchair Graphene Nanoribbons 66 3.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.4.2 Computational Approach . . . . . . . . . . . . . . . . . . . . . 68 3.4.3 Band Structures and Neargap Wavefunctions . . . . . . . . . 70 3.4.4 Structure Parameters . . . . . . . . . . . . . . . . . . . . . . . 73 3.4.5 Tightbinding Model . . . . . . . . . . . . . . . . . . . . . . . 75 3.4.6 Band Gap Change Under Twist . . . . . . . . . . . . . . . . . 77 3.5 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 81 4 SURFACE PASSIVATION EFFECTS IN SILICON NANOWIRES 82 4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Computational Approaches . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3 Structure Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.4 Mulliken Population Analysis . . . . . . . . . . . . . . . . . . . . . . 91 4.5 Band Gap Change . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 vi 4.6 Electronic Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.7 Neargap States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.8 E ective Mass of Charge Carriers . . . . . . . . . . . . . . . . . . . . 107 4.9 Surface Defects E ect on Transport . . . . . . . . . . . . . . . . . . . 109 4.9.1 Computational Approach . . . . . . . . . . . . . . . . . . . . . 109 4.9.2 Transport Conductance . . . . . . . . . . . . . . . . . . . . . . 111 4.10 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 113 5 CONCLUSIONS 114 vii LIST OF TABLES Table Page 3.1 Total energies divided by 2N and energy di erences between di erent spin con gurations of zigzag GNRs with various ribbon widths . . . . 39 3.2 Energy gaps of zigzag GNRs with various ribbon widths . . . . . . . 53 3.3 Band gaps Egap of armchair GNRs with ribbon width N . . . . . . . 60 3.4 Binding states of edge C pz orbitals HOMO and LUMO. " (#) indicates the energy point shifts upward (downward) with increasing twist angle. + () indicates the bonding (antibonding) states. . . . . . . . . . . 76 3.5 Band gaps (eV) as a function of twist angle (Deg) in armchair graphene nanoribbons with various widths N . . . . . . . . . . . . . . . . . . . 80 4.1 Band gaps (eV) as a function of diameter for h100i silicon nanowires with various passivations. X represents H, OH, or CH3. . . . . . . . . 95 4.2 Band gaps (eV) as a function of diameter for h110i silicon nanowires with various passivations. X represents H, OH, or CH3. . . . . . . . . 95 4.3 E ective mass of silicon nanowires along h100i and h110i directions with various passivations . . . . . . . . . . . . . . . . . . . . . . . . . 108 viii LIST OF FIGURES Figure Page 3.1 Sample zigzag GNRs with ribbon width N = 6. The numbering of the unit cells is shown at the left side of the ribbon with only the unit cell labeled 0 shown in their entirety. Each unit cell is composed of 2N C atoms and 2 H atoms. The 2N C atoms in unit cell 0 are numbered as shown. The index of zigzag chains is indicated at bottom from left to right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 (a) Total energy di erences divided by 2N between nonspinpolarized and ferromagnetic states in zigzag GNRs as a function of ribbon width N. (b) Total energy di erences divided by 2N between ferromagnetic and antiferromagnetic states in zigzag GNRs as a function of ribbon width N. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4 Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 42 3.5 Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 43 ix 3.6 Calculated HOMO band orbital densities of zigzag GNRs of width N = 6, calculated at wavevectors (a) , (b) 0.875 , (c) 0.75 , and (d) 0.625 , respectively. Calculated orbital densities for the LUMO band are given for wavevectors (e) , (f) 0.875 , (g) 0.75 , and (h) 0.625 , respectively. The isovalue is taken to be 0.008 for all plots. . . . . . . 44 3.7 Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.8 Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.9 Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.10 Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 51 3.12 Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 52 3.13 Magnitude of the HOMOLUMO antiferromagnetic gap of the zigzag GNRs as a function of ribbon width N. . . . . . . . . . . . . . . . . . 53 x 3.14 (a) Symmetric armchair GNRs with ribbon width N = 7. Each unit cell is composed of 2N C atoms and 4 H atoms. (b) Staggered armchair GNRs with ribbon width N = 8. Each unit cell is composed of N C atoms and 2 H atoms. For (a) and (b) the numbering of the unit cells is shown at the left edge of each gure with only the unit cells labeled 0 shown in their entirety. The C atoms in unit cell 0 are numbered as shown. The index of dimer lines is indicated at bottom from left to right. 57 3.15 Calculated band structures of armchair GNRs of various widths. (a), (b), and (c) correspond to ribbon with N = 7, 8, and 9, respectively. (X) indicates the center (edge) of the rst Brillouin zone. . . . . . . 59 3.16 LDF energy gaps of the armchair GNRs as a function of width N. . . 60 3.17 Magnitude of the HOMOLUMO gap of the N armchair GNRs from the TBM calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.18 LDF results for Egap compared to the corresponding tightbinding model (TBM) results obtained from Eq. (3.3) with V1 = 3:2 eV and V3 = 0:3 eV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.19 (a) Armchair GNR with width N = 8. Each helical unit cell is com posed of one zigzag chain terminated by H atoms. The N C atoms in a helical cell are numbered as shown. The bonds are labeled as a1a11. (b) Axial view of 10 twisted armchair GNR with width N = 8. . . . 69 3.20 Electronic band structure of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . . . . . . . . . . . . . . . 71 3.21 The LUMO and HOMO wave function of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . The ribbons are vertically oriented. . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.22 CC bond lengths of the twisted armchair GNRs with width N = 8 as a function of twist angle . . . . . . . . . . . . . . . . . . . . . . . . 74 xi 3.23 Variation of energy gaps of armchair GNRs as a function of twist angle . (a), (b), and (c) represent the armchair GNRs of width N = 3m, 3m + 1, and 3m + 2, respectively. . . . . . . . . . . . . . . . . . . . . 79 4.1 Crosssection views of the studied HSiNWs. h100i SiNWs of diameter (a) 0.43 nm, (b) 0.87 nm, (c) 1.3 nm, and (d) 1.73 nm, respectively. h110i SiNWs of diameter (e) 0.73 nm, (f) 1.09 nm, and (g) 1.46 nm, respectively. The goldenyellow and silver balls represent Si and H atoms, respectively. In OHSiNWs and CH3SiNWs, all H atoms are replaced with OH and CH3 groups, respectively. . . . . . . . . . . . . 87 4.2 Band gaps as a function of diameter with hydrogen passivation for silicon nanowires along h100i. The square and blue line indicate the LDF results and the uptriangle and red line indicate the PBE results. 88 4.3 SiSi bond length as a function of the distance from bond center to wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm passivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4 Net charge as a function of distance from the atom to the wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm pas sivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. Si, H, O, and C are indicated by +, , , and , respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.5 Band gaps as a function of diameter with various surface passivations for silicon nanowires along (a) h100i, (b) h110i. . . . . . . . . . . . . 94 xii 4.6 Electronic band structures of h100i SiNWs (diameter d = 1.73 nm) passivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. . . . . . . . . . . . . 97 4.7 Electronic band structures of h110i SiNWs (diameter d = 1.46 nm) passivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. . . . . . . . . . . . . 98 4.8 The HOMO and LUMO orbital density of the h100i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respec tively. Here red and blue represent HOMO and LUMO density, respec tively. The contour is at 10% of the maximum value. . . . . . . . . . 100 4.9 The HOMO and LUMO orbital density of the h110i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respec tively. Here red and blue represent HOMO and LUMO density, respec tively. The contour is at 10% of the maximum value. . . . . . . . . . 101 4.10 The HOMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 103 4.11 The LUMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 104 4.12 The HOMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 105 4.13 The LUMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 106 xiii 4.14 Model system for the conductance calculation. The system is composed of three regions: lefthand lead region, central conductor region, and righthand lead region. . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.15 The transport conductance through the OHpassivated region: The solid line is the conductance of the pure HSiNWs. The dashed line shows the conductance of the HSiNWs with OH defects on one helical cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 xiv CHAPTER 1 INTRODUCTION 1.1 Overview Nanotechnology deals with structures between 1 nanometer and 100 nanometer in size. The materials at this scale can display properties very di erent from their bulk counterparts. These unique nanoscale properties originate from size e ects such as the enhanced surfacetovolume ratio that causes surface atoms to experience po tentials di erent from those in the bulk. In recent years the rapid development of fabricating and manipulating materials at nanometer scale has opened up a great deal of opportunities in a variety of scienti c research e orts and practical applications. We can consider di erent numbers of spatial dimensions at nanometer scale. De pending on how many dimensions are con ned, the nanomaterials are termed as quantum well, quantum wire, and quantum dot for constraints on one, two, and three dimensions, respectively. Of these di erent systems, the quantum wire of one dimensional extension represents the simplest transport path and is of especially great signi cance in electronic applications. Two types of onedimensional nanostructures were investigated in this dissertation: graphene nanoribbons and silicon nanowires. Our study on graphene nanoribbons was motivated by the discovery of new exper imental methods for graphene, a single sheet of graphite. Novoselov, et al.,[1, 2] using a micromechanical cleavage method to extract a single sheet from threedimensional graphite, obtained the twodimensional form of carbon, i.e., graphene, which is stable under ambient conditions. The quality of the samples prepared were good enough so that the ballistic transport[1] and quantum Hall e ect[3] were observed easily. The 1 ballistic transport makes this new material a promising candidate for future electronic applications such as ballistic elde ect transistors (FETs). Although experimental work on graphene systems is recent, the theoretical study of graphene started long ago, with Wallace's tightbinding study in 1947 as an early work. [4] After that many properties of graphite and nanotubes were derived using graphene as starting point. The experimental achievement of graphene in 2004 attracted great attention immedi ately. Twodimensional graphene has a zero band gap with linear energy dispersion near the Fermi level and crossing at the Dirac point between valence and conduc tion bands. For practical electronic application in semiconductor industry a band gap could be induced by fabricating graphene nanoribbons. The ribbon structures have already been successfully obtained with various techniques such as lithographic patterning,[5] chemical vapor deposition,[6] Joule heating,[7] and unzipping the car bon nanotubes (CNTs).[8, 9] The other onedimensional structures investigated, silicon nanowires, have a wide variety of applications from use in batteries to use in energy conversion. Because of the geometry of silicon nanowire arrays, they are capable of in ating 4 times their normal size when absorbing lithium ions, which enables the new kind of battery to hold 10 times the charge of existing lithiumion batteries.[10] The onedimensional structure of silicon nanowires can allow pin coaxial solar cells consisting of a ptype silicon nanowire core surrounded by i and ntype silicon shells.[11] An advantage of this core/shell architecture is that carrier separation takes place in the radial rather than the axial direction, with a carrier collection distance smaller or comparable to the mi nority carrier di usion length. Si nanowire arrays show promise as highperformance, scalable thermoelectric materials because of the high electron conductivity and low phonon conductivity in silicon nanowires with rough surface.[12, 13] The promising applications of these nanostructures require a better theoretical understanding of their electronic properties. Our computational simulations of elec 2 tronic properties of graphene nanoribbons and silicon nanowires are of great signi  cance from both fundamental and practical points of view. 3 1.2 Organization of the Dissertation Chapter 2 presents a brief introduction of the theoretical background and com putational approaches utilized throughout this work. We brie y review the theo retical formalism of density functional theory, helical band structure methods, and the transport conductance approach with Landauer formula implemented by using nonequilibrium Green's function method. Chapter 3 is dedicated to the simulation results from a study carried out on graphene nanoribbons. This theoretical study was primarily motivated by the exper imental implementation of isolated graphene sheets which are a single layers of carbon atoms connected with hexagonal patterns. In this work, we investigated two kinds of graphene nanoribbons, namely zigzag and armchair graphene nanoribbons according to the edge shape. For zigzag nanoribbons, the calculations on energetic analysis and corresponding band structures were carried out. For armchair nanoribbons, we studied the electronic structure dependence on the ribbon width and on the torsional deformation. We have four papers related to the ribbon studies. [14{17] Chapter 4 is devoted to the investigation carried out involving surface passivation e ect in silicon nanowires. Because of the enhanced surfacetovolume ratio in silicon nanowires compared to bulk silicon, the electronic properties of the nanowires are strongly dependent on the surface substituents. We studied three kinds of passivating functional groups such as hydrogen, hydroxyl, and methyl. We also investigated the hydroxyl surface defects e ect on transport property of pure hydrogenpassivated silicon nanowires. We have two papers based on the simulations presented in this chapter.[18, 19] Chapter 5 summarizes this dissertation, focusing on the simulation ndings and the possible further investigation based on the present calculations. 4 CHAPTER 2 THEORETICAL BACKGROUND AND APPROACHES 2.1 Overview All band structure calculations reported in this work were carried out by us ing the Helical Nanostructures (HENS) parallax code developed by Dr. Mintmire at Oklahoma State University. Di erent from regular band structure codes that use translational symmetry, this approach takes advantage of the helical symme try in onedimensional nanostructures such as polymers, nanoribbons, nanotubes, and nanowires. Taking account of the helical symmetry can greatly reduce the unit cell size, and therefore can achieve major computational savings. Fujita and Imamura[20, 21] rst proposed this idea of utilizing helical symmetry when studying the electronic structures of polymers with tight binding model. And this approach was developed for rstprinciples methods by Dr. Mintmire while at the Naval Research Laboratory. Mintmire, et al.,[22] using this helical band structure code, successfully predicted that the armchair singlewalled carbon nanotubes would be metallic four years prior to the experimental veri cation. Now we apply this code in the study of graphene nanoribbons and silicon nanowires. The zerobias conductance was calculated within a Landauer approach.[23, 24] The transport code was developed by using the nonequilibrium Green's function (NEGF) formalism within a linear response regime. The approach is capable of taking any general Hamiltonian that can be described with a localized orbital basis, such as orthogonal, nonorthogonal tightbinding models, and rstprinciples results by making use of Gaussian orbitals as in the study of surface defects e ect in silicon 5 nanowires. The e ective Hamiltonian matrix elements and orbital overlap matrix elements used in the investigation of silicon nanowires were generated by the HENS parallax selfconsistent eld (SCF) calculations. In this chapter, we will brie y describe the background on which the codes are based, including the density functional theory (DFT), the helical band structure method and the NEGF approach. Before that the HartreeFock method is reviewed in order to gain some idea on the exchange interaction and the advantage of DFT when dealing with manyelectron system. The discussion in this chapter is mainly adapted from Refs. [25{27]. 6 2.2 HartreeFock Approximation For an isolated Nelectron atomic or molecular system, in the BornOppenheimer approximation the timeindependent nonrelativistic Schr odinger equation is given by ^H = E ; (2.1) where E is the electronic energy, = (x1; x2; :::; xN) is the Nelectron wavefunction, and ^H is the Hamiltonian operator composed of three terms, ^H = ^ T + ^ Vext + ^ Vee; (2.2) where ^ T = XN i=1 1 2 r2i (2.3) is the kinetic energy operator, ^ Vext = XN i=1 XM =1 Zi jr R j (2.4) is the electronnucleus attraction energy operator, and ^ Vee = XN i6=j 1 jri rj j (2.5) is the electronelectron repulsion energy operator. When additional elds are present, ^ T and ^ Vee remain unchanged and extra terms appear in ^ Vext. The total energy ETOT of the system is composed of the electronic energy E and the nucleusnucleus repulsion energy Vnn = XM < Z Z R : (2.6) So we may have ETOT = E + Vnn: (2.7) 7 The coordinate xi of electron i comprises spatial coordinates ri and spin freedom si. Atomic units are employed here and throughout this dissertation (unless otherwise speci ed): the length unit is the Bohr radius, a0 = 0:5292 A, the charge unit is the charge of the electron, e, and the mass unit is the mass of the electron, me. The ^ T and ^ Vext in Eq. (2.2) are singleparticle operators and the ^ Vee is a two particle operator from which the complications of manybody problem arises. We need to make approximation on the manybody wavefunction (x1; x2; :::; xN). The simplest way to approximate the wavefunction is through the Hartree approximation, where the manybody wavefunction is replaced by a product of singleparticle orbitals, i(xi), (x1; x2; :::; xN) = 1 p N 1(x1) 2(x2) N(xN); (2.8) where i(xi) is a combination of spatial function i(ri), and spin function (si), such that i(xi) = i(ri) (si); (2.9) where = ; indicate spinup and spindown, respectively. However, the electrons are Fermions and should obey the Pauli exclusion principle which states that one quantum state only can be occupied by at most one electron. The wavefunction should be antisymmetric with respect to the interchange of any two electron coordi nates xi and xj . Eq. (2.8), however, does not satisfy (x1; :::; xi; :::; xj ; :::; xN) = (x1; :::; xj ; :::; xi; :::; xN): (2.10) To account for Pauli principle the HartreeFock approximation was proposed by writing the wavefunction as an antisymmetrized product of orbitals explicitly. The HartreeFock wavefunction HF amounts to a linear combination of the terms in Eq. (2.8) including all permutations of the electron coordinates with the corresponding 8 weights 1, which can be expressed as a Slater determinant, HF = 1 p N! 1(x1) 1(x2) 1(xN) 2(x1) 2(x2) 2(xN) ... ... ... ... N(x1) N(x2) N(xN) : (2.11) The HartreeFock energy can be evaluated by taking the expectation value of the Hamiltonian with respect to the Slater determinant Eq. (2.11). This yields EHF = h HF j ^H j HF i = XN i=1 Z i (r) 1 2 r2 + vext(r) i(r) dr + 1 2 XN i=1 XN j=1 Z Z j i(r1)j2j j(r2)j2 jr1 r2j dr1dr2 1 2 XN i=1 XN j=1 Z Z i (r1) i(r2) j (r2) j(r1) jr1 r2j si;sj dr1dr2: (2.12) The last term is of signi cant interest since it arises from the antisymmetric nature of the HartreeFock wavefunction. It vanishes when si 6= sj , which is a result of the Pauli exclusion principle. Therefore it is termed as exchange energy. In order to obtain the total energy of the system an extra term coming from repulsion energy between the nucleus must be added to Eq. (2.12). In most cases, we are more concerned with the ground state of a system. Even from the ground state property of system, we could get a great deal of information on the property of excited states in room temperature since the system deviates not too much from the ground states. To obtain the HartreeFock ground state energy EHF 0 we could minimize Eq. (2.12) with respect to the orbitals, subject to the constraint that the orbitals remain orthonormal (h ij ji = ij). The minimization procedure carried out with the Euler 9 Lagrange method yields the corresponding stationary condition given by EHF 0 XN i=1 XN j=1 "i (h ij ji 1) ! = 0; (2.13) where "i is the Lagrange multiplier. The corresponding Euler equations are written as 1 2 r2 + vext(r) + XN j=1 Z j j(r0)j2 jr r0j dr0 ! i(r) XN j=1 Z i(r0) j (r0) j(r) jr r0j si;sj dr0 = "i i(r); (2.14) which is called HartreeFock equations. In general the HartreeFock equations can not be solved analytically. One excep tion is for the homogeneous electron gas, where the solutions are plane wave functions because of the constant external potential everywhere. Otherwise, the HartreeFock equations are solved using an iterative method. Because the desired orbitals are re quired to construct the oneelectron e ective potential, this process is known as the selfconsistent eld procedure. The selfconsistent procedure starts with an initial guess for the orbitals, and successive iterations are performed with new orbitals gen erating the new potentials until the convergence is achieved. The converged orbitals are the ground state orbitals for that system within the HartreeFock approximation. HartreeFock theory is not an exact theory because it only considers a single determinant for the electron wavefunction. The only case where a single determinant is the exact solution is for a noninteracting system of electrons. In real systems the motions of electrons are more correlated than the mean eld description provided by HartreeFock approximation. The interaction energy missed in HartreeFock approximation is termed as the correlation energy EC, EC = E0 EHF ; (2.15) 10 where E0 is the exact ground state energy. Since HartreeFock is a variational ap proach, EHF E0 always holds and the correlation energy takes nonpositive value. A natural way to include correlation e ects beyond the HartreeFock approxima tion is to represent the manyelectron wavefunction as a linear combination of Slater determinants corresponding to ground state and excited states. These post Hartree Fock methods, such as con guration interaction, coupledcluster and M llerPlesset theory have been extensively developed in quantum chemistry. Although the pre cision may be systematically improved by including more and more excited Slater determinants, the computational cost increases dramatically with the number of ex citation levels. Consequently, these post HartreeFock approaches are limited to small systems such as atoms and small molecules. 11 2.3 DFT Methods 2.3.1 Overview The main idea of density functional theory (DFT) is to describe a system composed of N interacting electrons via the electron density rather than the manyelectron wavefunction like the Slater determinant used in HartreeFock approximation. This means that the basic variable of the system depends only on 3 spatial coordinates (x, y, and z) rather than 3N degrees of freedom as in HartreeFock approximation. Additionally, the computational costs are relatively low compared to those methods based on the complicated manyelectron wavefunction such as HartreeFock theory and its descendants. In 1927 Thomas and Fermi rst proposed a model expressing the electronic energy as a functional of electron density.[28, 29] In the original idea they derived a di erential equation for the electron density without resorting to oneelectron orbitals. The DFT has its conceptual roots in the ThomasFermi model, but the rm theoretical foundation was set up by two HohenbergKohn theorems 37 years after the proposal of ThomasFermi method. 2.3.2 ThomasFermi Model Thomas and Fermi proposed that the electronic energy of a system can be ex pressed as the functional of electron density. In this idea the kinetic, exchange, and correlation contributions can be constructed from the model study on the homoge neous electron gas but should be dependent on the position. The electron density (r) is the central variable, given by (r) = N Z Z dr2 drN (r; r2; :::; rN) (r; r2; :::; rN): (2.16) 12 The total energy of a system is expressed as a functional ETF [ (r)], which is given by ETF [ (r)] = CF Z (r)5=3 dr + Z (r)vext(r) dr + 1 2 Z Z (r1) (r2) jr1 r2j dr1dr2: (2.17) The rst term in Eq. (2.17) is the electronic kinetic energy associated with non interacting homogeneous electron gas. This form could be obtained by integrating the kinetic energy density of a homogeneous electron gas t0[ (r)], TTF [ (r)] = Z t0[ (r)] dr; (2.18) where t0[ (r)] is obtained by summing all of the free electron energy states "k = k2 2 , up to the Fermi wave vector kF = [3 2 (r)]1=3, t0[ (r)] = 2 8 3 Z k2 2 nk dk = 1 2 Z kF 0 k4 dk (2.19) with nk is the density of states in reciprocal space. This leads to the form given in Eq. (2.17) with coe cient CF = 3 10(3 2)2=3. The second term is the classical electrostatic energy of attraction between the nucleus and the electrons, where vext(r) is the static Coulomb potential arising from the nucleus vext(r) = XM =1 Z jr R j : (2.20) Finally, the third term in Eq. (2.17) represents the electronelectron interactions of the system. In the ThomasFermi model this term contains only the classical Coulomb repulsion energy between electrons, known as the Hartree energy. In order to obtain the ground state density and energy of a system, we may min imize the ThomasFermi energy functional of Eq. (2.17) with constraint of conserved total number of electrons N. Applying the EulerLagrange method to Eq. (2.17) leads to the stationary condition ETF [ (r)] Z (r) dr N = 0; (2.21) which yields the socalled ThomasFermi equations 5 3 CF (r)2=3 + vext(r) + Z (r0) jr r0j dr0 = 0; (2.22) 13 which can be solved directly to obtain the ground state density. ThomasFermi theory su ers from many di culties. One of them is that it does not predict bonding between atoms, so molecules and solids can not exist in this theory. The main source of error is a crude approximation for the kinetic energy. The kinetic energy represents a substantial portion of the total energy of a system and so even small errors prove disastrous. Another shortcoming is the oversimpli ed description of the electronelectron interactions, which are treated classically and so do not take into account quantum mechanical e ects such as the exchange interaction. Shortly after the introduction of ThomasFermi theory, Dirac[30] developed an approximation for the exchange interaction based on the homogeneous electron gas. The resulting formula is simple, and is also a local functional of the density, EX[ (r)] = CX Z (r)4=3 dr (2.23) with CX = 3 4 3 1=3 . When exchange interaction is treated like Eq. (2.23), the theory is called ThomasFermiDirac Model. Correlation can also be easily included by using any local approximation derived from homogeneous electron gas. One commonly used was proposed by Wigner[31], EC[ (r)] = 0:056 Z (r)4=3 0:079 + (r)1=3 dr: (2.24) The ThomasFermi model was actually too crude, mainly because the approxima tion used for the kinetic energy of the electrons was unable to sustain bound states. However, it set up the basis for the later development of density functional theory. 14 2.3.3 HohenbergKohn Theorems The ThomasFermi approach was developed in the hopes that the energy can be expressed exclusively in terms of the electron density. This idea, however, was intuitive at the time. Until 1964, Hohenberg and Kohn[32] proved two theorems that put the ThomasFermi model on solid mathematical grounds. Theorem 1. The external potential vext(r) is determined, within a trivial additive constant, by the electron density (r). Let us suppose there are two di erent external potentials, vext;1(r) and vext;2(r) such that they correspond to same ground state electron density (r). Let 1 and E1 = h 1j ^H 1j 1i be the ground state wavefunction and ground state energy of the Hamiltonian ^H 1 = ^ T + ^ Vext;1 + ^ Vee. Let 2 and E2 = h 2j ^H 2j 2i be the ground state wave function and ground state energy of the Hamiltonian ^H 2 = ^ T + ^ Vext;2 + ^ Vee. We assumed that di erent Hamiltonians correspond to di erent ground state wave functions, i.e., 1 6= 2. Applying the variational principle we may have: E1 < h 2j ^H 1j 2i = h 2j ^H 2j 2i + h 2j ^H 1 ^H 2j 2i (2.25) = E2 + Z (r)[vext;1(r) vext;2(r)] dr (2.26) Similarly, we have E2 < h 1j ^H2j 1i = h 1j ^H 1j 1i h 1j ^H 1 ^H 2j 1i (2.27) = E1 Z (r)[vext;1(r) vext;2(r)] dr (2.28) Therefore adding the two inequalities leads to the result, E1 + E2 < E2 + E1; (2.29) which is a contradiction, and as a result there can not be two di erent external potentials that correspond to the same electron density for the ground state, unless they di er by a trivial additive constant. 15 Theorem 2. The ground state energy can be obtained variationally : the density that minimizes the total energy is the exact ground state density. As just shown, (r) determines vext(r), N and vext(r) determine ^H and therefore . This ultimately means is a functional of (r), and so the expectation value of any operator ^O is also a functional of (r), i.e., O[ (r)] = h [ (r)]j^O j [ (r)]i: (2.30) Assume we have 0(r), which determines its own v0 ext(r) and wavefunction 0. Then 0 could be taken as a trial function for the system of interest having external potential vext and corresponding (r). Thus, h 0j ^H j 0i = h 0j ^ Tj 0i + h 0j ^ Veej 0i + h 0j ^ Vextj 0i = E[ 0(r)] E[ (r)] (2.31) based on the variational argument. 16 2.3.4 KohnSham Formulation Although the HohenbergKohn theorems are extremely powerful, they do not o er a way to calculate the ground state density of a system in practice. About one year after the seminal DFT paper by Hohenberg and Kohn, Kohn and Sham (KS)[33] devised a simple method for carrying out DFT calculations, that retains the exact nature of DFT. Since the electron density is the central variable rather than the wavefunction, we could construct our electron density from orbitals which are obtained from easily solved system such as the noninteracting system. The KohnSham formulation centers on mapping the full interacting system with the real potential, onto a ctitious noninteracting system where the electrons are moving in an e ective singleparticle KS potential vKS(r). The KohnSham method is still exact since it yields the same ground state electron density as the real system, but greatly facilitates the calculation. First consider the variational problem presented in the second HohenbergKohn theorem. The ground state energy of a many electron system can be obtained by min imizing the energy functional, subject to the constraint that the number of electrons N is conserved, which leads to F[ (r)] + Z vext(r) (r) dr Z (r) dr N = 0; (2.32) where F[ (r)] = T[ (r)] + Vee[ (r)] is called universal functional since they do not depend on the external potentials. The Euler equation is given by = F[ (r)] (r) + vext(r); (2.33) where is the Lagrange multiplier associated with the constraint of conserved N. The idea of Kohn and Sham was to set up a system where the kinetic energy could be determined exactly, since this was a major problem in ThomasFermi theory. This 17 was achieved by resorting to a noninteracting system of electrons. The correspond ing ground state wavefunction KS for this type of system is given exactly by a determinant of singleparticle orbitals i(ri), KS = 1 p N! det[ 1(r1) 2(r2)::: N(rN)]: (2.34) The universal functional F[ (r)] was then partitioned into three terms, the rst two of which are known exactly and constitute the majority of the energy, the third being a small unknown quantity, F[ (r)] = TS[ (r)] + EH[ (r)] + EXC[ (r)]: (2.35) TS[ (r)] is the kinetic energy of a noninteracting electron gas of density (r), EH[ (r)] is the Hartree energy of the electrons EH[ (r)] = 1 2 Z Z (r1) (r2) jr1 r2j dr1dr2; (2.36) and EXC[ (r)] is the exchange correlation energy, which contains the di erence be tween the exact and noninteracting kinetic energies and also the nonclassical contri bution to the electronelectron interactions, of which the exchange energy is a part. In the KohnSham prescription the Euler equation now becomes = TS[ (r)] (r) + vKS(r) (2.37) where the e ective KS potential vKS(r) is given by vKS(r) = vext(r) + vH(r) + vXC(r) (2.38) with the Hartree potential vH(r), vH(r) = EH[ (r)] (r) = Z (r0) jr r0j dr0; (2.39) and the exchange correlation potential vXC(r), vXC(r) = EXC[ (r)] (r) : (2.40) 18 In KohnSham theory the several potential terms were just rearranged to make up vKS. So the density obtained when solving the ctitious noninteracting KohnSham system is the same as the exact ground state density. The ground state density is obtained in turn by solving the N oneelectron Schr odinger equations, 1 2 r2 + vKS(r) i(r) = "i i(r); (2.41) where "i are Lagrange multiplier corresponding to the orthonormality of the N single particle states i(r), and the density is constructed from (r) = XN i=1 j i(r)j2: (2.42) The noninteracting kinetic energy TS[ (r)] is therefore given by TS[ (r)] = 1 2 XN i=1 Z i (r)r2 i(r) dr: (2.43) Since vKS(r) depends on the density through the exchange correlation and Hartree potentials, the KohnSham equations must be solved selfconsistently as in the Hartree Fock method. In order to handle the kinetic energy in an exact manner, N equations have to be solved in KohnSham theory to obtain the set of Lagrange multiplier "i, as opposed to one equation that determines when solving for the density directly, as in the ThomasFermi approach. However an advantage of the KohnSham method is that as the complexity of a system increases, with N increasing, the problem becomes no more di cult, only the number of singleparticle equations to be solved increases. Although exact in principle, KohnSham theory is approximate in practice because of the unknown exchange correlation functional EXC[ (r)]. An implicit de nition of EXC[ (r)] can be given by, EXC[ (r)] = T[ (r)] TS[ (r)] + Eee[ (r)] EH[ (r)] (2.44) where T[ (r)] and Eee[ (r)] are the exact kinetic and electronelectron interaction energies respectively. The intention of Kohn and Sham was to make the unknown 19 contribution to the total energy of the noninteracting system as small as possible, and this is indeed the case with the exchange correlation energy, however it is still an important contribution since the binding energy of many systems is about the same size as EXC[ (r)], so an accurate description of exchange and correlation is crucial for the prediction of binding properties. 20 2.4 Exchange Correlation Approximations 2.4.1 Overview While DFT is an exact theory of ground state properties, practical applications of DFT must be based on approximations for the unknown exchange correlation po tential which describes the e ects of Pauli exclusion principle and Coulomb potential beyond the pure electrostatic interaction between electrons in an average sense. If the exact exchange correlation potential is obtained, the manybody problem can be solved exactly. All exchange correlation functionals can be written in a general form EXC[ (r)] = Z (r)"XC(r) dr; (2.45) where "XC(r) is the exchange correlation energy density. In DFT, the exchange correlation functionals can be approximated in di erent levels by the number and kind of their local ingredients.[34] The simplest one is the local density approximation in which only the local electron density is considered. In the high level approximation such as generalized gradient approximation (GGA), the gradient of electron density r (r) is also included. We will discuss local density approximation (LDA) and GGA especially PBE with more details because these two kinds of functionals are extensively used in our simulations. 2.4.2 Local Density Approximation In 1951, Slater[35] suggested a model in which the kinetic energy would be treated as in the HartreeFock model, but where the exchange term was replaced by a func tional of electron density vX[ (r)]. This work lead to what is known as the X SCF method. In the derivation of Slater's oneelectron exchange potential, we follow the 21 procedure of Kohn and Sham.[33] We may begin by writing the HartreeFock exchange operator in the form of an equivalent potential acting on the kth wave function, vxk(x) = XN k0=1 Z k(x) k0(x0) k0(x) k(x0) jr r0j dx0 k(x) k(x): (2.46) Next we make an approximation and simpli cation assuming the wave functions can be approximated by plane waves as in free electron gas. This leads to vxk(r) = kF (r) 1 + k2F (r) k2 2kkF (r) ln k + kF (r) k kF (r) (2.47) with kF (r) = (3 2 (r))1=3. Then, we average vxk over the occupied state k, which results in vX[ (r)] = 3 2 (3 2)1=3 (r)1=3; (2.48) which is the original Slater's exchange potential. Since electron density adjustments are mainly a ected by redistribution of the electrons near the Fermi level, it is reasonable to take k = kF (r) in Eq. (2.47), which is equivalent to taking the e ective exchange potential for a state at the top of the Fermi distributions. This leads to vX[ (r)] = 1 (3 2)1=3 (r)1=3; (2.49) which is the expression for exchange potential that Kohn and Sham derived in 1965. Eq. (2.49) is di erent from Slater's by a factor of 2 3 . The paper published by Kohn and Sham in 1965 is not the rst addressing this di erence of twothirds. Actually in 1954, G asp ar already obtained the same kind of dependence of the exchange energy on the electron density as that of Kohn and Sham. The model, in which vX[ (r)] is proportional to (r)1=3, became known as the X method, where the exchange term is written in the form, 3 2 (3 2)1=3 (r)1=3 (2.50) 22 with = 1 for Slater's model and = 2=3 for that of G asp ar, Kohn, and Sham. In Eq. (2.50), we only include the exchange e ect. Together with this exchange potential, a commonly used correlation formula is that of Perdew and Zunger[36] which makes use of accurate quantum Monte Carlo data for the homogeneous electron gas generated by Ceperley and Alder[37] to x the coe cients in the interpolation formula. Our local density functional (LDF) approach employs the G asp arKohnSham exchange potential and neglects the correlation e ect. This approach is also commonly referred to as the local density approximation (LDA). 2.4.3 PBE Generalized Gradient Approximation For exchange correlation energy the LDA makes use of the result obtained from homogeneous electron gas at each point irrespective of the nonhomogeneity of the real electron density and is the simplest approximation for exchange and correla tion functionals. For real system of nonhomogeneous electron density the exchange correlation energy can be signi cantly di erent from that of homogeneous electron gas. The gradient and higher spatial derivatives of the electron density are needed to account for this deviation. Perdew and colleges made great contributions on the development of GGA func tionals. [34, 36, 38{42] They introduced an analytic function known as the enhance ment factor, FXC[ (r);r (r)], that accounts for the nonhomogeneity of electron density and modi es the LDA energy density through EGGA XC [ (r)] = Z (r)"hom XC [ (r)]FXC[ (r);r (r)] dr: (2.51) Usually the GGA enhancement factor is written in terms of the Seitz radius, rs, which is related to the electron density (r) by rs = 3 4 (r) 1=3 ; (2.52) 23 and the dimensionless reduced density gradient s(r), s(r) = jr (r)j 2kF (r) (r) (2.53) with kF being the Fermi wave vector, kF (r) = 3 2 (r) 1=3 : (2.54) At present the most popular GGA functional in physics community is PBE.[38] The PBE functional EPBE XC is composed of two parts : PBE exchange EPBE X and PBE correlation EPBE C functionals. The exchange PBE functional is written in the following form, EPBE X [ (r)] = Z d3r (r)"hom X ( )FX(s) (2.55) with FX(s) = 1 + k k 1 + s2=k (2.56) where k = 0:804; = 0:21951. This exchange energy obeys the spinscaling relation ship,[43] EPBE X [ ; ] = 1 2 EPBE X [2 ] + 1 2 EPBE X [2 ]: (2.57) And the PBE correlation energy functional is given by EPBE C [ ; ] = Z d3r (r)["hom C (rs; ) + H(rs; ; t)]; (2.58) where the nonlocal part H(rs; ; t) depends on the parameter t including the density gradient via H = 3 ln 1 + t2 1 + At2 1 + At2 + A2t4 (2.59) where = 0:066725; = 0:03191, and = ( )= (r) t = ( =4)1=2(9 =4)1=6 2 ( )r1=4 s ( ) = 1 2 (1 + )2=3 + (1 )2=3 A = ( = ) e"C(rs; )= 3 1 : 24 The GGA takes account of the gradient of the electron density. For systems of slowly varying electron density, the GGA has proved to be an improvement over LDA. For example, GGAs lead to improvement on total energies,[44] atomization energies,[44, 45] energy barriers, and structural energy di erences.[46, 47] 25 2.5 Gaussian Basis Sets When molecular calculations are performed, basis functions are employed to con struct the molecular orbitals, which are expanded as a linear combination of such functions with the coe cients to be determined in the self consistent procedure. A natural choice for basis functions would be Slatertype atomic orbitals decaying exponentially with distance from the nuclei. Boys[48] pointed out that these Slater type orbitals could be approximated as linear combinations of Gaussiantype orbitals (GTOs) instead. Because it is easier to calculate overlap and other integrals involved in the construction of Hamiltonian matrix with Gaussian basis functions, this led to signi cant computational savings. In Cartesian coordinate system, Gaussian primitive basis functions are de ned as products of the integer powers of Cartesian coordinates and a Gaussian function: (r; ; ^n) = xnxynyznze r2 (2.60) and the corresponding normalization constant is given by N( ; ^n) = 2 3=4 (4 )n=2 p (2nx 1)!!(2ny 1)!!(2nz 1)!! (2.61) with n = nx + ny + nz. At present, there are various basis sets composed of Gaussiantype orbitals. The smallest one are called minimal basis sets typically composed of a minimum number of basis functions required to represent all of the electrons on each atom. The most common minimal basis set is STOnG, where n is an integer. The value of n represents the number of Gaussian primitive functions comprising a single basis function. In these basis sets, the same number of Gaussian primitives are used to construct core and valence orbitals. Minimal basis sets are usually used in testing calculations to get the rough idea. 26 Other more complicated basis sets are proposed to improve the quality of the minimal basis set. Based on the observation that the valence orbitals have greater e ects on chemical properties, a splitvalence basis set was proposed such that the number of functions used to describe the valence electrons is doubled but a single function is used for the inner shells. The notation for the splitvalence basis sets is typically XYZG. In this convention, X represents the number of primitive Gaussians comprising each core atomic orbital basis function. The Y and Z indicate that the valence orbitals are composed of two basis functions each, the rst one composed of a linear combination of Y primitive Gaussian functions, the other composed of a linear combination of Z primitive Gaussian functions. For example, in the 321G basis set, three Gaussian functions describe the core orbitals and the valence electrons are also represented by three Gaussians of which two are used for the contracted part and one for the di use part. Other improvement is the addition of polarization functions, indicated by an as terisk, . Two asterisks, , indicate that polarization functions are also added to light atoms (hydrogen and helium). These polarization functions have one more node. For example, the only basis function located on a hydrogen atom in a minimal basis set would be a function approximating the 1s atomic orbital. When polarization is added to this basis set, a pfunction is also added to the basis set. This adds some extra exibility within the basis set, e ectively allowing molecular orbitals involving the hydrogen atoms to be more asymmetric about the hydrogen nucleus. Throughout this work, we use 631G* split valence basis sets to model the twisted armchair graphene nanoribbons in Chapter 3. The 321G split valence set is used in the simulations of silicon nanowires discussed in Chapter 4. And we also use uncontracted basis sets[49] such as 7s3p for C atoms and 3s for H atoms in the study of planar graphene nanoribbons presented in Chapter 3. 27 2.6 Helical Band Structure Methods For a onedimensional system with helical symmetry such as carbon nanotubes, we can de ne a screw operation in terms of a translation l down the z axis in conjunction with a righthanded rotation about the z axis. That is, S r = S(a; ) r = 0 BBBB@ x cos y sin x sin + y cos z + l 1 CCCCA : (2.62) Because the symmetry group generated by the screw operation S is isomorphic with the onedimensional translation group, we may have a generalized Bloch's theorem in which the oneelectron wavefunctions will transform according to Sm i(r; ) = ei m i(r; ): (2.63) The quantity is a dimensionless quantity which is conventionally restricted to < . For the case of pure translation, i.e., = 0, corresponds to a normalized quasimomentum in terms of = kl, where k is the traditional reciprocal wavevector. The oneelectron wavefunctions i are constructed from a linear combination of Bloch functions 'j , which are in turn constructed from a linear combination of nuclearcentered Gaussiantype orbitals j(r) in forms of i(r; ) = X j cji( ) 'j(r; ) (2.64) 'j(r; ) = X m ei m Sm j(r): (2.65) The detailed description and formulism derivation of this helical method can be found in Ref. [50]. 28 2.7 Landauer Approach for Quantum Conductance Calculation The conductance calculations have been extensively discussed in literature. This section is mainly adapted from Refs. [51, 52] with more details carried out. The transport calculations are based on phasecoherent transport of electrons through the device region from the left semiin nite lead to the right semiin nite lead of a linear system. In Landauer approach,[23, 24] the conductance is proportional to the transmission function as follows Gcon = 2e2 h T; (2.66) where both the conductance Gcon and the transmission function T are functions of energy E. T represents the probability that an electron injected at one end of the conductor will emit at the other end. The transmission function T of the device can be expressed in terms of the Green's functions of the conductor and the coupling of the conductor to the leads, T = Tr[LGRRGA]; (2.67) where the advanced Green's function GA is the Hermitian conjugate of the retarded Green's function GR of the conductor describing the dynamics of the electrons inside the conductor, and the fL;Rg accounts for the coupling of the conductor to the leads. To compute the Green's function of the conductor we could start with the retarded Green's function of the whole system including the conductor and two semiin nite leads, (" H)G = I; (2.68) where " = E +i ( is an arbitrarily small positive real number) is a complex energy, I is the identity matrix, and we neglect the superscript in G for simplicity. Since the whole system can be conceptually divided into three distinct regions : a conductor 29 region, a lefthand lead, and a righthand lead, Eq. (2.68) can be expressed in terms of submatrices that correspond to di erent subsystems, 0 BBBB@ (" HL) hLC 0 hy LC (" HC) hCR 0 hy CR (" HR) 1 CCCCA 0 BBBB@ GL GLC GLCR GCL GC GCR GLRC GRC GR 1 CCCCA = I: (2.69) Since we are more concerned with the Green's function of conductor, the equations involving the GC can be written explicitly as follows : (" HL)GLC + hLCGC + 0 = 0 (2.70) hy LCGLC + (" HC)GC + hCRGRC = I (2.71) 0 + hy CRGC + (" HR)GRC = 0 (2.72) From Eq. (2.72) we have GRC = (" HR)1hy CRGC: (2.73) From Eq. (2.71) we get GLC = (" HL)1hLCGC: (2.74) Substituting the expressions for GRC and GLC into Eq. (2.72) leads to hy LC(" HL)1hLCGC + (" HC)GC hCR(" HR)1hy CRGC = I; (2.75) which suggests GC = h (" HC) hy LC(" HL)1hLC hCR(" HR)1hy CR i1 : (2.76) Then we can write the expression of the retarded Green's function of a system as GR = [(" HC) RL RR ]1 (2.77) where Rf L;Rg is the retarded selfenergy terms describing the coupling between the conductor and semiin nite leads, and is given by RL = hy LCgR LhLC; RR = hCRgR Rhy CR (2.78) 30 with gR L and gR R are the Green's functions of the semiin nite leads ("HL)1 and (" HR)1, respectively. The selfenergy term can be viewed as an e ective Hamiltonian term coming from the interaction of the conductor with leads. The coupling matrices fL;Rg can be easily obtained as fL;Rg = i[ Rf L;Rg Af L;Rg]; (2.79) where the advanced selfenergy Af L;Rg is the Hermitian conjugate of the retarded selfenergy Rf L;Rg. The core of the problem lies in the calculation of the Green's functions of the semi in nite leads. The lead Green's function in an orthogonal localizedorbital Hamil tonian can be computed with an e cient principal layer method and we take the righthand side lead as an example. With increasing the distance between two localized orbitals, the overlap and in teraction matrix elements are smaller and smaller. Beyond some length, the elements could be considered as zero. So the semiin nite lead can be viewed as an in nite stack of principal layers with only nonzero nearestneighbor interactions. This corresponds to transforming the original system into a linear chain of principal layers. Within this approach, the Green's function equation for righthand lead can be expressed as 0 BBBBBBBBBBB@ " H00 H01 0 0 Hy 01 " H11 H12 0 0 Hy 12 " H22 H23 0 0 Hy 32 " H33 ... ... ... ... 1 CCCCCCCCCCCA 0 BBBBBBBBBBB@ G00 G01 G02 G03 G10 G11 G12 G13 G20 G21 G22 G23 G30 G31 G32 G33 ... ... ... ... 1 CCCCCCCCCCCA 31 = 0 BBBBBBBBBBB@ I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I ... ... ... ... ... 1 CCCCCCCCCCCA (2.80) which can be expanded explicitly as follows (" H00)G00 = I + H01G10 (" H00)G10 = Hy 01G00 + H01G20 (" H00)Gn0 = Hy 01Gn1;0 + H01Gn+1;0 (2.81) where Hnm and Gnm are the matrix elements of the Hamiltonian and Green's function between the layer orbitals, and we assume that in a bulk system H00 = H11 = and H01 = H12 = . In the iterative method proposed by LopezSancho et al.,[53, 54] this chain can be transformed such than the Green's function of an individual layer can be expressed in terms of Green's function of the preceding (or following) one. This is achieved by introducing transfer matrices T and T, de ned such that G10 = TG00 and G00 = TG10. The transfer matrix can be easily computed from the Hamiltonian matrix elements via an iterative procedure. T and T can be written as T = t0 + ~t0t1 + ~t0~t1t2 + + ~t0~t1~t2 tn T = ~t0 + t0 + t0t1~t2 + + t0t1t2 ~tn where ti and ~ti are de ned via the recursion formulas ti = (I ti1~ti1 ~ti1ti1)1t2i 1 ~ti = (I ti1~ti1 ~ti1ti1)1~t2i 1 32 with initial values t0 = (" H00)1Hy 01 ~t0 = (" H00)1H01 The process is repeated until tn; ~tn with arbitrarily small. In our calculations, we xed at 106. We always order the basis function from left to right in the matrix. Since there is only nonnegligible interaction between the conductor and the surface layer close to the conductor, only the topleft block G00 in Eq. (2.80) is needed. So the nonzero part of surface Green's function of right lead can be written as gR;00 = G00 = [" H00 H01T]1 : (2.82) and similarly the nonzero part of the surface Green's function of left lead can be written as gL;00 = h " H00 Hy 01 T i1 (2.83) In the case of general nonorthogonal orbitals, the above derivations still apply except the following substitutions " H00 ! "S00 H00 H01 ! ("S01 H01) where the matrix S represents the overlap between the localized orbitals. 33 CHAPTER 3 GRAPHENE NANORIBBONS 3.1 Introduction Graphene has attracted a great deal of research interest since its isolation via me chanical exfoliation[1, 2] because graphene exhibits many unusual properties such as massless fermions, minimum quantum conductance, and anomalous integer quantum Hall e ect occurring at halfinteger lling factors.[55] In addition graphene sheets show promise for future nanoelectronic applications, such as ultrafast transistors[56] because of the high carrier mobility.[3] Graphene materials are expected to play an important role in future electronic applications and even replace silicon some day. Twodimensional graphene, however, has a zero band gap with linear energy disper sion near the Fermi level. For practical application in the semiconductor industry a band gap could be induced by making graphene nanoribbons (GNRs) via a variety of methods such as lithographic patterning,[5] chemical vapor deposition,[6] Joule heating,[7] and unzipping the carbon nanotubes (CNTs).[8, 9] Graphene has two major types of high symmetry structures, namely those with edges that have either armchair or zigzag con gurations. Fujita, et al., [57{59] carried out tightbinding studies of these systems using them as model structures to study edge defects in CNTs. It was predicted that in zigzag GNRs there were strongly localized edge states near the Fermi level, suggesting electronic structure and low bias transport very sensitive to di erent edge passivations. The localized edge states were not expected in the armchair GNRs. The electronic energy gaps of armchair GNRs were calculated to have strong dependence on the ribbon widths. When we 34 started the project on GNRs in summer 2006, there were very few simulations done with density functional theory.[60, 61] Most studies were focused on the planar GNRs so far. there were only several investigations addressing the torsional deformation in GNRs. Bets and Yakobson,[62] using Molecular Dynamics Simulations, found that the narrow bare GNRs prefer the twisted structure rather than the planar geometry. Hod and Scuseria[63] reported the density functional theory calculations on twisted nanoribbons of nite length and found the gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is tunable by torsional deformation. In this chapter, we studied the electronic structures of planar zigzag and arm chair graphene nanoribbons. In addition, we investigated how the armchair graphene nanoribbons respond under applied twist about the ribbon axis. 35 3.2 Planar Zigzag Graphene Nanoribbons 3.2.1 Overview Figure 3.1 depicts the model structure of a zigzag graphene nanoribbon with ribbon width N = 6. Here the width N refers to the number of zigzag chains along the transverse direction. In each unit cell there are 2N C atoms and 2 H atoms. The in nite ribbon could be obtained by repeating the unit cell framed by dashed line in Figure 3.1 periodically along two sides with translation length l = 2.46 A. In our calculations we xed the CC bond length at 1.42 A and CH bond length at 1.08 A. We de ne the ribbon extension direction as the z axis perpendicular to the transverse direction. The righthanded rotation of screw operation about the z axis is zero in this case. We investigated ribbons with even values of N from 6  20. We used a 7s3p Gaussian basis for C, a 3s basis for H,[49] and 32 points were evenly sampled over the central Brillouin zone. For the zigzag graphene nanoribbons we carried out the nonspinpolarized (NSP), ferromagnetic (FM), and antiferromagnetic (AFM) calculations. The energetic trend on these states corresponding to di erent spin con gurations was analyzed. We found the antiferromagnetic states are energetically favorable. The zigzag ribbons exhibited strongly localized edge states near the Fermi level. Our results with Gaussian basis sets are consistent with others' results obtained by using tightbinding model and density functional theory simulations with numerical atomic orbital methods. [60, 64, 65] 36 Figure 3.1: Sample zigzag GNRs with ribbon width N = 6. The numbering of the unit cells is shown at the left side of the ribbon with only the unit cell labeled 0 shown in their entirety. Each unit cell is composed of 2N C atoms and 2 H atoms. The 2N C atoms in unit cell 0 are numbered as shown. The index of zigzag chains is indicated at bottom from left to right. 37 3.2.2 Energetic Results Figure 3.2 depicts the total energies divided by 2N and energy di erences between di erent spin con gurations as a function of ribbon width N. The corresponding numerical values are presented in Table 3.1. For the ribbon of same width, the antiferromagnetic state has the lowest energy and is energetically favorable. For example, for ribbon of N = 6, the energy of nonspinpolarized state is 0.20 meV higher than that of ferromagnetic state, which is further stabilized by about 0.032 meV in antiferromagnetic ribbon. As the ribbon width increases, the edges play a less important role and the total energies have more contributions from the bulk C atoms in the middle. So the energy di erence is decreasing with increasing ribbon width. When ribbon width N approaches 20, the energy di erence between FM and AFM states has already reduced to 0.00016 meV, suggesting for even bigger ribbons the preference to one of these two states is eliminated. In order to observe the AFM ground state, low temperature is needed since the thermal energy at room temperature is KBT = 26 meV with KB = 8:617343 105 eV/K and T = 300K, which will be able to excite the system and mix the various spin con gurations. 38 (a) (b) Figure 3.2: (a) Total energy di erences divided by 2N between nonspinpolarized and ferromagnetic states in zigzag GNRs as a function of ribbon width N. (b) Total energy di erences divided by 2N between ferromagnetic and antiferromagnetic states in zigzag GNRs as a function of ribbon width N. Table 3.1: Total energies divided by 2N and energy di erences between di erent spin con gurations of zigzag GNRs with various ribbon widths width N NSP (eV) FM (eV) AFM (eV) NSPFM (meV) FMAFM (meV) 6 37.4358823900 37.4360743104 37.4361059109 0.191920458334 0.0316004916669 8 37.4147775171 37.4149432678 37.4149598679 0.165750693746 0.0166001437520 10 37.4021060503 37.4022471735 37.4022567277 0.141123230001 0.0095541400071 12 37.3936847190 37.3938081844 37.3938136308 0.123465429169 0.0054464166652 14 37.3876676725 37.3877772624 37.3877815346 0.109589875002 0.0042721500009 16 37.3831570254 37.3832537549 37.3832574688 0.096729465625 0.0037139687521 18 37.3796450210 37.3797317062 37.3797335467 0.086685177777 0.0018405499986 20 37.3768360202 37.3769151528 37.3769153107 0.079132679993 0.0001578675040 39 3.2.3 Nonspinpolarized LDF Band Structures In the nonspinpolarized calculations, we assumed the spinup and spindown states are degenerate. The typical calculated band structures are presented in Fig ures 3.3, 3.4, and 3.5, corresponding to the ribbon of widths N = 6, 8, and 10, respectively. They all exhibited nearly at bands crossing the Fermi level. The at region in Brillouin zone is ranging from about = 2 =3 to the zone edge, consis tent with the tightbinding calculations.[57] Because the density of states is inversely proportional to the slope of the energy dispersion, the atness in lowerlying bands generates high density of states around the Fermi level. We will show these states are localized on the ribbon edges by plotting the orbital density in real space. 40 Figure 3.3: Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. 41 Figure 3.4: Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. 42 Figure 3.5: Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. 43 Figure 3.6 depicts the orbital density of zigzag GNRs of width N = 6. The wavevectors are sampled at 0.625 , 0.75 , 0.875 , and in the Brillouin zone. For state at the zone edge with = , the wavefunctions are completely localized at the edge C sites. When is moving away from the zone edge, the orbital density extends and decays towards the ribbon center. If = 0:625 , the states become extended across the transverse direction. Our simulation results are consistent with the previous tightbinding prediction in which the neargap state is completely localized at the edge sites when = , and starts to gradually penetrate into the inner sites as deviates from , reaching the extended state at = 2 =3.[57] Figure 3.6: Calculated HOMO band orbital densities of zigzag GNRs of width N = 6, calculated at wavevectors (a) , (b) 0.875 , (c) 0.75 , and (d) 0.625 , respectively. Calculated orbital densities for the LUMO band are given for wavevectors (e) , (f) 0.875 , (g) 0.75 , and (h) 0.625 , respectively. The isovalue is taken to be 0.008 for all plots. 44 3.2.4 Ferromagnetic LDF Band Structures If we do not apply constrains of spin degeneracy, spinpolarized calculations gave us ferromagnetic states. The typical band structures are shown in Figures 3.7, 3.8, and 3.9, corresponding to the ribbon of widths N = 6, 8, and 10, respectively. In this case the spin degeneracy was lifted. Compared to the nonspinpolarized results, the at bands and the corresponding high density of states around Fermi level were not observed. Instead two peaks in density of states appeared with centers at energy points 0.3 eV above and below the Fermi level. For electron excitation, the spin down states dominate the density of states. For hole excitation, the spinup states contribute more. 45 Figure 3.7: Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. 46 Figure 3.8: Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. 47 Figure 3.9: Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. 48 3.2.5 Antiferromagnetic LDF Band Structures We also carried out the calculations with constraint of applying inversion sym metry about the z axis in the exchange correlation potential. From this, we got the antiferromagnetic ground states. The typical band structures are depicted in Fig ures 3.10, 3.11, and 3.12, corresponding to the ribbon of widths N = 6, 8, and 10, respectively. In Figure 3.13 we depicts the energy gaps as a function of ribbon width N. The corresponding numerical values of gaps are listed in Table 3.2. The ribbons exhibited direct gap nature. However the gaps open up at 0:70 rather than the zone center or edge. When increasing ribbon width, the band gaps are decreasing monotonically. Our results are consistent with that of Son, et al.[64] In their studies, the e ect of electric eld was also addressed. When an electric eld is applied along the transverse direction, spinup and spindown electrons respond in di erent way. The spinup electrons are moving toward the Fermi level, leading to smaller gap even gap closure for big enough eld. The spindown electrons are moving away from the Fermi level and will not play a role in the electron transport within some voltages. The response of electrons under applied electric eld renders the antiferromagnetic graphene nanoribbons a half metal. 49 Figure 3.10: Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. 50 Figure 3.11: Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. 51 Figure 3.12: Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. 52 Figure 3.13: Magnitude of the HOMOLUMO antiferromagnetic gap of the zigzag GNRs as a function of ribbon width N. Table 3.2: Energy gaps of zigzag GNRs with various ribbon widths Ribbon width N Energy gap (eV) Ribbon width N Energy gap (eV) 4 0.5732 14 0.3008 6 0.4972 16 0.2790 8 0.4341 18 0.2618 10 0.3901 20 0.2539 12 0.3404 53 3.2.6 Summary We calculated the electronic structures of zigzag graphene nanoribbons for dif ferent spin con gurations. It was found that there were localized edge states, which make the electronic properties of zigzag GNRs tunable by edge chemistry. [16, 66, 67] The zigzag GNRs exhibited ground states with two ferromagnetic edges antiferomag netically coupled. With transverse electric eld applied, the zigzag GNRs were found to be half metals,[64] which provides a suitable platform for the spintronics. This halfmetallicity is further enhanced by edgeoxidization[68] and edge modi cations with NO2 and CH3 terminations on opposite sides.[69] 54 3.3 Planar Armchair Graphene Nanoribbons 3.3.1 Overview Nakada et al.,[58] with nearestneighbor tightbinding model, predicted that the energy gaps of armchair GNRs are strongly dependent on the ribbon width N. If the ribbon width satis es N = 3m + 2 with m being an integer, the ribbon is metallic of zero gap. Otherwise, the semiconducting energy gaps monotonically decrease with increasing ribbon widths. Our rst principles results with local density functional (LDF) do not agree with the TBM prediction. Especially in our simulations, there is no metallic ribbon. In order to elucidate the discrepancy between TBM and LDF results, we introduced the thirdnearestneighbor interaction as perturbation correc tion to the nearestneighbor TBM and show that the longer range interactions play an important role in describing the band gaps of armchair graphene nanoribbons. 55 3.3.2 Model Structure Figure 3.14 depicts the model structure of armchair graphene nanoribbons with ribbon width N = 7 and 8. Here the width N is referred to as the number of dimer lines along the transverse direction. Depending on the way of constructing the whole in nite ribbon from the basic helical cell as framed by dashed line, the armchair graphene nanoribbons are classi ed into two categories : symmetric and staggered ribbons. For symmetric ribbons as shown in Figure 3.14 (a), the translation length l and rotation angle is 4.26 A and 0 degree, respectively, and in each unit cell there are 2N C atoms and 4 H atoms. For staggered ribbons as shown in Figure 3.14 (b), the translation length l and rotation angle is 2.13 A and 180 degrees, respectively, and each unit cell contains N C atoms and 2 H atoms. In our calculations the CC bond length and CH bond length are xed at 1.42 A and 1.08 A, respectively. 56 (a) (b) Figure 3.14: (a) Symmetric armchair GNRs with ribbon width N = 7. Each unit cell is composed of 2N C atoms and 4 H atoms. (b) Staggered armchair GNRs with ribbon width N = 8. Each unit cell is composed of N C atoms and 2 H atoms. For (a) and (b) the numbering of the unit cells is shown at the left edge of each gure with only the unit cells labeled 0 shown in their entirety. The C atoms in unit cell 0 are numbered as shown. The index of dimer lines is indicated at bottom from left to right. 57 3.3.3 Local Density Functional Results We carried out the local density functional calculations on armchair graphene nanoribbons with various widths ranging from 4 to 24. We used a 7s3p Gaussian basis for C, a 3s basis for H, and 32 points were evenly sampled over the central Brillouin zone. Figure 3.15 depicts the band structures for armchair GNRs of widths N = 7, 8, and 9. and in Figure 3.16 we present the energy gap as a function of ribbon width N using the gap values listed in Table 3.3. In the calculations of staggered armchair GNRs, we were making use of the helical symmetry of rotation. The presented band structures, however, have already been folded to be within the regular Brillouin zone corresponding to the translational operation. First, all ribbons are semiconductors and exhibit direct band gap with HOMO and LUMO occurring at point, irrespective of the ribbon width in our study. On the other hand, the gap decreases with increasing ribbon widths. However, the variations in energy gaps have alternating pattern rather than monotonic decrease. We have relative relation of Egap; 3m+1 > Egap; 3m > Egap; 3m+2 with m being integer. 58 Figure 3.15: Calculated band structures of armchair GNRs of various widths. (a), (b), and (c) correspond to ribbon with N = 7, 8, and 9, respectively. (X) indicates the center (edge) of the rst Brillouin zone. 59 6 9 12 15 18 21 24 N 0 0.5 1 1.5 2 2.5 Egap (eV) LDF mod (N + 1, 3) = 2 LDF mod (N + 1, 3) = 1 LDF mod (N + 1, 3) = 0 Figure 3.16: LDF energy gaps of the armchair GNRs as a function of width N. Table 3.3: Band gaps Egap of armchair GNRs with ribbon width N N Egap(eV) N Egap(eV ) N Egap(eV) N Egap(eV) 4 2.4106 10 1.0734 16 0.6874 22 0.5055 5 0.2982 11 0.1512 17 0.1006 23 0.0752 6 1.0158 12 0.5692 18 0.3953 24 0.3030 7 1.4872 13 0.8375 19 0.5821 8 0.2012 14 0.1210 20 0.0868 9 0.7298 15 0.4669 21 0.3437 60 3.3.4 Consideration of Longer Range Interactions Tightbinding models were frequently used to analyze the electronic properties of semiconductors. With only nearestneighbor interactions considered tightbinding models have been very successful in describing many properties of singlewall carbon nanotubes (SWCNTs).[22, 70{83] The geometries of SWCNTs and GNRs have great similarity because both of them can be obtained from graphene by rolling over a chiral vector or cutting along speci c edges. This similarity, along with the success of TBM in SWCNTs, allow us to expect the applicability of TBM in studying GNRs. However, the LDF results are quite di erent from the previous predictions based on nearestneighbor tightbinding model. In the tightbinding model, if the ribbon width N satis es N = 3m + 2, then the ribbon is metallic with zero gap at point where and bands crossing occurs. In addition, the band gaps of other ribbons are decreasing monotonically with increasing ribbon width. The discrepancy motivated us to think over the validity of nearestneighbor tight binding model in describing the electronic property of armchair GNRs. In the previous calculations, only nearestneighbor interactions are included. The gap openings in armchair GNRs of width N = 3m+ 2 from our rstprinciples calculations, however, indicated that there were some interactions which are important but ignored in the tightbinding description. Since the armchair GNRs can be obtained by cutting graphene along speci c direction, we could apply proper boundary condition on the energy dispersion for 2D graphene, discretize the wavevector along the transverse direction, and nd the energy dependence on the 1D wavevector in armchair GNRs. However, we are more concerned with the band gap instead of the energy values in whole Brillouin zone and will take a di erent method. The underlying procedure is described below. 61 The ribbons are direct band gap materials with HOMO and LUMO states occur ring at k = 0 (see Figure 3.15). Therefore, to study the gap dependence on ribbon width we need only to consider the eigenstates of nearestneighbor Hamiltonian ^H 0 at k = 0. In the following discussion, we make use of the regular Bloch theorem to construct the wavefunctions based on the translational symmetry down the ribbon axis. So each unit cell contains two zigzag chains for symmetric and staggered arm chair GNRs. We can apply the the boundary condition with vanishing nodes at the ends of zigzag chains. Then the eigenstates of ^H 0 for the armchair GNRs with ribbon width N are given by j p >= 1 p Ncell 1 p N + 1 X l XN n=1 sin( pn N + 1 )(jn; l > +sjn + N; l >); (3.1) where 1 p N, jn; l > denotes the jpz > orbital associated with the nth C atom in the unit cell labeled by l in the notation of Figure 3.14, and the s = 1 (1) indicates the bonding (antibonding) state between the wavefunctions belonging to neighboring zigzag chains. The calculated energy gap with nearestneighbor interaction V1 only is Egap = 2V1[2 cos( p N + 1 ) + s] (3.2) where we have : if mod(N+1,3) = 0, then p = (N +1)=3 and s = 1; if mod(N+1,3) = 1, then p = (2N + 3)=3, and s = 1; and, if mod(N+1,3) = 2, then p = (N + 2)=3, and s = 1. Then the band gap strongly depends on the ribbon width. The system is metallic if N = 3q 1, and is semiconducting otherwise, see Figure 3.17 in which V1 = 2:6 eV is assumed. Except the nearestneighbor interaction, longer range interactions between C atoms are neglected and may modify the band structure results of TBM. These longer range interactions can be incorporated into the TBM by adding ^U V2 and ^U V3 to ^H 0, the Hamiltonian in TBM, where ^U V2 (^U V3) includes all secondnearestneighbor (thirdnearestneighbor) interactions V2 (V3) between C atoms in the armchair GNRs. Because V2 and V3 are small in magnitude compared to V1, ^U = ^U V2 + ^U V3 can be 62 6 9 12 15 18 21 24 N 0 0.5 1 1.5 2 2.5 Egap(eV) TBM mod (N + 1, 3) = 2 TBM mod (N + 1, 3) = 1 TBM mod (N + 1, 3) = 0 Figure 3.17: Magnitude of the HOMOLUMO gap of the N armchair GNRs from the TBM calculations. treated as perturbation to ^H 0. Also, because ^U does not break the translational symmetry of the ribbon, it couples only states with the same k, which allows us to consider the states at k = 0 only. We carried out the timeindependent perturbation calculations to rst order with the V2 and V3 as perturbation terms to V1 and found the new expression for Egap given by Egap = 2V1[2 cos( p N + 1 ) + s] + s 2V3 N + 1 [3 + N + 2N cos( 2p N + 1 )] (3.3) where we have : if mod(N+1,3) = 0, then p = (N +1)=3 and s = 1; if mod(N+1,3) = 1, then p = (2N + 3)=3, and s = 1; and, if mod(N+1,3) = 2, then p = (N + 2)=3, and s = 1. The rst term on the righthand side of Eq. (3.3) is the gap arising from ^H 0. The second term gives the lowest order correction because of ^U . Although ^U V2 contributes to ^U, V2 does not appear in Eq. (3.3) because to rst order it shifts the HOMO and LUMO levels by the same amount 3NV2=(N + 1). If N + 1 = 3m, the rst term on the righthand side of Eq. (3.3) vanishes leaving 63 only the second term which reduces to Egap = 6V3 N + 1 : (3.4) Eq. (3.4) can be used to leastsquares t to the LDF data for N +1 = 3q to determine V3 and then with V3 xed at this value V1 can be determined by using Eq. (3.3) to leastsquares t to the LDF data for N +1 6= 3q. Implementing this procedure yields the physically reasonable results: V1 = 3:2 eV and V3 = 0:3 eV, which provide an excellent t to the LDF results for all N as shown in Figure 3.18. 6 9 12 15 18 21 24 N 0 0.5 1 1.5 2 2.5 Egap (eV) LDF mod (N + 1, 3) = 2 LDF mod (N + 1, 3) = 1 LDF mod (N + 1, 3) = 0 TBM mod (N + 1, 3) = 2 TBM mod (N + 1, 3) = 1 TBM mod (N + 1, 3) = 0 Figure 3.18: LDF results for Egap compared to the corresponding tightbinding model (TBM) results obtained from Eq. (3.3) with V1 = 3:2 eV and V3 = 0:3 eV. 64 3.3.5 Summary In our LDF simulations we xed all CC bond lengths at 1.42 A and found that the armchair GNRs can be separated into three families according to the band gap dependence on ribbon width. Then we proposed an explanation by using TBM with perturbation correction from the second and thirdnearestneighbor interactions. The nonzero semiconducting gaps in GNRs of width N = 3q + 2 and three families of armchair GNRs were also observed by Son, et al.[84] In their calculations, they carried out geometry optimization and found that the CC distances on edges decrease leading to 12% increase of the hopping integrals between pz orbitals of edge C atoms. After taking matrix elements change into account, they also got a TBM results which agree with their LDA calculations. Gunlycke, et al.[85] found that both edge contraction and thirdnearestneighbor interactions play an important role in producing good tightbinding band structures which t the rstprinciples results of optimized GNRs. Experimentally, great progress has been made in preparing GNRs with widths less than 10 nm,[86, 87] which makes our theoretical calculation of more fundamental and practical signi cance. 65 3.4 Twist E ect in Electronic Properties of Armchair Graphene Nanoribbons 3.4.1 Overview In the study of planar armchair GNRs, we show that the band gap could be engineered by controlling the ribbon width. The band gaps of nanostructures also could be tuned by other ways such as edge passivation, doping, and strain. Strain is a powerful approach to modify the electronic structure of nanomate rial to improve device performance. The carrier mobility of silicon nanowires could be signi cantly enhanced by applied strain.[88] Firstprinciples calculations have shown that for smaller diameters the band structures and gaps can be signi cantly modi ed.[89, 90] The strain e ect in graphene also attracted much attention. Ni, et al.[91] predicted a bandgap opening of 300 meV for graphene under 1% uniaxial tensile strain. Teague, et al.[92] found that the strain could modulate the local conductance of graphene. Lu and Guo,[93] using a tightbinding model, investigated armchair GNRs and found uniaxial weak strain changes the band gap in a linear fashion and large strain results in periodic oscillation of the band gap. Sun, et al.,[94] using GGA exchange correlation potentials, predicted that the similar oscillating behavior with linear change in each period. For the outofplane strain, Hod and Scuseria[63] studied the torsional deformation e ect in graphene nanoribbons of nite length. With traditional band structure codes using translational symmetry, it is hard even almost impossible to study the torsional deformation for in nitelong ribbon structure because of the huge unit cell and hence the Hamiltonian matrix of enormous size. Our HENS parallax code takes advantage of the helical symmetry and is well suited to the calculation of electronic structures 66 of the twisted quasionedimensional structures. In this section we investigated the band gaps of armchair GNRs as a function of twist angle. We found that the electronic properties of armchair GNRs are highly sensitive to the applied twist, indicating the potential application of armchair GNRs as building blocks in nanoelectromechanical devices. 67 3.4.2 Computational Approach We investigated the armchair GNRs of widths N = 8; 10; 12; :::; 18 in this study. The representative structure of the armchair GNRs is shown in Figure 3.19. The whole ribbon can be constructed with a unit cell consisting of one zigzag chain as framed by dashed lines in Figure 3.19(a) and a screw operation acting on that unit cell, where the screw operation combines a translation l down the ribbon axis with a righthanded rotation about that same axis. In our calculations we xed l at 2.13 A. For the planar structure, the is 180 . For the twisted armchair GNRs as shown in Figure 3.19(b), the torsional deformation of armchair GNRs can be represented by the relative rotation angle , de ned as = 180 . The planar ribbons have symmetry operations corresponding to a re ection perpendicular to the ribbon axis and passing through the center of CC bonds parallel to the ribbon axis. When applying the twist, the re ection symmetry of armchair GNRs is broken and only C2 symmetry is left. All edge C atoms are passivated by H atoms to remove the dangling bands near the Fermi level. Following conventional notation, the ribbon width of armchair GNRs is de ned as the number of dimer lines along the ribbon forming the width as in the preceeding section. For example, the ribbon width in Figure 3.19 is 8. The LDF band structures were obtained by using 32 discrete points evenly sampled in the central Brillouin zone and a 631G* basis set. [95] Although the LDF is well known to underestimate the band gap, we are primarily concerned with the relative change of the energy gaps under di erent twists. 68 (a) (b) Figure 3.19: (a) Armchair GNR with width N = 8. Each helical unit cell is composed of one zigzag chain terminated by H atoms. The N C atoms in a helical cell are numbered as shown. The bonds are labeled as a1a11. (b) Axial view of 10 twisted armchair GNR with width N = 8. 69 3.4.3 Band Structures and Neargap Wavefunctions We depict the electronic structures for armchair GNRs with width N = 8 under a series of twist angles in Figure 3.20, in which (a), (b), (c), and (d) correspond to twist angle equal to 0 , 8:5 , 9:0 , and 10 , respectively. In this ribbon, the top of the valence band and the bottom of the conduction band occur at the center of the one dimensional Brillouin zone, resulting in a direct band gap at the zone center point. The band gap for ribbon of width N = 8 is 0.40 eV for = 0. When we apply on the ribbon torsional strain about the ribbon axis, the HOMO is raised while the LUMO is lowered if < 9:0 . therefore the band gap is reduced in this stage. In general two bands cannot cross because they belong to the same irreducible representation for wave vectors not at the center or edge of the Brillouin zone. However, gap closure (accidental degeneracy) might take place at a point of high symmetry in the Brillouin zone such as the zone center or end. After the critical twist angle, the HOMO and LUMO do the inverse way, opening the band gap. In addition to the band gap change, we also observed the twist e ect on the characteristics of the bonding or antibonding states in HOMO and LUMO. Figure 3.21 depicts the wave functions of HOMO and LUMO under four di erent twist angles corresponding to that of Figure 3.20. According to the bonding states between the pz orbitals of neighboring zigzag chains, there are two kinds of bonds: bonding states (BS) and antibonding states (AS). When the twist angle is smaller than 9:0 , the HOMO is BS state while LUMO is AS state. When the twist angle gets bigger, the characteristics of HOMO and LUMO exchanged. The energy of BS and AS states respond di erently to the geometry change. Increasing bond length will lower the energy of AS state and raise the energy of BS state. 70 Figure 3.20: Electronic band structure of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . 71 Figure 3.21: The LUMO and HOMO wave function of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . The ribbons are vertically oriented. 72 3.4.4 Structure Parameters In our simulations we carried out geometry optimization for all structures. Fig ure 3.22 depicts the variation of bond lengths in armchair GNRs of width N = 8 as a function of the twist angle. Because of symmetry in this ribbon, there are only eight nonequivalent bonds labeled by a1, a4, a7, a10 along the ribbon extension direction and by a2, a3, a5, a6 along the transverse direction. By symmetry the bonds a8, a9, and a11 are equivalent to a5, a3, and a2, respectively. Depending on the position, di erent bonds in a di erent way respond to the torsional deformation. The bonds along the transverse direction change a little bit with maximum change about 0.01 A. The bonds connecting neighboring zigzag chains exhibit bigger change. Generally the change is increasing as bonds location move towards the ribbon edges. But the biggest bond length change of 0.045 A occurs on a10 instead of a1 which is furthest from the ribbon center. As the twist angle approaches 9 , some bonds exhibit abrupt change, which corresponds to the transition between the antibonding and bonding HOMOs. In other ribbons of di erent widths, we observed similar bond length change as a function of twist angle and the bond equivalent to a10 has the biggest change. The bonds underlying bigger change play more important role in determining the shift of one electron energy of HOMO and LUMO. So only from the bonding (antibonding) states between interchain bonds we could acquire some idea about how the twist a ects the electronic band gap of armchair GNRs. Tightbinding model could be used to carry out this analysis. 73 (a) (b) Figure 3.22: CC bond lengths of the twisted armchair GNRs with width N = 8 as a function of twist angle . 74 3.4.5 Tightbinding Model We start with the wavefunction as introduced in previous section, 1 p Ncell(N + 1) X l XN n=1 sin np N + 1 (jn; li + sjn + N; li) where 1 p N, and jn; li denotes the pz orbital associated with nth C atom in the unit cell labeled by l in the notation of Figure 3.19. And the eigen energy at can be expressed as follows: E(s; p) = 2V1[cos( p N + 1 ) + s]; where s = 1 indicates the BS and s = 1 indicates the AS. In this nearestneighbor tightbinding model, the ribbon of N = 3m + 2 is metallic with HOMO and LUMO touching each other. We cannot distinguish the HOMO and LUMO. In this case, we introduced the thirdnearestneighbor interactions as perturbation. The perturbation shifted the energy by 3=(N + 1)V3 upward for LUMO and the same amount down ward for HOMO. The third nearestneighbor interactions are required to open the gap to give results consistent with our LDA results. After introducing these interac tions, we determine that LUMO is at p = (N + 1)=3 and antibonding, HOMO is at 2(N + 1)=3 and bonding. Because in some ribbons, the second HOMO (HOMO1) and second LUMO (LUMO+1) are moving towards the Fermi level under twist and will play a role in determining the band gap beyond certain twist angle. We also analyzed the HOMO1 and LUMO+1 states and present the results in Table 3.4. 75 Family HOMO LUMO HOMO1 LUMO+1 N = 3m # +" +" # N = 3m + 1 +" # # +" N = 3m + 2 +" # # +" Table 3.4: Binding states of edge C pz orbitals HOMO and LUMO. " (#) indicates the energy point shifts upward (downward) with increasing twist angle. + () indicates the bonding (antibonding) states. 76 3.4.6 Band Gap Change Under Twist In Figure 3.23 and Table 3.5 we present the local density functional results for the variations of energy gaps of three family structures with di erent widths (N = 3m, 3m+1, 3m+2, where m = 2; 3; 4) as a function of the twist angle . In combination with the tightbinding model in previous section, we will discuss in details how the twist a ects the energy gaps in armchair graphene nanoribbons. Depending on the bonding (antibonding) states of near Fermi level bands, the response of the armchair GNRs to the twist can be classi ed into two categories. One is of width N = 3m, and the other two belong to the second family. For N = 3m, the energy gap will go through three stages with increasing twist. In the rst stage, the bonding LUMO goes up and antibonding HOMO goes down, resulting in the raising of band gap. The HOMO1 is raising and LUMO+1 is lowering. After certain twist angle, HOMO1 and LUMO+1 will take over to determine the band gap. This begins the second stage, in which the band gap is reduced with increasing twist angle. When the twist is increased further, the HOMO1 and LUMO+1 may touch and exchange bonding characteristics and then follow their previous way. This will start the third stage, in which the band gap is increased again. For N = 12 we only observed the rst stage. For N = 18, we observed all these three stages because the reduced quantum con nement and bigger change in length of bonds parallel to the ribbon axis resulted from the bigger ribbon width. For N = 3m + 1 and N = 3m + 2, the HOMO1 and LUMO+1 are moving away from the Fermi level and will not play a role in determining the band gap. Within the twist angle studied, we only have to consider the HOMO and LUMO orbitals. Just as we have seen in armchair GNRs with width N = 8, the band gap will have two stages. The rst stage is decreasing the band gap and the second would be increasing 77 the gap. In N = 10, we only have the rst stage because of the relatively bigger band gap and not big enough ribbon width. Generally the turning point is moving toward smaller twist angle with increasing ribbon width for same family. At each turning point corresponding to the minimum band gap, the gap closure may occur because of the additional symmetry at the zone center or edge. We can not have real band crossing because the eigen vectors of bands belong to the same irreducible representation group. Our bonding and antibonding analysis applies for other strain e ect in the arm chair GNRs as well. Our tight binding analysis are qualitatively consistent with the previous calculations on uniaxial strain e ect on armchair GNRs by others. But the linear relation by TB model[93] or density functional theory calculations[94] is not observed in our calculations. This may be due to the hybridization between the s, px, py, pz orbitals coming from the twist between neighboring zigzag chains. 78 Figure 3.23: Variation of energy gaps of armchair GNRs as a function of twist angle . (a), (b), and (c) represent the armchair GNRs of width N = 3m, 3m + 1, and 3m + 2, respectively. 79 Table 3.5: Band gaps (eV) as a function of twist angle (Deg) in armchair graphene nanoribbons with various widths N Twist angle N=8 N=10 N=12 N=14 N=16 N=18 0.0 0.3958 1.2364 0.4648 0.2830 0.8166 0.2795 0.5 0.3904 1.2382 0.4677 0.2800 0.8170 0.2850 1.0 0.3842 1.2281 0.4732 0.2717 0.8159 0.2949 1.5 0.3860 1.2202 0.4816 0.2613 0.7983 0.3147 2.0 0.3727 1.2125 0.4951 0.2406 0.7764 0.3415 2.5 0.3632 1.1979 0.5109 0.2268 0.7486 0.3744 3.0 0.3579 1.1869 0.5276 0.1972 0.7125 0.4146 3.5 0.3456 1.1658 0.5498 0.1651 0.6714 0.4619 4.0 0.3313 1.1440 0.5738 0.1276 0.6239 0.5165 4.5 0.3116 1.1226 0.6051 0.0868 0.5839 0.5737 5.0 0.2909 1.0954 0.6329 0.0217 0.5215 0.6404 5.5 0.2739 1.0688 0.6642 0.0697 0.4568 0.7115 6.0 0.2541 1.0399 0.7004 0.1190 0.3907 0.7197 6.5 0.2309 1.0085 0.7374 0.1732 0.3200 0.6293 7.0 0.2062 0.9763 0.7749 0.2302 0.2435 0.5341 7.5 0.1821 0.9423 0.8180 0.2885 0.1625 0.4341 8.0 0.1541 0.9063 0.8610 0.3533 0.0780 0.3273 8.5 0.1274 0.8689 0.9070 0.4214 0.0761 0.2141 9.0 0.1990 0.8293 0.9532 0.4891 0.1746 0.0994 9.5 0.0505 0.7895 1.0027 0.5673 0.2764 0.0702 10.0 0.0815 0.7492 1.0560 0.6476 0.3847 0.1582 80 3.5 Summary and Conclusions In this chapter we studied two kinds of GNRs: zigzag and armchair GNRs. We found that the zigzag GNRs have magnetically ordered insulating ground states that are ferromagnetically coupled along each edge and antiferromagnetically coupled across the edges. And the AFM band gap decreases with increasing ribbon width. The electronic states near Fermi level are localized at edge carbon sites, suggesting the edge sensitivity to passivating functional groups and therefore zigzag GNRs could be used in chemical sensor applications. Our LDF results with Gaussian basis sets agree with previous calculation by using plane wave basis, numerical atomic orbitals, and also are consistent with the tightbinding calculations on the zigzag GNRs. armchair GNRs do not have spinpolarized ground state and localized edge states. The band gaps of armchair GNRs are strongly dependent on the ribbon width. The complex alternating relation between gaps and ribbon width is di erent from the nearestneighbor tight binding prediction. The introduction of thirdnearestneighbor interaction across the hexagons resolved this discrepancy. Besides the planar GNRs, we investigated the twisted armchair GNRs and found the response to the applied torsional deformation could be classi ed into two categories. The tunable band gap upon applied twist is very useful in nanomechanical devices. It is wellknown that the LDF underestimated the band gaps. Yang, et al.[96] pre sented calculations of the quasiparticle energies and band gaps of planar graphene nanoribbons carried out using a rstprinciples manyelectron Green's function ap proach within the GW[97{100] approximation. The selfenergy corrections are di er ent for ribbons of di erent widths because of the screening e ect of various magni tudes. However, the characteristics of three families in band gaps remain the same in armchair graphene nanoribbons. 81 CHAPTER 4 SURFACE PASSIVATION EFFECTS IN SILICON NANOWIRES 4.1 Overview In recent years Si nanowires (SiNWs) have attracted much attention because of their potential applications in electronic, thermoelectric, and optoelectronic devices. They can be used as building blocks of nanoscale electrical devices such as eld e ect transistors (FETs).[101{104] They are capable of in ating 4 times their normal size when absorbing lithium ions, which enables the new battery to hold 10 times the charge of existing lithiumion batteries.[10] SiNWs can be used as thermoelectric materials because of their low thermal but high electrical conductivity.[12, 13] SiNWs also have potential application as solar cells to convert light into electricity. [11, 105{ 110] The broad range of potential applications of SiNWs makes it critically important to understand how to tailor their electronic properties. The electronic properties of SiNWs can be modi ed by varying a range of structural properties such as surface passivation,[101] doping,[102] and how the nanowires are mechanically processed.[111] Because of the enhanced surfacetovolume ratio of the SiNWs compared to bulk silicon, the surface of the SiNWs is of special status and the electronic properties of the SiNWs are strongly dependent on the characteristics of the SiNW surfaces. Cui, et al.[101] showed that exposure to 4nitrophenyloctadecanoate or tetraethylammo nium bromide can improve the conductance and ono ratios in oxidized SiNWbased FETs. Haick, et al.[112] achieved the methyl passivation on SiNWs via a twostep chlorination/alkylation method and showed that the SiNWs have enhanced air sta 82 bility and high hole mobility. Although the question of how the dangling bonds are passivated at the silicon nanowire surface is of great importance, most of the theoretical studies have been carried out on Hterminated[99, 113{119] and bare SiNWs,[120{122] focusing on quan tum con nement e ects, orientation, and crosssection dependence of the electronic properties of SiNWs, doping e ects, and surface reconstructions. There are very few investigations focusing on the surface e ects. For examples, Blase, et al.[123] calcu lated the e ects of a single alkyl chain on the zerobias conductance of SiNWs and found that the Landauer conductance near Fermi level is almost not a ected. Nolan, et al.,[124] using a combination of plane wave basis sets and pseudopotentials, ob served that passivation with OH or NH2 reduces the band gap of h100i SiNWs. Aradi, et al.[125] reported that the OH passivation also can induce the band gap reduction of h110i SiNWs. Leu, et al.[126] calculated the band structure of SiNWs with halogens surface substituents passivating the SiNW surfaces, and found that passivation with more weakly interacting surface species appears to lead to more marked reductions in the SiNW band gaps. The band gap reduction comes from the weakly interacting surface species which can not pull surface states of bare SiNWs out of midgap. In our studies we choose the silicon nanowires along h110i and h100i as objects of study based on the experimental observation on the relations between the wire diam eter and orientations. When synthesizing silicon nanowires using gold nanocluster catalyzed 1D growth method, Wu, et al.,[127] found that the nanowires of diameter less than 10 nm prefer the h110i directions. Most bigger nanowires were h112i and h111ioriented. Schmidt, et al.,[128] using epitaxial vaporliquidsolid (VLS) growth method, observed similar orientation dependence on the diameter. To my knowl edge, no h100i SiNWs were obtained with the nontemplated methods. Because only small silicon nanowires can be treated with the density functional theory, most com putational simulations are focused on h110i SiNWs such that the calculation results 83 could be compared with experimental measurements. The h100ioriented nanowires are highly desirable since the conventional CMOS microelectronics are built on (100) silicon wafer. The h100i SiNWs can be used to fabricate the vertical eld e ect tran sistors on the (100) wafer. Researchers, with anodic alumina membranes template directed growth technique, have obtained h100i silicon nanowires of diameter ranging from 60 nm to 200 nm.[129{131] In this chapter we rst present the rstprinciples study of the electronic band structures of SiNWs oriented along h100i and h110i directions with the surfaces pas sivated by hydrogen, hydroxyl, and methyl substituent groups. For convenience, sometimes we simply refer to the SiNWs of these three di erent surface passivations as HSiNWs, OHSiNWs, and CH3SiNWs, respectively. We observed the reduced band gaps with increasing wire diameter irrespective of passivating groups, mani festing the quantum con nement e ects. Band gap reductions in OHSiNWs with reference to HSiNWs are also observed in our simulations, in agreement with others' results.[124, 125] CH3 not only reduces the band gap but also leads to indirect band gaps for all h100i SiNWs studied. The passivations on surface not only can change the band gap nature and gap value, but also have e ect on the e ective mass of carriers in SiNWs. Our simulations suggest that passivation with CH3 surface substituents sub stantially increases the electron e ective mass for the h100i wires, while h110i SiNWs have small electron and hole e ective masses for all three passivations studied. The change on electronic band structures by surface passivations motivated us to study their e ect on the transport conductance in SiNWs. As a preliminary study, we investigated the electron conductance of a h110i SiNW having a defect region with hydroxyl passivation rather than hydrogen passivation. 84 4.2 Computational Approaches The Si nanowires studied herein were initially constructed from silicon in the di amond structure using the bulk lattice constant, a = 5.43 A. All h100i nanowires were constructed with square crosssections and f110g faces, and all h110i nanowires were constructed with diamond crosssections and f111g faces. We carried out geom etry optimization for all structures. The optimized structure for the h100i and h110i SiNWs are depicted in Figure 4.1. The h100i wires calculated have a 4fold screw symmetry about the wire axis with translation l = a=4 and screw angle = =2, and the h110i wires have a 2fold screw symmetry about the wire axis with translation l = p 2a=4 and screw angle = . The di erent crosssections of the h100i and h110i oriented SiNWs make a direct comparison of the radial scale of the wires somewhat arbitrary. We have chosen to de ne an e ective nanowire diameter, d, given by the expression d = p a3N=2 l, where a is the lattice constant of bulk crystalline silicon, N is the number of Si atoms in the unit cell, and l is the length of the unit cell along the nanowire axis. This produces an e ective diameter that scales as the square root of the crosssectional area. In the localdensity functional (LDF) calculations, we evenly sampled 64 discrete points over the central Brillouin zone and used the 321G basis set. [132{134] Since our code takes advantage of the helical symmetry, the calculated band structures look quit di erent from those obtained by using other codes, such as VASP, SIESTA, GAUSSIAN, although they are essentially the same. To be able to compare with results reported by others using translational symmetry, all band structures shown in this chapter have been folded to be within the traditional translational Brillouin zone. Although the LDF is wellknown to underestimate the band gap, we are primarily 85 concerned with the relative change of the energy gaps here. We also carried out some testing calculations with PBE functional on the silicon nanowires of small diameter. We numerically calculated the e ective mass of charge carriers in silicon nanowires using the standard expression from solidstate theory m = ~2=(d2E=dk2), which is equal to ~2=(l2d2E=d 2) if we use helical symmetry, where l is the translation length along wire axis and ~ is the reduced Planck's constant. 86 Figure 4.1: Crosssection views of the studied HSiNWs. h100i SiNWs of diameter (a) 0.43 nm, (b) 0.87 nm, (c) 1.3 nm, and (d) 1.73 nm, respectively. h110i SiNWs of diameter (e) 0.73 nm, (f) 1.09 nm, and (g) 1.46 nm, respectively. The golden yellow and silver balls represent Si and H atoms, respectively. In OHSiNWs and CH3SiNWs, all H atoms are replaced with OH and CH3 groups, respectively. 87 Figure 4.2: Band gaps as a function of diameter with hydrogen passivation for silicon nanowires along h100i. The square and blue line indicate the LDF results and the uptriangle and red line indicate the PBE results. For Hterminated h100i SiNWs, we compared the calculated band gaps with LDF and PBE functionals. Figure 4.2 depicts the band gaps as a function of diameter. The consideration of electron density gradient in PBE functional shifted up the LDF band gaps. The variation in the shifted amount, however, is small with maximum 0.05 eV. The LDF functional is su cient to capture the quantum con nement in this nanomaterial and also have advantage of computional cost with comparison to PBE. The LDF was utilized for the study of silicon nanowires in the following discussions. 88 4.3 Structure Parameters The surface passivating groups can a ect the structure parameters of silicon nanowires. The calculated SiSi bond lengths as a function of the distance from the bond center to the wire axis with di erent passivations are depicted in Figure 4.3. After geometry optimization the Hterminated h100i SiNWs exhibited small varia tions in SiSi bond length from 2.37 A to 2.39 A, compared to the equilibrium value of 2.35 A in the bulk silicon. In the OHpassivated h100i wire of diameter d = 1.73 nm, the SiSi distances have larger variations from 2.36 A to 2.49 A. The SiSi dis tances are ranging from 2.34 A to 2.41 A in the CH3terminated h100i wire of same diameter. The variations of SiSi bond length in the core part are smaller than those near the surface. In the h110i SiNWs, we found that the SiSi bond lengths of HSiNWs are similar to those of h100i HSiNWs. The passivation with with OH or CH3 induces smaller variations in SiSi distances because the distances between surfacepassivating groups on f111g faces of h110i wires are larger than those on f110g faces of h100i wires. The CH3 groups even give rise to bond length variations comparable to that of the Hpassivated wires along h110i direction. 89 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 SiSI (Å) H H OH OH CH3 CH3 <100> 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 <110> 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 SiSi (Å) 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 0 1 2 3 4 5 6 7 8 9 10 Radius (Å) 2.2 2.3 2.4 2.5 2.6 SiSi (Å) 0 1 2 3 4 5 6 7 8 9 10 Radius (Å) 2.2 2.3 2.4 2.5 2.6 Figure 4.3: SiSi bond length as a function of the distance from bond center to wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm passivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. 90 4.4 Mulliken Population Analysis We also did a Mulliken population analysis to see how the functional groups a ect the charge distribution in the silicon nanowires. The calculated net charge of atoms as a function of the distance from atom to wire axis are presented in Figure 4.4. For elements Si, H, C, and O, the electronegativity (EN) is 1.9, 2.1, 2.5, and 3.5, respectively.[135] So C and O could pull more electrons away from the neighboring Si atoms than H and weaken the bond length of SiSi. Although the EN of C is smaller than that of O, there are three H atoms attached to each C atom such that C atoms have more places to get charges. This results in that the charge is about 0.9 for C and is about 0.55 for O in the h100i wires of diameter d = 1.73 nm. The similar values were observed in the h110i wires of diameter d = 1.46 nm. The net charge of C and O are almost independent of the position and the wire orientation. The independence of position makes the corner Si atoms donate more charges to C or O and leave itself more positively charged because there are two groups attached to the corner Si. We found that the surface Si atoms are positive and the innermost Si atoms are nearly neutral. In the middle part, the net charge of Si atoms exhibited oscillations instead of gradual transition. 91 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 Net charge (e) H OH CH3 H OH CH3 <100> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 <110> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 Net charge (e) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Radius (Å) 1 0.5 0 0.5 1 Net charge (e) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Radius (Å) 1 0.5 0 0.5 1 Figure 4.4: Net charge as a function of distance from the atom to the wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm passivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. Si, H, O, and C are indicated by +, , , and , respectively. 92 4.5 Band Gap Change All of the passivations (H, OH, CH3) studied in the present paper give indirect band gap for the h100i nanowire of diameter d = 0.43 nm. CH3 even gives indirect band gap for all h100i nanowires investigated. Figure 4.5(a) shows the h100i Si nanowire band gaps as a function of diameter. With increasing diameter the band gap is decreasing because of the reduced quantum con nement. Because surface atoms have smaller and smaller percentage when increasing diameter from 0.43 nm to 1.73 nm, surface e ect becomes weak and energy gap di erence of di erent passivations is reduced from 2.75 eV to 0.64 eV for OHSiNWs and from 1.72 eV to 0.54 eV for CH3 SiNWs with HSiNWs as the reference. We observed similar trend in band gap change of h110i SiNWs as shown in Figure 4.5(b). But di erent from the h100i SiNWs, the h110i wires have direct band gaps independent of the surface substituents. The band gap values corresponding to Figures 4.5(a) and 4.5(b) are listed in Tables 4.1 and 4.2. 93 (a) (b) Figure 4.5: Band gaps as a function of diameter with various surface passivations for silicon nanowires along (a) h100i, (b) h110i. 94 Table 4.1: Band gaps (eV) as a function of diameter for h100i silicon nanowires with various passivations. X represents H, OH, or CH3. Stoichiometry Diameter (nm) H OH CH3 Si2X2 0.43 5.61 2.86 3.89 Si4X4 0.87 3.38 1.75 2.42 Si9X6 1.30 2.38 1.48 1.91 Si16X8 1.73 1.89 1.25 1.35 Table 4.2: Band gaps (eV) as a function of diameter for h110i silicon nanowires with various passivations. X represents H, OH, or CH3. Stoichiometry Diameter (nm) H OH CH3 Si4X4 0.73 2.38 1.01 2.07 Si9X6 1.09 1.79 0.37 1.65 Si16X8 1.46 1.49 0.30 1.39 95 4.6 Electronic Structures We show the band structures of h100i silicon nanowires of diameter d = 1.73 nm with di erent surface passivations in Figure 4.6 and in Figure 4.7 we present the band structures of h110i silicon nanowires of diameter d = 1.46 nm. The Fermi levels in HSiNWs are used as the reference and taken to be zero. Both OH and CH3 change the band structure a lot and OH gives rise to more reduction of band gap. In addition to the band gap change, the Fermi levels also vary with the passivations. The CH3 groups shift up the Fermi levels in h100i and h110i wires compared to H SiNWs. The OH groups induce lower Fermi level in h100i SiNWs and higher Fermi level in h110i SiNWs. It is known that the Fermi levels will be aligned via charge redistribution when systems of di erent Fermi levels are put together. If the silicon nanowires are selectively passivated such that one segment is covered with H atoms and neighboring segment is passivated by CH3 groups. Then the electron in CH3 terminated part will move to Hterminated part and hole will move in the inverse way to align the Fermi levels and lower the total energy of the system. The ability of separating charge carriers in selectively passivated SiNWs could nd potential application in photovoltaic devices. Growth techniques for semiconductor nanowires have developed rapidly in recent years. Not only can the diameter and direction be controlled during growth, but nanowires can also be selectively functionalized. This has been achieved by vari ous e orts such as electrochemical methods,[136] localized nanoscale Joule heating, [137] and adsorption and removal of selfassembled monolayers of the polymer (3 mercaptopropyl)trimethoxysilane.[138] 96 Figure 4.6: Electronic band structures of h100i SiNWs (diameter d = 1.73 nm) pas sivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. 97 Figure 4.7: Electronic band structures of h110i SiNWs (diameter d = 1.46 nm) pas sivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. 98 4.7 Neargap States The neargap states are mostly responsible for the transport and optical properties of the system. In order to see how the surface passivations will a ect those proper ties in detail, we plot the orbital density of the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) for the h100i SiNWs of diameter d = 1.73 nm in Figure 4.8 and for the h110i SiNWs of diameter d = 1.46 nm in Figure 4.9 with sub gures (a), (b), and (c) corresponding to the passivating substituents H, OH, and CH3, respectively. First look at the h100i SiNWs. For the HSiNWs, both HOMO and LUMO are concentrated in the interior of the nanowire. Upon OH passivation HOMO goes to the corner O atoms while LUMO leaves the central region and corners blank. In the CH3SiNWs HOMO and LUMO spread out over the whole wire. In the h110i SiNWs, the LUMO orbitals of HSiNWs and CH3SiNWs have similar distribution with more contribution from the central region while the OH groups make the LUMO concentrated in a rectangular region. For the HOMO orbitals of h110i SiNWs, CH3 groups make the orbitals more on the interior atoms compared to the HSiNWs. The OH groups change it dramatically and the HOMO moves to the top and bottom Si and O atoms. And more importantly, the HOMO and LUMO have considerable concentration on O atoms in both h100i and h110i SiNWs, indicating that the OHSiNWs are more sensitive to the external environment than the other two kinds of passivations. 99 Figure 4.8: The HOMO and LUMO orbital density of the h100i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respectively. Here red and blue represent HOMO and LUMO density, respectively. The contour is at 10% of the maximum value. 100 Figure 4.9: The HOMO and LUMO orbital density of the h110i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respectively. Here red and blue represent HOMO and LUMO density, respectively. The contour is at 10% of the maximum value. 101 Usually the surface substituents such as hydroxyl, halogens[126, 139] and methyl in h100i SiNWs from our simulations will induce big change in the spatial distribution of neargap orbitals. However, it is interesting to note that the CH3 substituents do not have signi cant e ect on the HOMO and LUMO orbitals in h110i SiNWs. In Figures 4.10 and 4.11 we depict the isosurface of the HOMO and LUMO wave functions in h110i HSiNWs, respectively. We can see that the wave functions have alternating positive and negative values along the wire axis and hence nodal surfaces perpendicular to the wire axis. When applying uniaxial tensile (compressive) strain, the energy values of the HOMO and LUMO will shift downward (upward) because the increase (decrease) of the distance between nodal planes results in the raising (lowering) of kinetic energy. So the conduction band minimum (CBM) and valence band maximum (VBM) will decrease or increase simultaneously. Similar e ect was predicted to be useful in separating electrons and holes when partial strain is applied in h110i SiNWs.[89, 140] Under 2% tension, LUMO is located in strained region while HOMO is in the regular part in a partially strained Silicon nanorod, suggesting photovoltaic application in terms of typeII homojunction solar cells. In our simulations, the CH3 passivation retains the HOMO and LUMO structures of alternating positive and negative values along the wire axis as shown in Figures 4.12 and 4.13. So the h110i CH3SiNWs are supposed to have similar capability of separat ing the positive and negative charge carriers under partial strain as in h110i HSiNWs. In addition to the charge carrier separation, the CH3SiNWs have better air stability and provide higher hole mobility and ono ratio than the HSiNWs [112] and do not a ect the conductance of the neargap channels, o ering quasiballistic transport within several subbands near the SiNWs band gap. [123] 102 Figure 4.10: The HOMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 103 Figure 4.11: The LUMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 104 Figure 4.12: The HOMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 105 Figure 4.13: The LUMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 106 4.8 E ective Mass of Charge Carriers The e ective mass of carriers was calculated from the band dispersions and sum marized in Table 4.3. The hole e ective mass mh of h100i HSiNWs studied is decreas ing with diameter in agreement with results of similar diameters reported by Yan, et al. [114] And in their calculations mh is much bigger than me for h100i HSiNWs. Although heavy hole is also observed in our simulations, the mh is greatly reduced with increased diameter and has comparable magnitude with me especially for the h100i wire of diameter d = 1.73 nm. The OH passivation has e ect on mh and me
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Firstprinciples Study of Electronic Properties of One Dimensional Nanostructures 
Date  20100501 
Author  Li, Junwen 
Keywords  DFT, Electronic structure, Graphene Nanoribbon, Helical symmetry, Quantum transport, Silicon nanowire 
Department  Physics 
Document Type  
Full Text Type  Open Access 
Abstract  In this work, we present simulation results for the electronic properties of two types of quasionedimensional nanostructures: graphene nanoribbons and silicon nanowires. The electronic structures throughout the dissertation were calculated within a firstprinciples, allelectron, selfconsistent local density functional (LDF) approach tailored for helical symmetry. Tightbinding models (TBM) were employed to study the energy gaps and bonding characteristics of armchair graphene nanoribbons. Within the study of silicon nanowires, we were interested in surface passivation effects and electronic transport properties of silicon nanowires with surface efects. The transport calculation makes use of the Landauer approach implemented with nonequilibrium Green's function formalism. For the study of graphene nanoribbons, we found that the zigzag ribbons exhibited spinpolarized ground states with ferromagnetic edges on opposite sides aligned antiferromagnetically. The electron states near the Fermi level in zigzag ribbons are localized on the edge carbon atoms. The armchair ribbons have ground states of degenerate spin configuration and the LDF energy gaps separate into three families, which is not consistent with the nearestneighbor TBM. The thirdnearestneighbor interaction was introduced to resolve the discrepancy between LDF and TBM in armchair ribbons. We found that armchair graphene ribbons have two different responses to twist deformation according to ribbon widths, suggesting potential nanomechanical application. In the investigation of silicon nanowires, we found that passivation using hydroxyl or methyl groups reduced the band gaps compared to hydrogenterminated Si nanowires, and passivation using methyl groups produced systems with indirect gaps for all <100> Si nanowires studied. All band gaps were direct in the <110> Si nanowires independent of passivation. The distributions of neargap orbitals were greatly affected by the surface substituents. We also found that the carrier effective masses of <100> Si nanowires were sensitive to the diameter and passivation, while those of <110> Si nanowires were not. We show that the hydroxyl defects can greatly reduce the conductance of hydrogenpassivated Si nanowires and can be used to tune the conductance of the silicon nanowires. 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  FIRSTPRINCIPLES STUDY OF ELECTRONIC PROPERTIES OF ONE DIMENSIONAL NANOSTRUCTURES By JUNWEN LI Bachelor of Science in Applied Physics Qingdao University Qingdao, Shandong, China 1999 Master of Science in Physics Institute of Physics, Chinese Academy of Sciences Beijing, China 2002 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial ful llment of the requirements for the Degree of DOCTOR OF PHILOSOPHY May, 2010 FIRSTPRINCIPLES STUDY OF ELECTRONIC PROPERTIES OF ONE DIMENSIONAL NANOSTRUCTURES Dissertation Approved: Dr. John W. Mintmire Dissertation Advisor Dr. Timothy M. Wilson Dr. Yin Guo Dr. Nicholas F. Materer Dr. A. Gordon Emslie Dean of the Graduate College ii ACKNOWLEDGMENTS I am pleased to acknowledge my appreciation for those who have helped this dissertation take shape. First and foremost, my sincere thanks go to my advisor, Dr. John Mintmire, for spending enormous time and e orts in guiding my research. Every week I obtain valuable feedback from our discussion at individual and group meetings. I will always remember the scenes when he derived formulas and elucidated my puzzles on paper or the board. This dissertation originates from these gradual accumulations. He provides strong support for my research projects with his insight and understanding of the research eld and grand knowledge. His rigorous scholarship is what I myself hope to achieve in my future research career. Dr. Mintmire also sets up an example for me to follow as an individual. I would like to thank my committee members: Dr. Timothy Wilson, Dr. Yin Guo, and Dr. Nicholas Materer for their time and patience in serving my committee and reviewing my research proposal and dissertation. I appreciate their valuable comments and advice on my work. I would like to thank Dr. Jacques Perk, Dr. Xincheng Xie, and Dr. Robert Hauen stein for their help and encouragement. I would also like to thank Dr. Carter White and Dr. Daniel Gunlycke at Naval Research Laboratory and Dr. Vincent Munier at Oak Ridge National Laboratory for our collaboration. My thanks must go to all the former and current members in the research group and fellow graduate students for their friendship and indispensable support, but I have to single out Dr. Shelly Elizondo, Dr. Thushari Jayasekera, Vasantha Jogireddy, iii and Pillalamarri Pavan for their help, friendship, and making a pleasant work envi ronment. Finally, I would like to thank my parents, my grandfather, my elder brother and his family for always being there for me, even though whole continents and oceans stand between us. Thanks to my wife, Jing Zhu, for her love, patience, and encouragement. She has always been an invaluable source of strength and support for me during the years we have been together. iv TABLE OF CONTENTS Chapter Page 1 INTRODUCTION 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . 4 2 THEORETICAL BACKGROUND AND APPROACHES 5 2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 HartreeFock Approximation . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 DFT Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.2 ThomasFermi Model . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.3 HohenbergKohn Theorems . . . . . . . . . . . . . . . . . . . 15 2.3.4 KohnSham Formulation . . . . . . . . . . . . . . . . . . . . . 17 2.4 Exchange Correlation Approximations . . . . . . . . . . . . . . . . . 21 2.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.2 Local Density Approximation . . . . . . . . . . . . . . . . . . 21 2.4.3 PBE Generalized Gradient Approximation . . . . . . . . . . . 23 2.5 Gaussian Basis Sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.6 Helical Band Structure Methods . . . . . . . . . . . . . . . . . . . . . 28 2.7 Landauer Approach for Quantum Conductance Calculation . . . . . . 29 3 GRAPHENE NANORIBBONS 34 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 v 3.2 Planar Zigzag Graphene Nanoribbons . . . . . . . . . . . . . . . . . . 36 3.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.2 Energetic Results . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.3 Nonspinpolarized LDF Band Structures . . . . . . . . . . . . 40 3.2.4 Ferromagnetic LDF Band Structures . . . . . . . . . . . . . . 45 3.2.5 Antiferromagnetic LDF Band Structures . . . . . . . . . . . . 49 3.2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3 Planar Armchair Graphene Nanoribbons . . . . . . . . . . . . . . . . 55 3.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.3.2 Model Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.3.3 Local Density Functional Results . . . . . . . . . . . . . . . . 58 3.3.4 Consideration of Longer Range Interactions . . . . . . . . . . 61 3.3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4 Twist E ect in Electronic Properties of Armchair Graphene Nanoribbons 66 3.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.4.2 Computational Approach . . . . . . . . . . . . . . . . . . . . . 68 3.4.3 Band Structures and Neargap Wavefunctions . . . . . . . . . 70 3.4.4 Structure Parameters . . . . . . . . . . . . . . . . . . . . . . . 73 3.4.5 Tightbinding Model . . . . . . . . . . . . . . . . . . . . . . . 75 3.4.6 Band Gap Change Under Twist . . . . . . . . . . . . . . . . . 77 3.5 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 81 4 SURFACE PASSIVATION EFFECTS IN SILICON NANOWIRES 82 4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Computational Approaches . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3 Structure Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.4 Mulliken Population Analysis . . . . . . . . . . . . . . . . . . . . . . 91 4.5 Band Gap Change . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 vi 4.6 Electronic Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.7 Neargap States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.8 E ective Mass of Charge Carriers . . . . . . . . . . . . . . . . . . . . 107 4.9 Surface Defects E ect on Transport . . . . . . . . . . . . . . . . . . . 109 4.9.1 Computational Approach . . . . . . . . . . . . . . . . . . . . . 109 4.9.2 Transport Conductance . . . . . . . . . . . . . . . . . . . . . . 111 4.10 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 113 5 CONCLUSIONS 114 vii LIST OF TABLES Table Page 3.1 Total energies divided by 2N and energy di erences between di erent spin con gurations of zigzag GNRs with various ribbon widths . . . . 39 3.2 Energy gaps of zigzag GNRs with various ribbon widths . . . . . . . 53 3.3 Band gaps Egap of armchair GNRs with ribbon width N . . . . . . . 60 3.4 Binding states of edge C pz orbitals HOMO and LUMO. " (#) indicates the energy point shifts upward (downward) with increasing twist angle. + () indicates the bonding (antibonding) states. . . . . . . . . . . 76 3.5 Band gaps (eV) as a function of twist angle (Deg) in armchair graphene nanoribbons with various widths N . . . . . . . . . . . . . . . . . . . 80 4.1 Band gaps (eV) as a function of diameter for h100i silicon nanowires with various passivations. X represents H, OH, or CH3. . . . . . . . . 95 4.2 Band gaps (eV) as a function of diameter for h110i silicon nanowires with various passivations. X represents H, OH, or CH3. . . . . . . . . 95 4.3 E ective mass of silicon nanowires along h100i and h110i directions with various passivations . . . . . . . . . . . . . . . . . . . . . . . . . 108 viii LIST OF FIGURES Figure Page 3.1 Sample zigzag GNRs with ribbon width N = 6. The numbering of the unit cells is shown at the left side of the ribbon with only the unit cell labeled 0 shown in their entirety. Each unit cell is composed of 2N C atoms and 2 H atoms. The 2N C atoms in unit cell 0 are numbered as shown. The index of zigzag chains is indicated at bottom from left to right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 (a) Total energy di erences divided by 2N between nonspinpolarized and ferromagnetic states in zigzag GNRs as a function of ribbon width N. (b) Total energy di erences divided by 2N between ferromagnetic and antiferromagnetic states in zigzag GNRs as a function of ribbon width N. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4 Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 42 3.5 Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 43 ix 3.6 Calculated HOMO band orbital densities of zigzag GNRs of width N = 6, calculated at wavevectors (a) , (b) 0.875 , (c) 0.75 , and (d) 0.625 , respectively. Calculated orbital densities for the LUMO band are given for wavevectors (e) , (f) 0.875 , (g) 0.75 , and (h) 0.625 , respectively. The isovalue is taken to be 0.008 for all plots. . . . . . . 44 3.7 Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.8 Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.9 Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.10 Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 51 3.12 Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. . . . . . . . . . . . . . . . . . . . . . . . . 52 3.13 Magnitude of the HOMOLUMO antiferromagnetic gap of the zigzag GNRs as a function of ribbon width N. . . . . . . . . . . . . . . . . . 53 x 3.14 (a) Symmetric armchair GNRs with ribbon width N = 7. Each unit cell is composed of 2N C atoms and 4 H atoms. (b) Staggered armchair GNRs with ribbon width N = 8. Each unit cell is composed of N C atoms and 2 H atoms. For (a) and (b) the numbering of the unit cells is shown at the left edge of each gure with only the unit cells labeled 0 shown in their entirety. The C atoms in unit cell 0 are numbered as shown. The index of dimer lines is indicated at bottom from left to right. 57 3.15 Calculated band structures of armchair GNRs of various widths. (a), (b), and (c) correspond to ribbon with N = 7, 8, and 9, respectively. (X) indicates the center (edge) of the rst Brillouin zone. . . . . . . 59 3.16 LDF energy gaps of the armchair GNRs as a function of width N. . . 60 3.17 Magnitude of the HOMOLUMO gap of the N armchair GNRs from the TBM calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.18 LDF results for Egap compared to the corresponding tightbinding model (TBM) results obtained from Eq. (3.3) with V1 = 3:2 eV and V3 = 0:3 eV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.19 (a) Armchair GNR with width N = 8. Each helical unit cell is com posed of one zigzag chain terminated by H atoms. The N C atoms in a helical cell are numbered as shown. The bonds are labeled as a1a11. (b) Axial view of 10 twisted armchair GNR with width N = 8. . . . 69 3.20 Electronic band structure of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . . . . . . . . . . . . . . . 71 3.21 The LUMO and HOMO wave function of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . The ribbons are vertically oriented. . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.22 CC bond lengths of the twisted armchair GNRs with width N = 8 as a function of twist angle . . . . . . . . . . . . . . . . . . . . . . . . 74 xi 3.23 Variation of energy gaps of armchair GNRs as a function of twist angle . (a), (b), and (c) represent the armchair GNRs of width N = 3m, 3m + 1, and 3m + 2, respectively. . . . . . . . . . . . . . . . . . . . . 79 4.1 Crosssection views of the studied HSiNWs. h100i SiNWs of diameter (a) 0.43 nm, (b) 0.87 nm, (c) 1.3 nm, and (d) 1.73 nm, respectively. h110i SiNWs of diameter (e) 0.73 nm, (f) 1.09 nm, and (g) 1.46 nm, respectively. The goldenyellow and silver balls represent Si and H atoms, respectively. In OHSiNWs and CH3SiNWs, all H atoms are replaced with OH and CH3 groups, respectively. . . . . . . . . . . . . 87 4.2 Band gaps as a function of diameter with hydrogen passivation for silicon nanowires along h100i. The square and blue line indicate the LDF results and the uptriangle and red line indicate the PBE results. 88 4.3 SiSi bond length as a function of the distance from bond center to wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm passivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4 Net charge as a function of distance from the atom to the wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm pas sivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. Si, H, O, and C are indicated by +, , , and , respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.5 Band gaps as a function of diameter with various surface passivations for silicon nanowires along (a) h100i, (b) h110i. . . . . . . . . . . . . 94 xii 4.6 Electronic band structures of h100i SiNWs (diameter d = 1.73 nm) passivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. . . . . . . . . . . . . 97 4.7 Electronic band structures of h110i SiNWs (diameter d = 1.46 nm) passivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. . . . . . . . . . . . . 98 4.8 The HOMO and LUMO orbital density of the h100i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respec tively. Here red and blue represent HOMO and LUMO density, respec tively. The contour is at 10% of the maximum value. . . . . . . . . . 100 4.9 The HOMO and LUMO orbital density of the h110i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respec tively. Here red and blue represent HOMO and LUMO density, respec tively. The contour is at 10% of the maximum value. . . . . . . . . . 101 4.10 The HOMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 103 4.11 The LUMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 104 4.12 The HOMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 105 4.13 The LUMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. . . . . . . 106 xiii 4.14 Model system for the conductance calculation. The system is composed of three regions: lefthand lead region, central conductor region, and righthand lead region. . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.15 The transport conductance through the OHpassivated region: The solid line is the conductance of the pure HSiNWs. The dashed line shows the conductance of the HSiNWs with OH defects on one helical cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 xiv CHAPTER 1 INTRODUCTION 1.1 Overview Nanotechnology deals with structures between 1 nanometer and 100 nanometer in size. The materials at this scale can display properties very di erent from their bulk counterparts. These unique nanoscale properties originate from size e ects such as the enhanced surfacetovolume ratio that causes surface atoms to experience po tentials di erent from those in the bulk. In recent years the rapid development of fabricating and manipulating materials at nanometer scale has opened up a great deal of opportunities in a variety of scienti c research e orts and practical applications. We can consider di erent numbers of spatial dimensions at nanometer scale. De pending on how many dimensions are con ned, the nanomaterials are termed as quantum well, quantum wire, and quantum dot for constraints on one, two, and three dimensions, respectively. Of these di erent systems, the quantum wire of one dimensional extension represents the simplest transport path and is of especially great signi cance in electronic applications. Two types of onedimensional nanostructures were investigated in this dissertation: graphene nanoribbons and silicon nanowires. Our study on graphene nanoribbons was motivated by the discovery of new exper imental methods for graphene, a single sheet of graphite. Novoselov, et al.,[1, 2] using a micromechanical cleavage method to extract a single sheet from threedimensional graphite, obtained the twodimensional form of carbon, i.e., graphene, which is stable under ambient conditions. The quality of the samples prepared were good enough so that the ballistic transport[1] and quantum Hall e ect[3] were observed easily. The 1 ballistic transport makes this new material a promising candidate for future electronic applications such as ballistic elde ect transistors (FETs). Although experimental work on graphene systems is recent, the theoretical study of graphene started long ago, with Wallace's tightbinding study in 1947 as an early work. [4] After that many properties of graphite and nanotubes were derived using graphene as starting point. The experimental achievement of graphene in 2004 attracted great attention immedi ately. Twodimensional graphene has a zero band gap with linear energy dispersion near the Fermi level and crossing at the Dirac point between valence and conduc tion bands. For practical electronic application in semiconductor industry a band gap could be induced by fabricating graphene nanoribbons. The ribbon structures have already been successfully obtained with various techniques such as lithographic patterning,[5] chemical vapor deposition,[6] Joule heating,[7] and unzipping the car bon nanotubes (CNTs).[8, 9] The other onedimensional structures investigated, silicon nanowires, have a wide variety of applications from use in batteries to use in energy conversion. Because of the geometry of silicon nanowire arrays, they are capable of in ating 4 times their normal size when absorbing lithium ions, which enables the new kind of battery to hold 10 times the charge of existing lithiumion batteries.[10] The onedimensional structure of silicon nanowires can allow pin coaxial solar cells consisting of a ptype silicon nanowire core surrounded by i and ntype silicon shells.[11] An advantage of this core/shell architecture is that carrier separation takes place in the radial rather than the axial direction, with a carrier collection distance smaller or comparable to the mi nority carrier di usion length. Si nanowire arrays show promise as highperformance, scalable thermoelectric materials because of the high electron conductivity and low phonon conductivity in silicon nanowires with rough surface.[12, 13] The promising applications of these nanostructures require a better theoretical understanding of their electronic properties. Our computational simulations of elec 2 tronic properties of graphene nanoribbons and silicon nanowires are of great signi  cance from both fundamental and practical points of view. 3 1.2 Organization of the Dissertation Chapter 2 presents a brief introduction of the theoretical background and com putational approaches utilized throughout this work. We brie y review the theo retical formalism of density functional theory, helical band structure methods, and the transport conductance approach with Landauer formula implemented by using nonequilibrium Green's function method. Chapter 3 is dedicated to the simulation results from a study carried out on graphene nanoribbons. This theoretical study was primarily motivated by the exper imental implementation of isolated graphene sheets which are a single layers of carbon atoms connected with hexagonal patterns. In this work, we investigated two kinds of graphene nanoribbons, namely zigzag and armchair graphene nanoribbons according to the edge shape. For zigzag nanoribbons, the calculations on energetic analysis and corresponding band structures were carried out. For armchair nanoribbons, we studied the electronic structure dependence on the ribbon width and on the torsional deformation. We have four papers related to the ribbon studies. [14{17] Chapter 4 is devoted to the investigation carried out involving surface passivation e ect in silicon nanowires. Because of the enhanced surfacetovolume ratio in silicon nanowires compared to bulk silicon, the electronic properties of the nanowires are strongly dependent on the surface substituents. We studied three kinds of passivating functional groups such as hydrogen, hydroxyl, and methyl. We also investigated the hydroxyl surface defects e ect on transport property of pure hydrogenpassivated silicon nanowires. We have two papers based on the simulations presented in this chapter.[18, 19] Chapter 5 summarizes this dissertation, focusing on the simulation ndings and the possible further investigation based on the present calculations. 4 CHAPTER 2 THEORETICAL BACKGROUND AND APPROACHES 2.1 Overview All band structure calculations reported in this work were carried out by us ing the Helical Nanostructures (HENS) parallax code developed by Dr. Mintmire at Oklahoma State University. Di erent from regular band structure codes that use translational symmetry, this approach takes advantage of the helical symme try in onedimensional nanostructures such as polymers, nanoribbons, nanotubes, and nanowires. Taking account of the helical symmetry can greatly reduce the unit cell size, and therefore can achieve major computational savings. Fujita and Imamura[20, 21] rst proposed this idea of utilizing helical symmetry when studying the electronic structures of polymers with tight binding model. And this approach was developed for rstprinciples methods by Dr. Mintmire while at the Naval Research Laboratory. Mintmire, et al.,[22] using this helical band structure code, successfully predicted that the armchair singlewalled carbon nanotubes would be metallic four years prior to the experimental veri cation. Now we apply this code in the study of graphene nanoribbons and silicon nanowires. The zerobias conductance was calculated within a Landauer approach.[23, 24] The transport code was developed by using the nonequilibrium Green's function (NEGF) formalism within a linear response regime. The approach is capable of taking any general Hamiltonian that can be described with a localized orbital basis, such as orthogonal, nonorthogonal tightbinding models, and rstprinciples results by making use of Gaussian orbitals as in the study of surface defects e ect in silicon 5 nanowires. The e ective Hamiltonian matrix elements and orbital overlap matrix elements used in the investigation of silicon nanowires were generated by the HENS parallax selfconsistent eld (SCF) calculations. In this chapter, we will brie y describe the background on which the codes are based, including the density functional theory (DFT), the helical band structure method and the NEGF approach. Before that the HartreeFock method is reviewed in order to gain some idea on the exchange interaction and the advantage of DFT when dealing with manyelectron system. The discussion in this chapter is mainly adapted from Refs. [25{27]. 6 2.2 HartreeFock Approximation For an isolated Nelectron atomic or molecular system, in the BornOppenheimer approximation the timeindependent nonrelativistic Schr odinger equation is given by ^H = E ; (2.1) where E is the electronic energy, = (x1; x2; :::; xN) is the Nelectron wavefunction, and ^H is the Hamiltonian operator composed of three terms, ^H = ^ T + ^ Vext + ^ Vee; (2.2) where ^ T = XN i=1 1 2 r2i (2.3) is the kinetic energy operator, ^ Vext = XN i=1 XM =1 Zi jr R j (2.4) is the electronnucleus attraction energy operator, and ^ Vee = XN i6=j 1 jri rj j (2.5) is the electronelectron repulsion energy operator. When additional elds are present, ^ T and ^ Vee remain unchanged and extra terms appear in ^ Vext. The total energy ETOT of the system is composed of the electronic energy E and the nucleusnucleus repulsion energy Vnn = XM < Z Z R : (2.6) So we may have ETOT = E + Vnn: (2.7) 7 The coordinate xi of electron i comprises spatial coordinates ri and spin freedom si. Atomic units are employed here and throughout this dissertation (unless otherwise speci ed): the length unit is the Bohr radius, a0 = 0:5292 A, the charge unit is the charge of the electron, e, and the mass unit is the mass of the electron, me. The ^ T and ^ Vext in Eq. (2.2) are singleparticle operators and the ^ Vee is a two particle operator from which the complications of manybody problem arises. We need to make approximation on the manybody wavefunction (x1; x2; :::; xN). The simplest way to approximate the wavefunction is through the Hartree approximation, where the manybody wavefunction is replaced by a product of singleparticle orbitals, i(xi), (x1; x2; :::; xN) = 1 p N 1(x1) 2(x2) N(xN); (2.8) where i(xi) is a combination of spatial function i(ri), and spin function (si), such that i(xi) = i(ri) (si); (2.9) where = ; indicate spinup and spindown, respectively. However, the electrons are Fermions and should obey the Pauli exclusion principle which states that one quantum state only can be occupied by at most one electron. The wavefunction should be antisymmetric with respect to the interchange of any two electron coordi nates xi and xj . Eq. (2.8), however, does not satisfy (x1; :::; xi; :::; xj ; :::; xN) = (x1; :::; xj ; :::; xi; :::; xN): (2.10) To account for Pauli principle the HartreeFock approximation was proposed by writing the wavefunction as an antisymmetrized product of orbitals explicitly. The HartreeFock wavefunction HF amounts to a linear combination of the terms in Eq. (2.8) including all permutations of the electron coordinates with the corresponding 8 weights 1, which can be expressed as a Slater determinant, HF = 1 p N! 1(x1) 1(x2) 1(xN) 2(x1) 2(x2) 2(xN) ... ... ... ... N(x1) N(x2) N(xN) : (2.11) The HartreeFock energy can be evaluated by taking the expectation value of the Hamiltonian with respect to the Slater determinant Eq. (2.11). This yields EHF = h HF j ^H j HF i = XN i=1 Z i (r) 1 2 r2 + vext(r) i(r) dr + 1 2 XN i=1 XN j=1 Z Z j i(r1)j2j j(r2)j2 jr1 r2j dr1dr2 1 2 XN i=1 XN j=1 Z Z i (r1) i(r2) j (r2) j(r1) jr1 r2j si;sj dr1dr2: (2.12) The last term is of signi cant interest since it arises from the antisymmetric nature of the HartreeFock wavefunction. It vanishes when si 6= sj , which is a result of the Pauli exclusion principle. Therefore it is termed as exchange energy. In order to obtain the total energy of the system an extra term coming from repulsion energy between the nucleus must be added to Eq. (2.12). In most cases, we are more concerned with the ground state of a system. Even from the ground state property of system, we could get a great deal of information on the property of excited states in room temperature since the system deviates not too much from the ground states. To obtain the HartreeFock ground state energy EHF 0 we could minimize Eq. (2.12) with respect to the orbitals, subject to the constraint that the orbitals remain orthonormal (h ij ji = ij). The minimization procedure carried out with the Euler 9 Lagrange method yields the corresponding stationary condition given by EHF 0 XN i=1 XN j=1 "i (h ij ji 1) ! = 0; (2.13) where "i is the Lagrange multiplier. The corresponding Euler equations are written as 1 2 r2 + vext(r) + XN j=1 Z j j(r0)j2 jr r0j dr0 ! i(r) XN j=1 Z i(r0) j (r0) j(r) jr r0j si;sj dr0 = "i i(r); (2.14) which is called HartreeFock equations. In general the HartreeFock equations can not be solved analytically. One excep tion is for the homogeneous electron gas, where the solutions are plane wave functions because of the constant external potential everywhere. Otherwise, the HartreeFock equations are solved using an iterative method. Because the desired orbitals are re quired to construct the oneelectron e ective potential, this process is known as the selfconsistent eld procedure. The selfconsistent procedure starts with an initial guess for the orbitals, and successive iterations are performed with new orbitals gen erating the new potentials until the convergence is achieved. The converged orbitals are the ground state orbitals for that system within the HartreeFock approximation. HartreeFock theory is not an exact theory because it only considers a single determinant for the electron wavefunction. The only case where a single determinant is the exact solution is for a noninteracting system of electrons. In real systems the motions of electrons are more correlated than the mean eld description provided by HartreeFock approximation. The interaction energy missed in HartreeFock approximation is termed as the correlation energy EC, EC = E0 EHF ; (2.15) 10 where E0 is the exact ground state energy. Since HartreeFock is a variational ap proach, EHF E0 always holds and the correlation energy takes nonpositive value. A natural way to include correlation e ects beyond the HartreeFock approxima tion is to represent the manyelectron wavefunction as a linear combination of Slater determinants corresponding to ground state and excited states. These post Hartree Fock methods, such as con guration interaction, coupledcluster and M llerPlesset theory have been extensively developed in quantum chemistry. Although the pre cision may be systematically improved by including more and more excited Slater determinants, the computational cost increases dramatically with the number of ex citation levels. Consequently, these post HartreeFock approaches are limited to small systems such as atoms and small molecules. 11 2.3 DFT Methods 2.3.1 Overview The main idea of density functional theory (DFT) is to describe a system composed of N interacting electrons via the electron density rather than the manyelectron wavefunction like the Slater determinant used in HartreeFock approximation. This means that the basic variable of the system depends only on 3 spatial coordinates (x, y, and z) rather than 3N degrees of freedom as in HartreeFock approximation. Additionally, the computational costs are relatively low compared to those methods based on the complicated manyelectron wavefunction such as HartreeFock theory and its descendants. In 1927 Thomas and Fermi rst proposed a model expressing the electronic energy as a functional of electron density.[28, 29] In the original idea they derived a di erential equation for the electron density without resorting to oneelectron orbitals. The DFT has its conceptual roots in the ThomasFermi model, but the rm theoretical foundation was set up by two HohenbergKohn theorems 37 years after the proposal of ThomasFermi method. 2.3.2 ThomasFermi Model Thomas and Fermi proposed that the electronic energy of a system can be ex pressed as the functional of electron density. In this idea the kinetic, exchange, and correlation contributions can be constructed from the model study on the homoge neous electron gas but should be dependent on the position. The electron density (r) is the central variable, given by (r) = N Z Z dr2 drN (r; r2; :::; rN) (r; r2; :::; rN): (2.16) 12 The total energy of a system is expressed as a functional ETF [ (r)], which is given by ETF [ (r)] = CF Z (r)5=3 dr + Z (r)vext(r) dr + 1 2 Z Z (r1) (r2) jr1 r2j dr1dr2: (2.17) The rst term in Eq. (2.17) is the electronic kinetic energy associated with non interacting homogeneous electron gas. This form could be obtained by integrating the kinetic energy density of a homogeneous electron gas t0[ (r)], TTF [ (r)] = Z t0[ (r)] dr; (2.18) where t0[ (r)] is obtained by summing all of the free electron energy states "k = k2 2 , up to the Fermi wave vector kF = [3 2 (r)]1=3, t0[ (r)] = 2 8 3 Z k2 2 nk dk = 1 2 Z kF 0 k4 dk (2.19) with nk is the density of states in reciprocal space. This leads to the form given in Eq. (2.17) with coe cient CF = 3 10(3 2)2=3. The second term is the classical electrostatic energy of attraction between the nucleus and the electrons, where vext(r) is the static Coulomb potential arising from the nucleus vext(r) = XM =1 Z jr R j : (2.20) Finally, the third term in Eq. (2.17) represents the electronelectron interactions of the system. In the ThomasFermi model this term contains only the classical Coulomb repulsion energy between electrons, known as the Hartree energy. In order to obtain the ground state density and energy of a system, we may min imize the ThomasFermi energy functional of Eq. (2.17) with constraint of conserved total number of electrons N. Applying the EulerLagrange method to Eq. (2.17) leads to the stationary condition ETF [ (r)] Z (r) dr N = 0; (2.21) which yields the socalled ThomasFermi equations 5 3 CF (r)2=3 + vext(r) + Z (r0) jr r0j dr0 = 0; (2.22) 13 which can be solved directly to obtain the ground state density. ThomasFermi theory su ers from many di culties. One of them is that it does not predict bonding between atoms, so molecules and solids can not exist in this theory. The main source of error is a crude approximation for the kinetic energy. The kinetic energy represents a substantial portion of the total energy of a system and so even small errors prove disastrous. Another shortcoming is the oversimpli ed description of the electronelectron interactions, which are treated classically and so do not take into account quantum mechanical e ects such as the exchange interaction. Shortly after the introduction of ThomasFermi theory, Dirac[30] developed an approximation for the exchange interaction based on the homogeneous electron gas. The resulting formula is simple, and is also a local functional of the density, EX[ (r)] = CX Z (r)4=3 dr (2.23) with CX = 3 4 3 1=3 . When exchange interaction is treated like Eq. (2.23), the theory is called ThomasFermiDirac Model. Correlation can also be easily included by using any local approximation derived from homogeneous electron gas. One commonly used was proposed by Wigner[31], EC[ (r)] = 0:056 Z (r)4=3 0:079 + (r)1=3 dr: (2.24) The ThomasFermi model was actually too crude, mainly because the approxima tion used for the kinetic energy of the electrons was unable to sustain bound states. However, it set up the basis for the later development of density functional theory. 14 2.3.3 HohenbergKohn Theorems The ThomasFermi approach was developed in the hopes that the energy can be expressed exclusively in terms of the electron density. This idea, however, was intuitive at the time. Until 1964, Hohenberg and Kohn[32] proved two theorems that put the ThomasFermi model on solid mathematical grounds. Theorem 1. The external potential vext(r) is determined, within a trivial additive constant, by the electron density (r). Let us suppose there are two di erent external potentials, vext;1(r) and vext;2(r) such that they correspond to same ground state electron density (r). Let 1 and E1 = h 1j ^H 1j 1i be the ground state wavefunction and ground state energy of the Hamiltonian ^H 1 = ^ T + ^ Vext;1 + ^ Vee. Let 2 and E2 = h 2j ^H 2j 2i be the ground state wave function and ground state energy of the Hamiltonian ^H 2 = ^ T + ^ Vext;2 + ^ Vee. We assumed that di erent Hamiltonians correspond to di erent ground state wave functions, i.e., 1 6= 2. Applying the variational principle we may have: E1 < h 2j ^H 1j 2i = h 2j ^H 2j 2i + h 2j ^H 1 ^H 2j 2i (2.25) = E2 + Z (r)[vext;1(r) vext;2(r)] dr (2.26) Similarly, we have E2 < h 1j ^H2j 1i = h 1j ^H 1j 1i h 1j ^H 1 ^H 2j 1i (2.27) = E1 Z (r)[vext;1(r) vext;2(r)] dr (2.28) Therefore adding the two inequalities leads to the result, E1 + E2 < E2 + E1; (2.29) which is a contradiction, and as a result there can not be two di erent external potentials that correspond to the same electron density for the ground state, unless they di er by a trivial additive constant. 15 Theorem 2. The ground state energy can be obtained variationally : the density that minimizes the total energy is the exact ground state density. As just shown, (r) determines vext(r), N and vext(r) determine ^H and therefore . This ultimately means is a functional of (r), and so the expectation value of any operator ^O is also a functional of (r), i.e., O[ (r)] = h [ (r)]j^O j [ (r)]i: (2.30) Assume we have 0(r), which determines its own v0 ext(r) and wavefunction 0. Then 0 could be taken as a trial function for the system of interest having external potential vext and corresponding (r). Thus, h 0j ^H j 0i = h 0j ^ Tj 0i + h 0j ^ Veej 0i + h 0j ^ Vextj 0i = E[ 0(r)] E[ (r)] (2.31) based on the variational argument. 16 2.3.4 KohnSham Formulation Although the HohenbergKohn theorems are extremely powerful, they do not o er a way to calculate the ground state density of a system in practice. About one year after the seminal DFT paper by Hohenberg and Kohn, Kohn and Sham (KS)[33] devised a simple method for carrying out DFT calculations, that retains the exact nature of DFT. Since the electron density is the central variable rather than the wavefunction, we could construct our electron density from orbitals which are obtained from easily solved system such as the noninteracting system. The KohnSham formulation centers on mapping the full interacting system with the real potential, onto a ctitious noninteracting system where the electrons are moving in an e ective singleparticle KS potential vKS(r). The KohnSham method is still exact since it yields the same ground state electron density as the real system, but greatly facilitates the calculation. First consider the variational problem presented in the second HohenbergKohn theorem. The ground state energy of a many electron system can be obtained by min imizing the energy functional, subject to the constraint that the number of electrons N is conserved, which leads to F[ (r)] + Z vext(r) (r) dr Z (r) dr N = 0; (2.32) where F[ (r)] = T[ (r)] + Vee[ (r)] is called universal functional since they do not depend on the external potentials. The Euler equation is given by = F[ (r)] (r) + vext(r); (2.33) where is the Lagrange multiplier associated with the constraint of conserved N. The idea of Kohn and Sham was to set up a system where the kinetic energy could be determined exactly, since this was a major problem in ThomasFermi theory. This 17 was achieved by resorting to a noninteracting system of electrons. The correspond ing ground state wavefunction KS for this type of system is given exactly by a determinant of singleparticle orbitals i(ri), KS = 1 p N! det[ 1(r1) 2(r2)::: N(rN)]: (2.34) The universal functional F[ (r)] was then partitioned into three terms, the rst two of which are known exactly and constitute the majority of the energy, the third being a small unknown quantity, F[ (r)] = TS[ (r)] + EH[ (r)] + EXC[ (r)]: (2.35) TS[ (r)] is the kinetic energy of a noninteracting electron gas of density (r), EH[ (r)] is the Hartree energy of the electrons EH[ (r)] = 1 2 Z Z (r1) (r2) jr1 r2j dr1dr2; (2.36) and EXC[ (r)] is the exchange correlation energy, which contains the di erence be tween the exact and noninteracting kinetic energies and also the nonclassical contri bution to the electronelectron interactions, of which the exchange energy is a part. In the KohnSham prescription the Euler equation now becomes = TS[ (r)] (r) + vKS(r) (2.37) where the e ective KS potential vKS(r) is given by vKS(r) = vext(r) + vH(r) + vXC(r) (2.38) with the Hartree potential vH(r), vH(r) = EH[ (r)] (r) = Z (r0) jr r0j dr0; (2.39) and the exchange correlation potential vXC(r), vXC(r) = EXC[ (r)] (r) : (2.40) 18 In KohnSham theory the several potential terms were just rearranged to make up vKS. So the density obtained when solving the ctitious noninteracting KohnSham system is the same as the exact ground state density. The ground state density is obtained in turn by solving the N oneelectron Schr odinger equations, 1 2 r2 + vKS(r) i(r) = "i i(r); (2.41) where "i are Lagrange multiplier corresponding to the orthonormality of the N single particle states i(r), and the density is constructed from (r) = XN i=1 j i(r)j2: (2.42) The noninteracting kinetic energy TS[ (r)] is therefore given by TS[ (r)] = 1 2 XN i=1 Z i (r)r2 i(r) dr: (2.43) Since vKS(r) depends on the density through the exchange correlation and Hartree potentials, the KohnSham equations must be solved selfconsistently as in the Hartree Fock method. In order to handle the kinetic energy in an exact manner, N equations have to be solved in KohnSham theory to obtain the set of Lagrange multiplier "i, as opposed to one equation that determines when solving for the density directly, as in the ThomasFermi approach. However an advantage of the KohnSham method is that as the complexity of a system increases, with N increasing, the problem becomes no more di cult, only the number of singleparticle equations to be solved increases. Although exact in principle, KohnSham theory is approximate in practice because of the unknown exchange correlation functional EXC[ (r)]. An implicit de nition of EXC[ (r)] can be given by, EXC[ (r)] = T[ (r)] TS[ (r)] + Eee[ (r)] EH[ (r)] (2.44) where T[ (r)] and Eee[ (r)] are the exact kinetic and electronelectron interaction energies respectively. The intention of Kohn and Sham was to make the unknown 19 contribution to the total energy of the noninteracting system as small as possible, and this is indeed the case with the exchange correlation energy, however it is still an important contribution since the binding energy of many systems is about the same size as EXC[ (r)], so an accurate description of exchange and correlation is crucial for the prediction of binding properties. 20 2.4 Exchange Correlation Approximations 2.4.1 Overview While DFT is an exact theory of ground state properties, practical applications of DFT must be based on approximations for the unknown exchange correlation po tential which describes the e ects of Pauli exclusion principle and Coulomb potential beyond the pure electrostatic interaction between electrons in an average sense. If the exact exchange correlation potential is obtained, the manybody problem can be solved exactly. All exchange correlation functionals can be written in a general form EXC[ (r)] = Z (r)"XC(r) dr; (2.45) where "XC(r) is the exchange correlation energy density. In DFT, the exchange correlation functionals can be approximated in di erent levels by the number and kind of their local ingredients.[34] The simplest one is the local density approximation in which only the local electron density is considered. In the high level approximation such as generalized gradient approximation (GGA), the gradient of electron density r (r) is also included. We will discuss local density approximation (LDA) and GGA especially PBE with more details because these two kinds of functionals are extensively used in our simulations. 2.4.2 Local Density Approximation In 1951, Slater[35] suggested a model in which the kinetic energy would be treated as in the HartreeFock model, but where the exchange term was replaced by a func tional of electron density vX[ (r)]. This work lead to what is known as the X SCF method. In the derivation of Slater's oneelectron exchange potential, we follow the 21 procedure of Kohn and Sham.[33] We may begin by writing the HartreeFock exchange operator in the form of an equivalent potential acting on the kth wave function, vxk(x) = XN k0=1 Z k(x) k0(x0) k0(x) k(x0) jr r0j dx0 k(x) k(x): (2.46) Next we make an approximation and simpli cation assuming the wave functions can be approximated by plane waves as in free electron gas. This leads to vxk(r) = kF (r) 1 + k2F (r) k2 2kkF (r) ln k + kF (r) k kF (r) (2.47) with kF (r) = (3 2 (r))1=3. Then, we average vxk over the occupied state k, which results in vX[ (r)] = 3 2 (3 2)1=3 (r)1=3; (2.48) which is the original Slater's exchange potential. Since electron density adjustments are mainly a ected by redistribution of the electrons near the Fermi level, it is reasonable to take k = kF (r) in Eq. (2.47), which is equivalent to taking the e ective exchange potential for a state at the top of the Fermi distributions. This leads to vX[ (r)] = 1 (3 2)1=3 (r)1=3; (2.49) which is the expression for exchange potential that Kohn and Sham derived in 1965. Eq. (2.49) is di erent from Slater's by a factor of 2 3 . The paper published by Kohn and Sham in 1965 is not the rst addressing this di erence of twothirds. Actually in 1954, G asp ar already obtained the same kind of dependence of the exchange energy on the electron density as that of Kohn and Sham. The model, in which vX[ (r)] is proportional to (r)1=3, became known as the X method, where the exchange term is written in the form, 3 2 (3 2)1=3 (r)1=3 (2.50) 22 with = 1 for Slater's model and = 2=3 for that of G asp ar, Kohn, and Sham. In Eq. (2.50), we only include the exchange e ect. Together with this exchange potential, a commonly used correlation formula is that of Perdew and Zunger[36] which makes use of accurate quantum Monte Carlo data for the homogeneous electron gas generated by Ceperley and Alder[37] to x the coe cients in the interpolation formula. Our local density functional (LDF) approach employs the G asp arKohnSham exchange potential and neglects the correlation e ect. This approach is also commonly referred to as the local density approximation (LDA). 2.4.3 PBE Generalized Gradient Approximation For exchange correlation energy the LDA makes use of the result obtained from homogeneous electron gas at each point irrespective of the nonhomogeneity of the real electron density and is the simplest approximation for exchange and correla tion functionals. For real system of nonhomogeneous electron density the exchange correlation energy can be signi cantly di erent from that of homogeneous electron gas. The gradient and higher spatial derivatives of the electron density are needed to account for this deviation. Perdew and colleges made great contributions on the development of GGA func tionals. [34, 36, 38{42] They introduced an analytic function known as the enhance ment factor, FXC[ (r);r (r)], that accounts for the nonhomogeneity of electron density and modi es the LDA energy density through EGGA XC [ (r)] = Z (r)"hom XC [ (r)]FXC[ (r);r (r)] dr: (2.51) Usually the GGA enhancement factor is written in terms of the Seitz radius, rs, which is related to the electron density (r) by rs = 3 4 (r) 1=3 ; (2.52) 23 and the dimensionless reduced density gradient s(r), s(r) = jr (r)j 2kF (r) (r) (2.53) with kF being the Fermi wave vector, kF (r) = 3 2 (r) 1=3 : (2.54) At present the most popular GGA functional in physics community is PBE.[38] The PBE functional EPBE XC is composed of two parts : PBE exchange EPBE X and PBE correlation EPBE C functionals. The exchange PBE functional is written in the following form, EPBE X [ (r)] = Z d3r (r)"hom X ( )FX(s) (2.55) with FX(s) = 1 + k k 1 + s2=k (2.56) where k = 0:804; = 0:21951. This exchange energy obeys the spinscaling relation ship,[43] EPBE X [ ; ] = 1 2 EPBE X [2 ] + 1 2 EPBE X [2 ]: (2.57) And the PBE correlation energy functional is given by EPBE C [ ; ] = Z d3r (r)["hom C (rs; ) + H(rs; ; t)]; (2.58) where the nonlocal part H(rs; ; t) depends on the parameter t including the density gradient via H = 3 ln 1 + t2 1 + At2 1 + At2 + A2t4 (2.59) where = 0:066725; = 0:03191, and = ( )= (r) t = ( =4)1=2(9 =4)1=6 2 ( )r1=4 s ( ) = 1 2 (1 + )2=3 + (1 )2=3 A = ( = ) e"C(rs; )= 3 1 : 24 The GGA takes account of the gradient of the electron density. For systems of slowly varying electron density, the GGA has proved to be an improvement over LDA. For example, GGAs lead to improvement on total energies,[44] atomization energies,[44, 45] energy barriers, and structural energy di erences.[46, 47] 25 2.5 Gaussian Basis Sets When molecular calculations are performed, basis functions are employed to con struct the molecular orbitals, which are expanded as a linear combination of such functions with the coe cients to be determined in the self consistent procedure. A natural choice for basis functions would be Slatertype atomic orbitals decaying exponentially with distance from the nuclei. Boys[48] pointed out that these Slater type orbitals could be approximated as linear combinations of Gaussiantype orbitals (GTOs) instead. Because it is easier to calculate overlap and other integrals involved in the construction of Hamiltonian matrix with Gaussian basis functions, this led to signi cant computational savings. In Cartesian coordinate system, Gaussian primitive basis functions are de ned as products of the integer powers of Cartesian coordinates and a Gaussian function: (r; ; ^n) = xnxynyznze r2 (2.60) and the corresponding normalization constant is given by N( ; ^n) = 2 3=4 (4 )n=2 p (2nx 1)!!(2ny 1)!!(2nz 1)!! (2.61) with n = nx + ny + nz. At present, there are various basis sets composed of Gaussiantype orbitals. The smallest one are called minimal basis sets typically composed of a minimum number of basis functions required to represent all of the electrons on each atom. The most common minimal basis set is STOnG, where n is an integer. The value of n represents the number of Gaussian primitive functions comprising a single basis function. In these basis sets, the same number of Gaussian primitives are used to construct core and valence orbitals. Minimal basis sets are usually used in testing calculations to get the rough idea. 26 Other more complicated basis sets are proposed to improve the quality of the minimal basis set. Based on the observation that the valence orbitals have greater e ects on chemical properties, a splitvalence basis set was proposed such that the number of functions used to describe the valence electrons is doubled but a single function is used for the inner shells. The notation for the splitvalence basis sets is typically XYZG. In this convention, X represents the number of primitive Gaussians comprising each core atomic orbital basis function. The Y and Z indicate that the valence orbitals are composed of two basis functions each, the rst one composed of a linear combination of Y primitive Gaussian functions, the other composed of a linear combination of Z primitive Gaussian functions. For example, in the 321G basis set, three Gaussian functions describe the core orbitals and the valence electrons are also represented by three Gaussians of which two are used for the contracted part and one for the di use part. Other improvement is the addition of polarization functions, indicated by an as terisk, . Two asterisks, , indicate that polarization functions are also added to light atoms (hydrogen and helium). These polarization functions have one more node. For example, the only basis function located on a hydrogen atom in a minimal basis set would be a function approximating the 1s atomic orbital. When polarization is added to this basis set, a pfunction is also added to the basis set. This adds some extra exibility within the basis set, e ectively allowing molecular orbitals involving the hydrogen atoms to be more asymmetric about the hydrogen nucleus. Throughout this work, we use 631G* split valence basis sets to model the twisted armchair graphene nanoribbons in Chapter 3. The 321G split valence set is used in the simulations of silicon nanowires discussed in Chapter 4. And we also use uncontracted basis sets[49] such as 7s3p for C atoms and 3s for H atoms in the study of planar graphene nanoribbons presented in Chapter 3. 27 2.6 Helical Band Structure Methods For a onedimensional system with helical symmetry such as carbon nanotubes, we can de ne a screw operation in terms of a translation l down the z axis in conjunction with a righthanded rotation about the z axis. That is, S r = S(a; ) r = 0 BBBB@ x cos y sin x sin + y cos z + l 1 CCCCA : (2.62) Because the symmetry group generated by the screw operation S is isomorphic with the onedimensional translation group, we may have a generalized Bloch's theorem in which the oneelectron wavefunctions will transform according to Sm i(r; ) = ei m i(r; ): (2.63) The quantity is a dimensionless quantity which is conventionally restricted to < . For the case of pure translation, i.e., = 0, corresponds to a normalized quasimomentum in terms of = kl, where k is the traditional reciprocal wavevector. The oneelectron wavefunctions i are constructed from a linear combination of Bloch functions 'j , which are in turn constructed from a linear combination of nuclearcentered Gaussiantype orbitals j(r) in forms of i(r; ) = X j cji( ) 'j(r; ) (2.64) 'j(r; ) = X m ei m Sm j(r): (2.65) The detailed description and formulism derivation of this helical method can be found in Ref. [50]. 28 2.7 Landauer Approach for Quantum Conductance Calculation The conductance calculations have been extensively discussed in literature. This section is mainly adapted from Refs. [51, 52] with more details carried out. The transport calculations are based on phasecoherent transport of electrons through the device region from the left semiin nite lead to the right semiin nite lead of a linear system. In Landauer approach,[23, 24] the conductance is proportional to the transmission function as follows Gcon = 2e2 h T; (2.66) where both the conductance Gcon and the transmission function T are functions of energy E. T represents the probability that an electron injected at one end of the conductor will emit at the other end. The transmission function T of the device can be expressed in terms of the Green's functions of the conductor and the coupling of the conductor to the leads, T = Tr[LGRRGA]; (2.67) where the advanced Green's function GA is the Hermitian conjugate of the retarded Green's function GR of the conductor describing the dynamics of the electrons inside the conductor, and the fL;Rg accounts for the coupling of the conductor to the leads. To compute the Green's function of the conductor we could start with the retarded Green's function of the whole system including the conductor and two semiin nite leads, (" H)G = I; (2.68) where " = E +i ( is an arbitrarily small positive real number) is a complex energy, I is the identity matrix, and we neglect the superscript in G for simplicity. Since the whole system can be conceptually divided into three distinct regions : a conductor 29 region, a lefthand lead, and a righthand lead, Eq. (2.68) can be expressed in terms of submatrices that correspond to di erent subsystems, 0 BBBB@ (" HL) hLC 0 hy LC (" HC) hCR 0 hy CR (" HR) 1 CCCCA 0 BBBB@ GL GLC GLCR GCL GC GCR GLRC GRC GR 1 CCCCA = I: (2.69) Since we are more concerned with the Green's function of conductor, the equations involving the GC can be written explicitly as follows : (" HL)GLC + hLCGC + 0 = 0 (2.70) hy LCGLC + (" HC)GC + hCRGRC = I (2.71) 0 + hy CRGC + (" HR)GRC = 0 (2.72) From Eq. (2.72) we have GRC = (" HR)1hy CRGC: (2.73) From Eq. (2.71) we get GLC = (" HL)1hLCGC: (2.74) Substituting the expressions for GRC and GLC into Eq. (2.72) leads to hy LC(" HL)1hLCGC + (" HC)GC hCR(" HR)1hy CRGC = I; (2.75) which suggests GC = h (" HC) hy LC(" HL)1hLC hCR(" HR)1hy CR i1 : (2.76) Then we can write the expression of the retarded Green's function of a system as GR = [(" HC) RL RR ]1 (2.77) where Rf L;Rg is the retarded selfenergy terms describing the coupling between the conductor and semiin nite leads, and is given by RL = hy LCgR LhLC; RR = hCRgR Rhy CR (2.78) 30 with gR L and gR R are the Green's functions of the semiin nite leads ("HL)1 and (" HR)1, respectively. The selfenergy term can be viewed as an e ective Hamiltonian term coming from the interaction of the conductor with leads. The coupling matrices fL;Rg can be easily obtained as fL;Rg = i[ Rf L;Rg Af L;Rg]; (2.79) where the advanced selfenergy Af L;Rg is the Hermitian conjugate of the retarded selfenergy Rf L;Rg. The core of the problem lies in the calculation of the Green's functions of the semi in nite leads. The lead Green's function in an orthogonal localizedorbital Hamil tonian can be computed with an e cient principal layer method and we take the righthand side lead as an example. With increasing the distance between two localized orbitals, the overlap and in teraction matrix elements are smaller and smaller. Beyond some length, the elements could be considered as zero. So the semiin nite lead can be viewed as an in nite stack of principal layers with only nonzero nearestneighbor interactions. This corresponds to transforming the original system into a linear chain of principal layers. Within this approach, the Green's function equation for righthand lead can be expressed as 0 BBBBBBBBBBB@ " H00 H01 0 0 Hy 01 " H11 H12 0 0 Hy 12 " H22 H23 0 0 Hy 32 " H33 ... ... ... ... 1 CCCCCCCCCCCA 0 BBBBBBBBBBB@ G00 G01 G02 G03 G10 G11 G12 G13 G20 G21 G22 G23 G30 G31 G32 G33 ... ... ... ... 1 CCCCCCCCCCCA 31 = 0 BBBBBBBBBBB@ I 0 0 0 0 I 0 0 0 0 I 0 0 0 0 I ... ... ... ... ... 1 CCCCCCCCCCCA (2.80) which can be expanded explicitly as follows (" H00)G00 = I + H01G10 (" H00)G10 = Hy 01G00 + H01G20 (" H00)Gn0 = Hy 01Gn1;0 + H01Gn+1;0 (2.81) where Hnm and Gnm are the matrix elements of the Hamiltonian and Green's function between the layer orbitals, and we assume that in a bulk system H00 = H11 = and H01 = H12 = . In the iterative method proposed by LopezSancho et al.,[53, 54] this chain can be transformed such than the Green's function of an individual layer can be expressed in terms of Green's function of the preceding (or following) one. This is achieved by introducing transfer matrices T and T, de ned such that G10 = TG00 and G00 = TG10. The transfer matrix can be easily computed from the Hamiltonian matrix elements via an iterative procedure. T and T can be written as T = t0 + ~t0t1 + ~t0~t1t2 + + ~t0~t1~t2 tn T = ~t0 + t0 + t0t1~t2 + + t0t1t2 ~tn where ti and ~ti are de ned via the recursion formulas ti = (I ti1~ti1 ~ti1ti1)1t2i 1 ~ti = (I ti1~ti1 ~ti1ti1)1~t2i 1 32 with initial values t0 = (" H00)1Hy 01 ~t0 = (" H00)1H01 The process is repeated until tn; ~tn with arbitrarily small. In our calculations, we xed at 106. We always order the basis function from left to right in the matrix. Since there is only nonnegligible interaction between the conductor and the surface layer close to the conductor, only the topleft block G00 in Eq. (2.80) is needed. So the nonzero part of surface Green's function of right lead can be written as gR;00 = G00 = [" H00 H01T]1 : (2.82) and similarly the nonzero part of the surface Green's function of left lead can be written as gL;00 = h " H00 Hy 01 T i1 (2.83) In the case of general nonorthogonal orbitals, the above derivations still apply except the following substitutions " H00 ! "S00 H00 H01 ! ("S01 H01) where the matrix S represents the overlap between the localized orbitals. 33 CHAPTER 3 GRAPHENE NANORIBBONS 3.1 Introduction Graphene has attracted a great deal of research interest since its isolation via me chanical exfoliation[1, 2] because graphene exhibits many unusual properties such as massless fermions, minimum quantum conductance, and anomalous integer quantum Hall e ect occurring at halfinteger lling factors.[55] In addition graphene sheets show promise for future nanoelectronic applications, such as ultrafast transistors[56] because of the high carrier mobility.[3] Graphene materials are expected to play an important role in future electronic applications and even replace silicon some day. Twodimensional graphene, however, has a zero band gap with linear energy disper sion near the Fermi level. For practical application in the semiconductor industry a band gap could be induced by making graphene nanoribbons (GNRs) via a variety of methods such as lithographic patterning,[5] chemical vapor deposition,[6] Joule heating,[7] and unzipping the carbon nanotubes (CNTs).[8, 9] Graphene has two major types of high symmetry structures, namely those with edges that have either armchair or zigzag con gurations. Fujita, et al., [57{59] carried out tightbinding studies of these systems using them as model structures to study edge defects in CNTs. It was predicted that in zigzag GNRs there were strongly localized edge states near the Fermi level, suggesting electronic structure and low bias transport very sensitive to di erent edge passivations. The localized edge states were not expected in the armchair GNRs. The electronic energy gaps of armchair GNRs were calculated to have strong dependence on the ribbon widths. When we 34 started the project on GNRs in summer 2006, there were very few simulations done with density functional theory.[60, 61] Most studies were focused on the planar GNRs so far. there were only several investigations addressing the torsional deformation in GNRs. Bets and Yakobson,[62] using Molecular Dynamics Simulations, found that the narrow bare GNRs prefer the twisted structure rather than the planar geometry. Hod and Scuseria[63] reported the density functional theory calculations on twisted nanoribbons of nite length and found the gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is tunable by torsional deformation. In this chapter, we studied the electronic structures of planar zigzag and arm chair graphene nanoribbons. In addition, we investigated how the armchair graphene nanoribbons respond under applied twist about the ribbon axis. 35 3.2 Planar Zigzag Graphene Nanoribbons 3.2.1 Overview Figure 3.1 depicts the model structure of a zigzag graphene nanoribbon with ribbon width N = 6. Here the width N refers to the number of zigzag chains along the transverse direction. In each unit cell there are 2N C atoms and 2 H atoms. The in nite ribbon could be obtained by repeating the unit cell framed by dashed line in Figure 3.1 periodically along two sides with translation length l = 2.46 A. In our calculations we xed the CC bond length at 1.42 A and CH bond length at 1.08 A. We de ne the ribbon extension direction as the z axis perpendicular to the transverse direction. The righthanded rotation of screw operation about the z axis is zero in this case. We investigated ribbons with even values of N from 6  20. We used a 7s3p Gaussian basis for C, a 3s basis for H,[49] and 32 points were evenly sampled over the central Brillouin zone. For the zigzag graphene nanoribbons we carried out the nonspinpolarized (NSP), ferromagnetic (FM), and antiferromagnetic (AFM) calculations. The energetic trend on these states corresponding to di erent spin con gurations was analyzed. We found the antiferromagnetic states are energetically favorable. The zigzag ribbons exhibited strongly localized edge states near the Fermi level. Our results with Gaussian basis sets are consistent with others' results obtained by using tightbinding model and density functional theory simulations with numerical atomic orbital methods. [60, 64, 65] 36 Figure 3.1: Sample zigzag GNRs with ribbon width N = 6. The numbering of the unit cells is shown at the left side of the ribbon with only the unit cell labeled 0 shown in their entirety. Each unit cell is composed of 2N C atoms and 2 H atoms. The 2N C atoms in unit cell 0 are numbered as shown. The index of zigzag chains is indicated at bottom from left to right. 37 3.2.2 Energetic Results Figure 3.2 depicts the total energies divided by 2N and energy di erences between di erent spin con gurations as a function of ribbon width N. The corresponding numerical values are presented in Table 3.1. For the ribbon of same width, the antiferromagnetic state has the lowest energy and is energetically favorable. For example, for ribbon of N = 6, the energy of nonspinpolarized state is 0.20 meV higher than that of ferromagnetic state, which is further stabilized by about 0.032 meV in antiferromagnetic ribbon. As the ribbon width increases, the edges play a less important role and the total energies have more contributions from the bulk C atoms in the middle. So the energy di erence is decreasing with increasing ribbon width. When ribbon width N approaches 20, the energy di erence between FM and AFM states has already reduced to 0.00016 meV, suggesting for even bigger ribbons the preference to one of these two states is eliminated. In order to observe the AFM ground state, low temperature is needed since the thermal energy at room temperature is KBT = 26 meV with KB = 8:617343 105 eV/K and T = 300K, which will be able to excite the system and mix the various spin con gurations. 38 (a) (b) Figure 3.2: (a) Total energy di erences divided by 2N between nonspinpolarized and ferromagnetic states in zigzag GNRs as a function of ribbon width N. (b) Total energy di erences divided by 2N between ferromagnetic and antiferromagnetic states in zigzag GNRs as a function of ribbon width N. Table 3.1: Total energies divided by 2N and energy di erences between di erent spin con gurations of zigzag GNRs with various ribbon widths width N NSP (eV) FM (eV) AFM (eV) NSPFM (meV) FMAFM (meV) 6 37.4358823900 37.4360743104 37.4361059109 0.191920458334 0.0316004916669 8 37.4147775171 37.4149432678 37.4149598679 0.165750693746 0.0166001437520 10 37.4021060503 37.4022471735 37.4022567277 0.141123230001 0.0095541400071 12 37.3936847190 37.3938081844 37.3938136308 0.123465429169 0.0054464166652 14 37.3876676725 37.3877772624 37.3877815346 0.109589875002 0.0042721500009 16 37.3831570254 37.3832537549 37.3832574688 0.096729465625 0.0037139687521 18 37.3796450210 37.3797317062 37.3797335467 0.086685177777 0.0018405499986 20 37.3768360202 37.3769151528 37.3769153107 0.079132679993 0.0001578675040 39 3.2.3 Nonspinpolarized LDF Band Structures In the nonspinpolarized calculations, we assumed the spinup and spindown states are degenerate. The typical calculated band structures are presented in Fig ures 3.3, 3.4, and 3.5, corresponding to the ribbon of widths N = 6, 8, and 10, respectively. They all exhibited nearly at bands crossing the Fermi level. The at region in Brillouin zone is ranging from about = 2 =3 to the zone edge, consis tent with the tightbinding calculations.[57] Because the density of states is inversely proportional to the slope of the energy dispersion, the atness in lowerlying bands generates high density of states around the Fermi level. We will show these states are localized on the ribbon edges by plotting the orbital density in real space. 40 Figure 3.3: Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. 41 Figure 3.4: Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. 42 Figure 3.5: Calculated nonspinpolarized band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. 43 Figure 3.6 depicts the orbital density of zigzag GNRs of width N = 6. The wavevectors are sampled at 0.625 , 0.75 , 0.875 , and in the Brillouin zone. For state at the zone edge with = , the wavefunctions are completely localized at the edge C sites. When is moving away from the zone edge, the orbital density extends and decays towards the ribbon center. If = 0:625 , the states become extended across the transverse direction. Our simulation results are consistent with the previous tightbinding prediction in which the neargap state is completely localized at the edge sites when = , and starts to gradually penetrate into the inner sites as deviates from , reaching the extended state at = 2 =3.[57] Figure 3.6: Calculated HOMO band orbital densities of zigzag GNRs of width N = 6, calculated at wavevectors (a) , (b) 0.875 , (c) 0.75 , and (d) 0.625 , respectively. Calculated orbital densities for the LUMO band are given for wavevectors (e) , (f) 0.875 , (g) 0.75 , and (h) 0.625 , respectively. The isovalue is taken to be 0.008 for all plots. 44 3.2.4 Ferromagnetic LDF Band Structures If we do not apply constrains of spin degeneracy, spinpolarized calculations gave us ferromagnetic states. The typical band structures are shown in Figures 3.7, 3.8, and 3.9, corresponding to the ribbon of widths N = 6, 8, and 10, respectively. In this case the spin degeneracy was lifted. Compared to the nonspinpolarized results, the at bands and the corresponding high density of states around Fermi level were not observed. Instead two peaks in density of states appeared with centers at energy points 0.3 eV above and below the Fermi level. For electron excitation, the spin down states dominate the density of states. For hole excitation, the spinup states contribute more. 45 Figure 3.7: Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. 46 Figure 3.8: Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. 47 Figure 3.9: Calculated ferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. 48 3.2.5 Antiferromagnetic LDF Band Structures We also carried out the calculations with constraint of applying inversion sym metry about the z axis in the exchange correlation potential. From this, we got the antiferromagnetic ground states. The typical band structures are depicted in Fig ures 3.10, 3.11, and 3.12, corresponding to the ribbon of widths N = 6, 8, and 10, respectively. In Figure 3.13 we depicts the energy gaps as a function of ribbon width N. The corresponding numerical values of gaps are listed in Table 3.2. The ribbons exhibited direct gap nature. However the gaps open up at 0:70 rather than the zone center or edge. When increasing ribbon width, the band gaps are decreasing monotonically. Our results are consistent with that of Son, et al.[64] In their studies, the e ect of electric eld was also addressed. When an electric eld is applied along the transverse direction, spinup and spindown electrons respond in di erent way. The spinup electrons are moving toward the Fermi level, leading to smaller gap even gap closure for big enough eld. The spindown electrons are moving away from the Fermi level and will not play a role in the electron transport within some voltages. The response of electrons under applied electric eld renders the antiferromagnetic graphene nanoribbons a half metal. 49 Figure 3.10: Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 6. (X) corresponds to the center (edge) of the rst Brillouin zone. 50 Figure 3.11: Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 8. (X) corresponds to the center (edge) of the rst Brillouin zone. 51 Figure 3.12: Calculated antiferromagnetic band structure and density of states of zigzag GNRs of width N = 10. (X) corresponds to the center (edge) of the rst Brillouin zone. 52 Figure 3.13: Magnitude of the HOMOLUMO antiferromagnetic gap of the zigzag GNRs as a function of ribbon width N. Table 3.2: Energy gaps of zigzag GNRs with various ribbon widths Ribbon width N Energy gap (eV) Ribbon width N Energy gap (eV) 4 0.5732 14 0.3008 6 0.4972 16 0.2790 8 0.4341 18 0.2618 10 0.3901 20 0.2539 12 0.3404 53 3.2.6 Summary We calculated the electronic structures of zigzag graphene nanoribbons for dif ferent spin con gurations. It was found that there were localized edge states, which make the electronic properties of zigzag GNRs tunable by edge chemistry. [16, 66, 67] The zigzag GNRs exhibited ground states with two ferromagnetic edges antiferomag netically coupled. With transverse electric eld applied, the zigzag GNRs were found to be half metals,[64] which provides a suitable platform for the spintronics. This halfmetallicity is further enhanced by edgeoxidization[68] and edge modi cations with NO2 and CH3 terminations on opposite sides.[69] 54 3.3 Planar Armchair Graphene Nanoribbons 3.3.1 Overview Nakada et al.,[58] with nearestneighbor tightbinding model, predicted that the energy gaps of armchair GNRs are strongly dependent on the ribbon width N. If the ribbon width satis es N = 3m + 2 with m being an integer, the ribbon is metallic of zero gap. Otherwise, the semiconducting energy gaps monotonically decrease with increasing ribbon widths. Our rst principles results with local density functional (LDF) do not agree with the TBM prediction. Especially in our simulations, there is no metallic ribbon. In order to elucidate the discrepancy between TBM and LDF results, we introduced the thirdnearestneighbor interaction as perturbation correc tion to the nearestneighbor TBM and show that the longer range interactions play an important role in describing the band gaps of armchair graphene nanoribbons. 55 3.3.2 Model Structure Figure 3.14 depicts the model structure of armchair graphene nanoribbons with ribbon width N = 7 and 8. Here the width N is referred to as the number of dimer lines along the transverse direction. Depending on the way of constructing the whole in nite ribbon from the basic helical cell as framed by dashed line, the armchair graphene nanoribbons are classi ed into two categories : symmetric and staggered ribbons. For symmetric ribbons as shown in Figure 3.14 (a), the translation length l and rotation angle is 4.26 A and 0 degree, respectively, and in each unit cell there are 2N C atoms and 4 H atoms. For staggered ribbons as shown in Figure 3.14 (b), the translation length l and rotation angle is 2.13 A and 180 degrees, respectively, and each unit cell contains N C atoms and 2 H atoms. In our calculations the CC bond length and CH bond length are xed at 1.42 A and 1.08 A, respectively. 56 (a) (b) Figure 3.14: (a) Symmetric armchair GNRs with ribbon width N = 7. Each unit cell is composed of 2N C atoms and 4 H atoms. (b) Staggered armchair GNRs with ribbon width N = 8. Each unit cell is composed of N C atoms and 2 H atoms. For (a) and (b) the numbering of the unit cells is shown at the left edge of each gure with only the unit cells labeled 0 shown in their entirety. The C atoms in unit cell 0 are numbered as shown. The index of dimer lines is indicated at bottom from left to right. 57 3.3.3 Local Density Functional Results We carried out the local density functional calculations on armchair graphene nanoribbons with various widths ranging from 4 to 24. We used a 7s3p Gaussian basis for C, a 3s basis for H, and 32 points were evenly sampled over the central Brillouin zone. Figure 3.15 depicts the band structures for armchair GNRs of widths N = 7, 8, and 9. and in Figure 3.16 we present the energy gap as a function of ribbon width N using the gap values listed in Table 3.3. In the calculations of staggered armchair GNRs, we were making use of the helical symmetry of rotation. The presented band structures, however, have already been folded to be within the regular Brillouin zone corresponding to the translational operation. First, all ribbons are semiconductors and exhibit direct band gap with HOMO and LUMO occurring at point, irrespective of the ribbon width in our study. On the other hand, the gap decreases with increasing ribbon widths. However, the variations in energy gaps have alternating pattern rather than monotonic decrease. We have relative relation of Egap; 3m+1 > Egap; 3m > Egap; 3m+2 with m being integer. 58 Figure 3.15: Calculated band structures of armchair GNRs of various widths. (a), (b), and (c) correspond to ribbon with N = 7, 8, and 9, respectively. (X) indicates the center (edge) of the rst Brillouin zone. 59 6 9 12 15 18 21 24 N 0 0.5 1 1.5 2 2.5 Egap (eV) LDF mod (N + 1, 3) = 2 LDF mod (N + 1, 3) = 1 LDF mod (N + 1, 3) = 0 Figure 3.16: LDF energy gaps of the armchair GNRs as a function of width N. Table 3.3: Band gaps Egap of armchair GNRs with ribbon width N N Egap(eV) N Egap(eV ) N Egap(eV) N Egap(eV) 4 2.4106 10 1.0734 16 0.6874 22 0.5055 5 0.2982 11 0.1512 17 0.1006 23 0.0752 6 1.0158 12 0.5692 18 0.3953 24 0.3030 7 1.4872 13 0.8375 19 0.5821 8 0.2012 14 0.1210 20 0.0868 9 0.7298 15 0.4669 21 0.3437 60 3.3.4 Consideration of Longer Range Interactions Tightbinding models were frequently used to analyze the electronic properties of semiconductors. With only nearestneighbor interactions considered tightbinding models have been very successful in describing many properties of singlewall carbon nanotubes (SWCNTs).[22, 70{83] The geometries of SWCNTs and GNRs have great similarity because both of them can be obtained from graphene by rolling over a chiral vector or cutting along speci c edges. This similarity, along with the success of TBM in SWCNTs, allow us to expect the applicability of TBM in studying GNRs. However, the LDF results are quite di erent from the previous predictions based on nearestneighbor tightbinding model. In the tightbinding model, if the ribbon width N satis es N = 3m + 2, then the ribbon is metallic with zero gap at point where and bands crossing occurs. In addition, the band gaps of other ribbons are decreasing monotonically with increasing ribbon width. The discrepancy motivated us to think over the validity of nearestneighbor tight binding model in describing the electronic property of armchair GNRs. In the previous calculations, only nearestneighbor interactions are included. The gap openings in armchair GNRs of width N = 3m+ 2 from our rstprinciples calculations, however, indicated that there were some interactions which are important but ignored in the tightbinding description. Since the armchair GNRs can be obtained by cutting graphene along speci c direction, we could apply proper boundary condition on the energy dispersion for 2D graphene, discretize the wavevector along the transverse direction, and nd the energy dependence on the 1D wavevector in armchair GNRs. However, we are more concerned with the band gap instead of the energy values in whole Brillouin zone and will take a di erent method. The underlying procedure is described below. 61 The ribbons are direct band gap materials with HOMO and LUMO states occur ring at k = 0 (see Figure 3.15). Therefore, to study the gap dependence on ribbon width we need only to consider the eigenstates of nearestneighbor Hamiltonian ^H 0 at k = 0. In the following discussion, we make use of the regular Bloch theorem to construct the wavefunctions based on the translational symmetry down the ribbon axis. So each unit cell contains two zigzag chains for symmetric and staggered arm chair GNRs. We can apply the the boundary condition with vanishing nodes at the ends of zigzag chains. Then the eigenstates of ^H 0 for the armchair GNRs with ribbon width N are given by j p >= 1 p Ncell 1 p N + 1 X l XN n=1 sin( pn N + 1 )(jn; l > +sjn + N; l >); (3.1) where 1 p N, jn; l > denotes the jpz > orbital associated with the nth C atom in the unit cell labeled by l in the notation of Figure 3.14, and the s = 1 (1) indicates the bonding (antibonding) state between the wavefunctions belonging to neighboring zigzag chains. The calculated energy gap with nearestneighbor interaction V1 only is Egap = 2V1[2 cos( p N + 1 ) + s] (3.2) where we have : if mod(N+1,3) = 0, then p = (N +1)=3 and s = 1; if mod(N+1,3) = 1, then p = (2N + 3)=3, and s = 1; and, if mod(N+1,3) = 2, then p = (N + 2)=3, and s = 1. Then the band gap strongly depends on the ribbon width. The system is metallic if N = 3q 1, and is semiconducting otherwise, see Figure 3.17 in which V1 = 2:6 eV is assumed. Except the nearestneighbor interaction, longer range interactions between C atoms are neglected and may modify the band structure results of TBM. These longer range interactions can be incorporated into the TBM by adding ^U V2 and ^U V3 to ^H 0, the Hamiltonian in TBM, where ^U V2 (^U V3) includes all secondnearestneighbor (thirdnearestneighbor) interactions V2 (V3) between C atoms in the armchair GNRs. Because V2 and V3 are small in magnitude compared to V1, ^U = ^U V2 + ^U V3 can be 62 6 9 12 15 18 21 24 N 0 0.5 1 1.5 2 2.5 Egap(eV) TBM mod (N + 1, 3) = 2 TBM mod (N + 1, 3) = 1 TBM mod (N + 1, 3) = 0 Figure 3.17: Magnitude of the HOMOLUMO gap of the N armchair GNRs from the TBM calculations. treated as perturbation to ^H 0. Also, because ^U does not break the translational symmetry of the ribbon, it couples only states with the same k, which allows us to consider the states at k = 0 only. We carried out the timeindependent perturbation calculations to rst order with the V2 and V3 as perturbation terms to V1 and found the new expression for Egap given by Egap = 2V1[2 cos( p N + 1 ) + s] + s 2V3 N + 1 [3 + N + 2N cos( 2p N + 1 )] (3.3) where we have : if mod(N+1,3) = 0, then p = (N +1)=3 and s = 1; if mod(N+1,3) = 1, then p = (2N + 3)=3, and s = 1; and, if mod(N+1,3) = 2, then p = (N + 2)=3, and s = 1. The rst term on the righthand side of Eq. (3.3) is the gap arising from ^H 0. The second term gives the lowest order correction because of ^U . Although ^U V2 contributes to ^U, V2 does not appear in Eq. (3.3) because to rst order it shifts the HOMO and LUMO levels by the same amount 3NV2=(N + 1). If N + 1 = 3m, the rst term on the righthand side of Eq. (3.3) vanishes leaving 63 only the second term which reduces to Egap = 6V3 N + 1 : (3.4) Eq. (3.4) can be used to leastsquares t to the LDF data for N +1 = 3q to determine V3 and then with V3 xed at this value V1 can be determined by using Eq. (3.3) to leastsquares t to the LDF data for N +1 6= 3q. Implementing this procedure yields the physically reasonable results: V1 = 3:2 eV and V3 = 0:3 eV, which provide an excellent t to the LDF results for all N as shown in Figure 3.18. 6 9 12 15 18 21 24 N 0 0.5 1 1.5 2 2.5 Egap (eV) LDF mod (N + 1, 3) = 2 LDF mod (N + 1, 3) = 1 LDF mod (N + 1, 3) = 0 TBM mod (N + 1, 3) = 2 TBM mod (N + 1, 3) = 1 TBM mod (N + 1, 3) = 0 Figure 3.18: LDF results for Egap compared to the corresponding tightbinding model (TBM) results obtained from Eq. (3.3) with V1 = 3:2 eV and V3 = 0:3 eV. 64 3.3.5 Summary In our LDF simulations we xed all CC bond lengths at 1.42 A and found that the armchair GNRs can be separated into three families according to the band gap dependence on ribbon width. Then we proposed an explanation by using TBM with perturbation correction from the second and thirdnearestneighbor interactions. The nonzero semiconducting gaps in GNRs of width N = 3q + 2 and three families of armchair GNRs were also observed by Son, et al.[84] In their calculations, they carried out geometry optimization and found that the CC distances on edges decrease leading to 12% increase of the hopping integrals between pz orbitals of edge C atoms. After taking matrix elements change into account, they also got a TBM results which agree with their LDA calculations. Gunlycke, et al.[85] found that both edge contraction and thirdnearestneighbor interactions play an important role in producing good tightbinding band structures which t the rstprinciples results of optimized GNRs. Experimentally, great progress has been made in preparing GNRs with widths less than 10 nm,[86, 87] which makes our theoretical calculation of more fundamental and practical signi cance. 65 3.4 Twist E ect in Electronic Properties of Armchair Graphene Nanoribbons 3.4.1 Overview In the study of planar armchair GNRs, we show that the band gap could be engineered by controlling the ribbon width. The band gaps of nanostructures also could be tuned by other ways such as edge passivation, doping, and strain. Strain is a powerful approach to modify the electronic structure of nanomate rial to improve device performance. The carrier mobility of silicon nanowires could be signi cantly enhanced by applied strain.[88] Firstprinciples calculations have shown that for smaller diameters the band structures and gaps can be signi cantly modi ed.[89, 90] The strain e ect in graphene also attracted much attention. Ni, et al.[91] predicted a bandgap opening of 300 meV for graphene under 1% uniaxial tensile strain. Teague, et al.[92] found that the strain could modulate the local conductance of graphene. Lu and Guo,[93] using a tightbinding model, investigated armchair GNRs and found uniaxial weak strain changes the band gap in a linear fashion and large strain results in periodic oscillation of the band gap. Sun, et al.,[94] using GGA exchange correlation potentials, predicted that the similar oscillating behavior with linear change in each period. For the outofplane strain, Hod and Scuseria[63] studied the torsional deformation e ect in graphene nanoribbons of nite length. With traditional band structure codes using translational symmetry, it is hard even almost impossible to study the torsional deformation for in nitelong ribbon structure because of the huge unit cell and hence the Hamiltonian matrix of enormous size. Our HENS parallax code takes advantage of the helical symmetry and is well suited to the calculation of electronic structures 66 of the twisted quasionedimensional structures. In this section we investigated the band gaps of armchair GNRs as a function of twist angle. We found that the electronic properties of armchair GNRs are highly sensitive to the applied twist, indicating the potential application of armchair GNRs as building blocks in nanoelectromechanical devices. 67 3.4.2 Computational Approach We investigated the armchair GNRs of widths N = 8; 10; 12; :::; 18 in this study. The representative structure of the armchair GNRs is shown in Figure 3.19. The whole ribbon can be constructed with a unit cell consisting of one zigzag chain as framed by dashed lines in Figure 3.19(a) and a screw operation acting on that unit cell, where the screw operation combines a translation l down the ribbon axis with a righthanded rotation about that same axis. In our calculations we xed l at 2.13 A. For the planar structure, the is 180 . For the twisted armchair GNRs as shown in Figure 3.19(b), the torsional deformation of armchair GNRs can be represented by the relative rotation angle , de ned as = 180 . The planar ribbons have symmetry operations corresponding to a re ection perpendicular to the ribbon axis and passing through the center of CC bonds parallel to the ribbon axis. When applying the twist, the re ection symmetry of armchair GNRs is broken and only C2 symmetry is left. All edge C atoms are passivated by H atoms to remove the dangling bands near the Fermi level. Following conventional notation, the ribbon width of armchair GNRs is de ned as the number of dimer lines along the ribbon forming the width as in the preceeding section. For example, the ribbon width in Figure 3.19 is 8. The LDF band structures were obtained by using 32 discrete points evenly sampled in the central Brillouin zone and a 631G* basis set. [95] Although the LDF is well known to underestimate the band gap, we are primarily concerned with the relative change of the energy gaps under di erent twists. 68 (a) (b) Figure 3.19: (a) Armchair GNR with width N = 8. Each helical unit cell is composed of one zigzag chain terminated by H atoms. The N C atoms in a helical cell are numbered as shown. The bonds are labeled as a1a11. (b) Axial view of 10 twisted armchair GNR with width N = 8. 69 3.4.3 Band Structures and Neargap Wavefunctions We depict the electronic structures for armchair GNRs with width N = 8 under a series of twist angles in Figure 3.20, in which (a), (b), (c), and (d) correspond to twist angle equal to 0 , 8:5 , 9:0 , and 10 , respectively. In this ribbon, the top of the valence band and the bottom of the conduction band occur at the center of the one dimensional Brillouin zone, resulting in a direct band gap at the zone center point. The band gap for ribbon of width N = 8 is 0.40 eV for = 0. When we apply on the ribbon torsional strain about the ribbon axis, the HOMO is raised while the LUMO is lowered if < 9:0 . therefore the band gap is reduced in this stage. In general two bands cannot cross because they belong to the same irreducible representation for wave vectors not at the center or edge of the Brillouin zone. However, gap closure (accidental degeneracy) might take place at a point of high symmetry in the Brillouin zone such as the zone center or end. After the critical twist angle, the HOMO and LUMO do the inverse way, opening the band gap. In addition to the band gap change, we also observed the twist e ect on the characteristics of the bonding or antibonding states in HOMO and LUMO. Figure 3.21 depicts the wave functions of HOMO and LUMO under four di erent twist angles corresponding to that of Figure 3.20. According to the bonding states between the pz orbitals of neighboring zigzag chains, there are two kinds of bonds: bonding states (BS) and antibonding states (AS). When the twist angle is smaller than 9:0 , the HOMO is BS state while LUMO is AS state. When the twist angle gets bigger, the characteristics of HOMO and LUMO exchanged. The energy of BS and AS states respond di erently to the geometry change. Increasing bond length will lower the energy of AS state and raise the energy of BS state. 70 Figure 3.20: Electronic band structure of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . 71 Figure 3.21: The LUMO and HOMO wave function of armchair GNRs (width N = 8) twisted by (a) 0 , (b) 8:5 , (c) 9:0 , and (d) 10:0 . The ribbons are vertically oriented. 72 3.4.4 Structure Parameters In our simulations we carried out geometry optimization for all structures. Fig ure 3.22 depicts the variation of bond lengths in armchair GNRs of width N = 8 as a function of the twist angle. Because of symmetry in this ribbon, there are only eight nonequivalent bonds labeled by a1, a4, a7, a10 along the ribbon extension direction and by a2, a3, a5, a6 along the transverse direction. By symmetry the bonds a8, a9, and a11 are equivalent to a5, a3, and a2, respectively. Depending on the position, di erent bonds in a di erent way respond to the torsional deformation. The bonds along the transverse direction change a little bit with maximum change about 0.01 A. The bonds connecting neighboring zigzag chains exhibit bigger change. Generally the change is increasing as bonds location move towards the ribbon edges. But the biggest bond length change of 0.045 A occurs on a10 instead of a1 which is furthest from the ribbon center. As the twist angle approaches 9 , some bonds exhibit abrupt change, which corresponds to the transition between the antibonding and bonding HOMOs. In other ribbons of di erent widths, we observed similar bond length change as a function of twist angle and the bond equivalent to a10 has the biggest change. The bonds underlying bigger change play more important role in determining the shift of one electron energy of HOMO and LUMO. So only from the bonding (antibonding) states between interchain bonds we could acquire some idea about how the twist a ects the electronic band gap of armchair GNRs. Tightbinding model could be used to carry out this analysis. 73 (a) (b) Figure 3.22: CC bond lengths of the twisted armchair GNRs with width N = 8 as a function of twist angle . 74 3.4.5 Tightbinding Model We start with the wavefunction as introduced in previous section, 1 p Ncell(N + 1) X l XN n=1 sin np N + 1 (jn; li + sjn + N; li) where 1 p N, and jn; li denotes the pz orbital associated with nth C atom in the unit cell labeled by l in the notation of Figure 3.19. And the eigen energy at can be expressed as follows: E(s; p) = 2V1[cos( p N + 1 ) + s]; where s = 1 indicates the BS and s = 1 indicates the AS. In this nearestneighbor tightbinding model, the ribbon of N = 3m + 2 is metallic with HOMO and LUMO touching each other. We cannot distinguish the HOMO and LUMO. In this case, we introduced the thirdnearestneighbor interactions as perturbation. The perturbation shifted the energy by 3=(N + 1)V3 upward for LUMO and the same amount down ward for HOMO. The third nearestneighbor interactions are required to open the gap to give results consistent with our LDA results. After introducing these interac tions, we determine that LUMO is at p = (N + 1)=3 and antibonding, HOMO is at 2(N + 1)=3 and bonding. Because in some ribbons, the second HOMO (HOMO1) and second LUMO (LUMO+1) are moving towards the Fermi level under twist and will play a role in determining the band gap beyond certain twist angle. We also analyzed the HOMO1 and LUMO+1 states and present the results in Table 3.4. 75 Family HOMO LUMO HOMO1 LUMO+1 N = 3m # +" +" # N = 3m + 1 +" # # +" N = 3m + 2 +" # # +" Table 3.4: Binding states of edge C pz orbitals HOMO and LUMO. " (#) indicates the energy point shifts upward (downward) with increasing twist angle. + () indicates the bonding (antibonding) states. 76 3.4.6 Band Gap Change Under Twist In Figure 3.23 and Table 3.5 we present the local density functional results for the variations of energy gaps of three family structures with di erent widths (N = 3m, 3m+1, 3m+2, where m = 2; 3; 4) as a function of the twist angle . In combination with the tightbinding model in previous section, we will discuss in details how the twist a ects the energy gaps in armchair graphene nanoribbons. Depending on the bonding (antibonding) states of near Fermi level bands, the response of the armchair GNRs to the twist can be classi ed into two categories. One is of width N = 3m, and the other two belong to the second family. For N = 3m, the energy gap will go through three stages with increasing twist. In the rst stage, the bonding LUMO goes up and antibonding HOMO goes down, resulting in the raising of band gap. The HOMO1 is raising and LUMO+1 is lowering. After certain twist angle, HOMO1 and LUMO+1 will take over to determine the band gap. This begins the second stage, in which the band gap is reduced with increasing twist angle. When the twist is increased further, the HOMO1 and LUMO+1 may touch and exchange bonding characteristics and then follow their previous way. This will start the third stage, in which the band gap is increased again. For N = 12 we only observed the rst stage. For N = 18, we observed all these three stages because the reduced quantum con nement and bigger change in length of bonds parallel to the ribbon axis resulted from the bigger ribbon width. For N = 3m + 1 and N = 3m + 2, the HOMO1 and LUMO+1 are moving away from the Fermi level and will not play a role in determining the band gap. Within the twist angle studied, we only have to consider the HOMO and LUMO orbitals. Just as we have seen in armchair GNRs with width N = 8, the band gap will have two stages. The rst stage is decreasing the band gap and the second would be increasing 77 the gap. In N = 10, we only have the rst stage because of the relatively bigger band gap and not big enough ribbon width. Generally the turning point is moving toward smaller twist angle with increasing ribbon width for same family. At each turning point corresponding to the minimum band gap, the gap closure may occur because of the additional symmetry at the zone center or edge. We can not have real band crossing because the eigen vectors of bands belong to the same irreducible representation group. Our bonding and antibonding analysis applies for other strain e ect in the arm chair GNRs as well. Our tight binding analysis are qualitatively consistent with the previous calculations on uniaxial strain e ect on armchair GNRs by others. But the linear relation by TB model[93] or density functional theory calculations[94] is not observed in our calculations. This may be due to the hybridization between the s, px, py, pz orbitals coming from the twist between neighboring zigzag chains. 78 Figure 3.23: Variation of energy gaps of armchair GNRs as a function of twist angle . (a), (b), and (c) represent the armchair GNRs of width N = 3m, 3m + 1, and 3m + 2, respectively. 79 Table 3.5: Band gaps (eV) as a function of twist angle (Deg) in armchair graphene nanoribbons with various widths N Twist angle N=8 N=10 N=12 N=14 N=16 N=18 0.0 0.3958 1.2364 0.4648 0.2830 0.8166 0.2795 0.5 0.3904 1.2382 0.4677 0.2800 0.8170 0.2850 1.0 0.3842 1.2281 0.4732 0.2717 0.8159 0.2949 1.5 0.3860 1.2202 0.4816 0.2613 0.7983 0.3147 2.0 0.3727 1.2125 0.4951 0.2406 0.7764 0.3415 2.5 0.3632 1.1979 0.5109 0.2268 0.7486 0.3744 3.0 0.3579 1.1869 0.5276 0.1972 0.7125 0.4146 3.5 0.3456 1.1658 0.5498 0.1651 0.6714 0.4619 4.0 0.3313 1.1440 0.5738 0.1276 0.6239 0.5165 4.5 0.3116 1.1226 0.6051 0.0868 0.5839 0.5737 5.0 0.2909 1.0954 0.6329 0.0217 0.5215 0.6404 5.5 0.2739 1.0688 0.6642 0.0697 0.4568 0.7115 6.0 0.2541 1.0399 0.7004 0.1190 0.3907 0.7197 6.5 0.2309 1.0085 0.7374 0.1732 0.3200 0.6293 7.0 0.2062 0.9763 0.7749 0.2302 0.2435 0.5341 7.5 0.1821 0.9423 0.8180 0.2885 0.1625 0.4341 8.0 0.1541 0.9063 0.8610 0.3533 0.0780 0.3273 8.5 0.1274 0.8689 0.9070 0.4214 0.0761 0.2141 9.0 0.1990 0.8293 0.9532 0.4891 0.1746 0.0994 9.5 0.0505 0.7895 1.0027 0.5673 0.2764 0.0702 10.0 0.0815 0.7492 1.0560 0.6476 0.3847 0.1582 80 3.5 Summary and Conclusions In this chapter we studied two kinds of GNRs: zigzag and armchair GNRs. We found that the zigzag GNRs have magnetically ordered insulating ground states that are ferromagnetically coupled along each edge and antiferromagnetically coupled across the edges. And the AFM band gap decreases with increasing ribbon width. The electronic states near Fermi level are localized at edge carbon sites, suggesting the edge sensitivity to passivating functional groups and therefore zigzag GNRs could be used in chemical sensor applications. Our LDF results with Gaussian basis sets agree with previous calculation by using plane wave basis, numerical atomic orbitals, and also are consistent with the tightbinding calculations on the zigzag GNRs. armchair GNRs do not have spinpolarized ground state and localized edge states. The band gaps of armchair GNRs are strongly dependent on the ribbon width. The complex alternating relation between gaps and ribbon width is di erent from the nearestneighbor tight binding prediction. The introduction of thirdnearestneighbor interaction across the hexagons resolved this discrepancy. Besides the planar GNRs, we investigated the twisted armchair GNRs and found the response to the applied torsional deformation could be classi ed into two categories. The tunable band gap upon applied twist is very useful in nanomechanical devices. It is wellknown that the LDF underestimated the band gaps. Yang, et al.[96] pre sented calculations of the quasiparticle energies and band gaps of planar graphene nanoribbons carried out using a rstprinciples manyelectron Green's function ap proach within the GW[97{100] approximation. The selfenergy corrections are di er ent for ribbons of di erent widths because of the screening e ect of various magni tudes. However, the characteristics of three families in band gaps remain the same in armchair graphene nanoribbons. 81 CHAPTER 4 SURFACE PASSIVATION EFFECTS IN SILICON NANOWIRES 4.1 Overview In recent years Si nanowires (SiNWs) have attracted much attention because of their potential applications in electronic, thermoelectric, and optoelectronic devices. They can be used as building blocks of nanoscale electrical devices such as eld e ect transistors (FETs).[101{104] They are capable of in ating 4 times their normal size when absorbing lithium ions, which enables the new battery to hold 10 times the charge of existing lithiumion batteries.[10] SiNWs can be used as thermoelectric materials because of their low thermal but high electrical conductivity.[12, 13] SiNWs also have potential application as solar cells to convert light into electricity. [11, 105{ 110] The broad range of potential applications of SiNWs makes it critically important to understand how to tailor their electronic properties. The electronic properties of SiNWs can be modi ed by varying a range of structural properties such as surface passivation,[101] doping,[102] and how the nanowires are mechanically processed.[111] Because of the enhanced surfacetovolume ratio of the SiNWs compared to bulk silicon, the surface of the SiNWs is of special status and the electronic properties of the SiNWs are strongly dependent on the characteristics of the SiNW surfaces. Cui, et al.[101] showed that exposure to 4nitrophenyloctadecanoate or tetraethylammo nium bromide can improve the conductance and ono ratios in oxidized SiNWbased FETs. Haick, et al.[112] achieved the methyl passivation on SiNWs via a twostep chlorination/alkylation method and showed that the SiNWs have enhanced air sta 82 bility and high hole mobility. Although the question of how the dangling bonds are passivated at the silicon nanowire surface is of great importance, most of the theoretical studies have been carried out on Hterminated[99, 113{119] and bare SiNWs,[120{122] focusing on quan tum con nement e ects, orientation, and crosssection dependence of the electronic properties of SiNWs, doping e ects, and surface reconstructions. There are very few investigations focusing on the surface e ects. For examples, Blase, et al.[123] calcu lated the e ects of a single alkyl chain on the zerobias conductance of SiNWs and found that the Landauer conductance near Fermi level is almost not a ected. Nolan, et al.,[124] using a combination of plane wave basis sets and pseudopotentials, ob served that passivation with OH or NH2 reduces the band gap of h100i SiNWs. Aradi, et al.[125] reported that the OH passivation also can induce the band gap reduction of h110i SiNWs. Leu, et al.[126] calculated the band structure of SiNWs with halogens surface substituents passivating the SiNW surfaces, and found that passivation with more weakly interacting surface species appears to lead to more marked reductions in the SiNW band gaps. The band gap reduction comes from the weakly interacting surface species which can not pull surface states of bare SiNWs out of midgap. In our studies we choose the silicon nanowires along h110i and h100i as objects of study based on the experimental observation on the relations between the wire diam eter and orientations. When synthesizing silicon nanowires using gold nanocluster catalyzed 1D growth method, Wu, et al.,[127] found that the nanowires of diameter less than 10 nm prefer the h110i directions. Most bigger nanowires were h112i and h111ioriented. Schmidt, et al.,[128] using epitaxial vaporliquidsolid (VLS) growth method, observed similar orientation dependence on the diameter. To my knowl edge, no h100i SiNWs were obtained with the nontemplated methods. Because only small silicon nanowires can be treated with the density functional theory, most com putational simulations are focused on h110i SiNWs such that the calculation results 83 could be compared with experimental measurements. The h100ioriented nanowires are highly desirable since the conventional CMOS microelectronics are built on (100) silicon wafer. The h100i SiNWs can be used to fabricate the vertical eld e ect tran sistors on the (100) wafer. Researchers, with anodic alumina membranes template directed growth technique, have obtained h100i silicon nanowires of diameter ranging from 60 nm to 200 nm.[129{131] In this chapter we rst present the rstprinciples study of the electronic band structures of SiNWs oriented along h100i and h110i directions with the surfaces pas sivated by hydrogen, hydroxyl, and methyl substituent groups. For convenience, sometimes we simply refer to the SiNWs of these three di erent surface passivations as HSiNWs, OHSiNWs, and CH3SiNWs, respectively. We observed the reduced band gaps with increasing wire diameter irrespective of passivating groups, mani festing the quantum con nement e ects. Band gap reductions in OHSiNWs with reference to HSiNWs are also observed in our simulations, in agreement with others' results.[124, 125] CH3 not only reduces the band gap but also leads to indirect band gaps for all h100i SiNWs studied. The passivations on surface not only can change the band gap nature and gap value, but also have e ect on the e ective mass of carriers in SiNWs. Our simulations suggest that passivation with CH3 surface substituents sub stantially increases the electron e ective mass for the h100i wires, while h110i SiNWs have small electron and hole e ective masses for all three passivations studied. The change on electronic band structures by surface passivations motivated us to study their e ect on the transport conductance in SiNWs. As a preliminary study, we investigated the electron conductance of a h110i SiNW having a defect region with hydroxyl passivation rather than hydrogen passivation. 84 4.2 Computational Approaches The Si nanowires studied herein were initially constructed from silicon in the di amond structure using the bulk lattice constant, a = 5.43 A. All h100i nanowires were constructed with square crosssections and f110g faces, and all h110i nanowires were constructed with diamond crosssections and f111g faces. We carried out geom etry optimization for all structures. The optimized structure for the h100i and h110i SiNWs are depicted in Figure 4.1. The h100i wires calculated have a 4fold screw symmetry about the wire axis with translation l = a=4 and screw angle = =2, and the h110i wires have a 2fold screw symmetry about the wire axis with translation l = p 2a=4 and screw angle = . The di erent crosssections of the h100i and h110i oriented SiNWs make a direct comparison of the radial scale of the wires somewhat arbitrary. We have chosen to de ne an e ective nanowire diameter, d, given by the expression d = p a3N=2 l, where a is the lattice constant of bulk crystalline silicon, N is the number of Si atoms in the unit cell, and l is the length of the unit cell along the nanowire axis. This produces an e ective diameter that scales as the square root of the crosssectional area. In the localdensity functional (LDF) calculations, we evenly sampled 64 discrete points over the central Brillouin zone and used the 321G basis set. [132{134] Since our code takes advantage of the helical symmetry, the calculated band structures look quit di erent from those obtained by using other codes, such as VASP, SIESTA, GAUSSIAN, although they are essentially the same. To be able to compare with results reported by others using translational symmetry, all band structures shown in this chapter have been folded to be within the traditional translational Brillouin zone. Although the LDF is wellknown to underestimate the band gap, we are primarily 85 concerned with the relative change of the energy gaps here. We also carried out some testing calculations with PBE functional on the silicon nanowires of small diameter. We numerically calculated the e ective mass of charge carriers in silicon nanowires using the standard expression from solidstate theory m = ~2=(d2E=dk2), which is equal to ~2=(l2d2E=d 2) if we use helical symmetry, where l is the translation length along wire axis and ~ is the reduced Planck's constant. 86 Figure 4.1: Crosssection views of the studied HSiNWs. h100i SiNWs of diameter (a) 0.43 nm, (b) 0.87 nm, (c) 1.3 nm, and (d) 1.73 nm, respectively. h110i SiNWs of diameter (e) 0.73 nm, (f) 1.09 nm, and (g) 1.46 nm, respectively. The golden yellow and silver balls represent Si and H atoms, respectively. In OHSiNWs and CH3SiNWs, all H atoms are replaced with OH and CH3 groups, respectively. 87 Figure 4.2: Band gaps as a function of diameter with hydrogen passivation for silicon nanowires along h100i. The square and blue line indicate the LDF results and the uptriangle and red line indicate the PBE results. For Hterminated h100i SiNWs, we compared the calculated band gaps with LDF and PBE functionals. Figure 4.2 depicts the band gaps as a function of diameter. The consideration of electron density gradient in PBE functional shifted up the LDF band gaps. The variation in the shifted amount, however, is small with maximum 0.05 eV. The LDF functional is su cient to capture the quantum con nement in this nanomaterial and also have advantage of computional cost with comparison to PBE. The LDF was utilized for the study of silicon nanowires in the following discussions. 88 4.3 Structure Parameters The surface passivating groups can a ect the structure parameters of silicon nanowires. The calculated SiSi bond lengths as a function of the distance from the bond center to the wire axis with di erent passivations are depicted in Figure 4.3. After geometry optimization the Hterminated h100i SiNWs exhibited small varia tions in SiSi bond length from 2.37 A to 2.39 A, compared to the equilibrium value of 2.35 A in the bulk silicon. In the OHpassivated h100i wire of diameter d = 1.73 nm, the SiSi distances have larger variations from 2.36 A to 2.49 A. The SiSi dis tances are ranging from 2.34 A to 2.41 A in the CH3terminated h100i wire of same diameter. The variations of SiSi bond length in the core part are smaller than those near the surface. In the h110i SiNWs, we found that the SiSi bond lengths of HSiNWs are similar to those of h100i HSiNWs. The passivation with with OH or CH3 induces smaller variations in SiSi distances because the distances between surfacepassivating groups on f111g faces of h110i wires are larger than those on f110g faces of h100i wires. The CH3 groups even give rise to bond length variations comparable to that of the Hpassivated wires along h110i direction. 89 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 SiSI (Å) H H OH OH CH3 CH3 <100> 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 <110> 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 SiSi (Å) 0 1 2 3 4 5 6 7 8 9 10 2.2 2.3 2.4 2.5 2.6 0 1 2 3 4 5 6 7 8 9 10 Radius (Å) 2.2 2.3 2.4 2.5 2.6 SiSi (Å) 0 1 2 3 4 5 6 7 8 9 10 Radius (Å) 2.2 2.3 2.4 2.5 2.6 Figure 4.3: SiSi bond length as a function of the distance from bond center to wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm passivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. 90 4.4 Mulliken Population Analysis We also did a Mulliken population analysis to see how the functional groups a ect the charge distribution in the silicon nanowires. The calculated net charge of atoms as a function of the distance from atom to wire axis are presented in Figure 4.4. For elements Si, H, C, and O, the electronegativity (EN) is 1.9, 2.1, 2.5, and 3.5, respectively.[135] So C and O could pull more electrons away from the neighboring Si atoms than H and weaken the bond length of SiSi. Although the EN of C is smaller than that of O, there are three H atoms attached to each C atom such that C atoms have more places to get charges. This results in that the charge is about 0.9 for C and is about 0.55 for O in the h100i wires of diameter d = 1.73 nm. The similar values were observed in the h110i wires of diameter d = 1.46 nm. The net charge of C and O are almost independent of the position and the wire orientation. The independence of position makes the corner Si atoms donate more charges to C or O and leave itself more positively charged because there are two groups attached to the corner Si. We found that the surface Si atoms are positive and the innermost Si atoms are nearly neutral. In the middle part, the net charge of Si atoms exhibited oscillations instead of gradual transition. 91 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 Net charge (e) H OH CH3 H OH CH3 <100> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 <110> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 Net charge (e) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0.5 0 0.5 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Radius (Å) 1 0.5 0 0.5 1 Net charge (e) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Radius (Å) 1 0.5 0 0.5 1 Figure 4.4: Net charge as a function of distance from the atom to the wire axis. (a), (b), (c) correspond to h100i SiNWs of diameter d = 1.73 nm passivated by H, OH, and CH3, respectively. (d), (e), and (f) correspond to h110i SiNWs of diameter d = 1.46 nm passivated by H, OH, and CH3, respectively. Si, H, O, and C are indicated by +, , , and , respectively. 92 4.5 Band Gap Change All of the passivations (H, OH, CH3) studied in the present paper give indirect band gap for the h100i nanowire of diameter d = 0.43 nm. CH3 even gives indirect band gap for all h100i nanowires investigated. Figure 4.5(a) shows the h100i Si nanowire band gaps as a function of diameter. With increasing diameter the band gap is decreasing because of the reduced quantum con nement. Because surface atoms have smaller and smaller percentage when increasing diameter from 0.43 nm to 1.73 nm, surface e ect becomes weak and energy gap di erence of di erent passivations is reduced from 2.75 eV to 0.64 eV for OHSiNWs and from 1.72 eV to 0.54 eV for CH3 SiNWs with HSiNWs as the reference. We observed similar trend in band gap change of h110i SiNWs as shown in Figure 4.5(b). But di erent from the h100i SiNWs, the h110i wires have direct band gaps independent of the surface substituents. The band gap values corresponding to Figures 4.5(a) and 4.5(b) are listed in Tables 4.1 and 4.2. 93 (a) (b) Figure 4.5: Band gaps as a function of diameter with various surface passivations for silicon nanowires along (a) h100i, (b) h110i. 94 Table 4.1: Band gaps (eV) as a function of diameter for h100i silicon nanowires with various passivations. X represents H, OH, or CH3. Stoichiometry Diameter (nm) H OH CH3 Si2X2 0.43 5.61 2.86 3.89 Si4X4 0.87 3.38 1.75 2.42 Si9X6 1.30 2.38 1.48 1.91 Si16X8 1.73 1.89 1.25 1.35 Table 4.2: Band gaps (eV) as a function of diameter for h110i silicon nanowires with various passivations. X represents H, OH, or CH3. Stoichiometry Diameter (nm) H OH CH3 Si4X4 0.73 2.38 1.01 2.07 Si9X6 1.09 1.79 0.37 1.65 Si16X8 1.46 1.49 0.30 1.39 95 4.6 Electronic Structures We show the band structures of h100i silicon nanowires of diameter d = 1.73 nm with di erent surface passivations in Figure 4.6 and in Figure 4.7 we present the band structures of h110i silicon nanowires of diameter d = 1.46 nm. The Fermi levels in HSiNWs are used as the reference and taken to be zero. Both OH and CH3 change the band structure a lot and OH gives rise to more reduction of band gap. In addition to the band gap change, the Fermi levels also vary with the passivations. The CH3 groups shift up the Fermi levels in h100i and h110i wires compared to H SiNWs. The OH groups induce lower Fermi level in h100i SiNWs and higher Fermi level in h110i SiNWs. It is known that the Fermi levels will be aligned via charge redistribution when systems of di erent Fermi levels are put together. If the silicon nanowires are selectively passivated such that one segment is covered with H atoms and neighboring segment is passivated by CH3 groups. Then the electron in CH3 terminated part will move to Hterminated part and hole will move in the inverse way to align the Fermi levels and lower the total energy of the system. The ability of separating charge carriers in selectively passivated SiNWs could nd potential application in photovoltaic devices. Growth techniques for semiconductor nanowires have developed rapidly in recent years. Not only can the diameter and direction be controlled during growth, but nanowires can also be selectively functionalized. This has been achieved by vari ous e orts such as electrochemical methods,[136] localized nanoscale Joule heating, [137] and adsorption and removal of selfassembled monolayers of the polymer (3 mercaptopropyl)trimethoxysilane.[138] 96 Figure 4.6: Electronic band structures of h100i SiNWs (diameter d = 1.73 nm) pas sivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. 97 Figure 4.7: Electronic band structures of h110i SiNWs (diameter d = 1.46 nm) pas sivated by (a) H, (b) OH, and (c) CH3. The Fermi level in HSiNW is shifted to zero eV and taken as the reference. 98 4.7 Neargap States The neargap states are mostly responsible for the transport and optical properties of the system. In order to see how the surface passivations will a ect those proper ties in detail, we plot the orbital density of the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) for the h100i SiNWs of diameter d = 1.73 nm in Figure 4.8 and for the h110i SiNWs of diameter d = 1.46 nm in Figure 4.9 with sub gures (a), (b), and (c) corresponding to the passivating substituents H, OH, and CH3, respectively. First look at the h100i SiNWs. For the HSiNWs, both HOMO and LUMO are concentrated in the interior of the nanowire. Upon OH passivation HOMO goes to the corner O atoms while LUMO leaves the central region and corners blank. In the CH3SiNWs HOMO and LUMO spread out over the whole wire. In the h110i SiNWs, the LUMO orbitals of HSiNWs and CH3SiNWs have similar distribution with more contribution from the central region while the OH groups make the LUMO concentrated in a rectangular region. For the HOMO orbitals of h110i SiNWs, CH3 groups make the orbitals more on the interior atoms compared to the HSiNWs. The OH groups change it dramatically and the HOMO moves to the top and bottom Si and O atoms. And more importantly, the HOMO and LUMO have considerable concentration on O atoms in both h100i and h110i SiNWs, indicating that the OHSiNWs are more sensitive to the external environment than the other two kinds of passivations. 99 Figure 4.8: The HOMO and LUMO orbital density of the h100i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respectively. Here red and blue represent HOMO and LUMO density, respectively. The contour is at 10% of the maximum value. 100 Figure 4.9: The HOMO and LUMO orbital density of the h110i SiNWs. (a), (b), and (c) correspond to SiNWs passivated by H, OH, and CH3, respectively. Here red and blue represent HOMO and LUMO density, respectively. The contour is at 10% of the maximum value. 101 Usually the surface substituents such as hydroxyl, halogens[126, 139] and methyl in h100i SiNWs from our simulations will induce big change in the spatial distribution of neargap orbitals. However, it is interesting to note that the CH3 substituents do not have signi cant e ect on the HOMO and LUMO orbitals in h110i SiNWs. In Figures 4.10 and 4.11 we depict the isosurface of the HOMO and LUMO wave functions in h110i HSiNWs, respectively. We can see that the wave functions have alternating positive and negative values along the wire axis and hence nodal surfaces perpendicular to the wire axis. When applying uniaxial tensile (compressive) strain, the energy values of the HOMO and LUMO will shift downward (upward) because the increase (decrease) of the distance between nodal planes results in the raising (lowering) of kinetic energy. So the conduction band minimum (CBM) and valence band maximum (VBM) will decrease or increase simultaneously. Similar e ect was predicted to be useful in separating electrons and holes when partial strain is applied in h110i SiNWs.[89, 140] Under 2% tension, LUMO is located in strained region while HOMO is in the regular part in a partially strained Silicon nanorod, suggesting photovoltaic application in terms of typeII homojunction solar cells. In our simulations, the CH3 passivation retains the HOMO and LUMO structures of alternating positive and negative values along the wire axis as shown in Figures 4.12 and 4.13. So the h110i CH3SiNWs are supposed to have similar capability of separat ing the positive and negative charge carriers under partial strain as in h110i HSiNWs. In addition to the charge carrier separation, the CH3SiNWs have better air stability and provide higher hole mobility and ono ratio than the HSiNWs [112] and do not a ect the conductance of the neargap channels, o ering quasiballistic transport within several subbands near the SiNWs band gap. [123] 102 Figure 4.10: The HOMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 103 Figure 4.11: The LUMO orbital isosurface of the h110i HSiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 104 Figure 4.12: The HOMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 105 Figure 4.13: The LUMO orbital isosurface of the h110i CH3SiNWs (diameter d = 1.46 nm). Here red and blue represent positive and negative values, respectively. The contour is at 10% of the maximum value. 106 4.8 E ective Mass of Charge Carriers The e ective mass of carriers was calculated from the band dispersions and sum marized in Table 4.3. The hole e ective mass mh of h100i HSiNWs studied is decreas ing with diameter in agreement with results of similar diameters reported by Yan, et al. [114] And in their calculations mh is much bigger than me for h100i HSiNWs. Although heavy hole is also observed in our simulations, the mh is greatly reduced with increased diameter and has comparable magnitude with me especially for the h100i wire of diameter d = 1.73 nm. The OH passivation has e ect on mh and me 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


