ATOMISTIC SIMULATION OF NANOSTRUCTURED MATERIALS

A Dissertation Presented to The Graduate Faculty of the University of Akron

In Partial Fulfillment of the Requirements for the Degree Doctor of Philosophy

Ronghua Zhu December, 2006

ATOMISTIC SIMULATION OF NANOSTRUCTURED MATERIALS

Ronghua Zhu Dissertation

Approved:

Accepted:

______________________________ Advisor Dr. Ernian Pan

______________________________ Department Chair Dr. Wieslaw K. Binienda

______________________________ Committee Member Dr. Wieslaw K. Binienda

______________________________ Dean of the College Dr. George K. Haritos

______________________________ Committee Member Dr. Ala R. Abbas

______________________________ Dean of the Graduate School Dr. George R. Newkome

______________________________ Committee Member Dr. Xiaosheng Gao

______________________________ Date

______________________________ Committee Member Dr. Alper Buldum ______________________________ Committee Member Dr. Peter W. Chung

ii

ABSTRACT

Atomistic based computer modeling and simulation of nanostructured materials has become an important subfield of materials research. Based on the multiresolution method, which combines the continuum mechanics, kinetic Monte Carlo method and molecular dynamics method, we study the nanostructured materials grown by quantum-dot selfassembly, mechanical properties of strained semiconductors, and mechanical properties of carbon nanotube reinforced composites. This thesis covers the following three main contributions. 1). Self-organization of semiconductors InAs/GaAs in Stranski-Krastanov growth mode is studied using kinetic Monte Carlo simulations method coupled with the Green’s function solution for the elastic strain energy distribution. The relevant growth parameters such as growth temperature, surface coverage, flux rate, and growth interruption time are investigated. It is shown clearly that when the long-range strain energy is included in the simulation, ordered uniform size distribution can be achieved. To address the effect of material anisotropy, the anisotropic substrates of GaAs with different growth orientations (001), (111), and (113) and an isotropic substrate Iso (001), reduced from cubic GaAs, are also investigated. Simulation results show that at selected growth

parameters

for

temperature, coverage, and growth interruption time, iii

strain energy field in the substrate is the key factor that controls the pattern of island distribution. Furthermore, layer-by-layer growth of quantum dots is also simulated briefly, and vertical alignment is observed that could lead to progressively uniform island sizes and spatial ordering. 2). Since the misfit strain will be induced during the quantum dots epitaxial growth, the mechanical property of the grown semiconductors will be influenced. In this thesis, utilizing the basic continuum mechanics, we present a molecular dynamic prediction for the elastic stiffness C11, C12 and C44 in strained silicon and InAs as functions of the volumetric (misfit) strain. The bulk modulus, effective elastic stiffness C11, C12 and C44 of the strained silicon, including also the effective Young’s modulus and Poisson’s ratio, are all calculated and presented. In general, our simulation indicates that the bulk moduli, C11 and C12 increase with increasing volumetric strain whilst C44 is almost independent of the volumetric strain. 3). Also using MD method, the carbon nanotube reinforced Epon 862 composite is studied. The stress-strain relations and the elastic Young’s moduli along the longitudinal direction (parallel to CNT) are simulated with the results being also compared with those from the rule-of-mixture. Our results show that, with increasing strain in the longitudinal direction, the Young’s modulus of CNT increases whilst that of the Epon 862 composite or matrix decreases. A long CNT can greatly improve the Young’s modulus of the Epon 862 composite (about 10 times stiffer), which is also consistent with the prediction based on the rule-of-mixture at low strain level. Furthermore, even a short CNT can also enhance the Young’s modulus of the Epon 862 composite, with an increment of 20% being observed as compared to that of the Epon 862 matrix. iv

ACKNOWLEDGEMENTS

This thesis is the end of my long journey in obtaining my doctor’s degree in Civil Engineering. I have not traveled in a vacuum in this journey. There are some people who made this journey easier with words of encouragement, with hands of help, with experience of guidance …… To my advisor Dr. Ernian Pan for his help, encouragement and suggestions during the past three and half years! Deeply indebted to Dr. Ernian Pan’s kindness and friendship. His diligence and unremittingness set me an excellent example, which will be helpful to me during my whole life. To my committee members: Dr. Wieslaw K. Binienda, Dr. Ala R. Abbas, Dr. Yu Qiao, Dr. Xiaosheng Gao, Dr. Alper Buldum and Dr. Peter W. Chung for their time, kind support and helpful suggestions. To my colleagues, Feng Han, Wael Alkasawn, Yan Zhang, Xu Wang, Mingkun Sun, Yuanguo Chen, James Ramsey and Kotchanu Jittava, for giving me the feeling of being at home at work. To my parents, my family and friends …… And, to Army Research Office for their support.

v

TABLE OF CONTENTS Page LIST OF FIGURES …….….……………………..…………………………………….. x CHAPTER I INTRODUCTION…….….……………………..……………………………………..... 1 1.1

Atomistic Simulation in Nanostructured Materials........................................... 1

1.2

Multiresolution in Atomistic Simulation........................................................... 2

1.3

Research Purpose .............................................................................................. 4

II MONTE CARLO SIMULATION OF SELF-ORGANIZED QUANTUM DOT GROWTH…………………………………………………………………………

6

2.1

Quantum Dots.................................................................................................... 6

2.2

Self-organized Growth ...................................................................................... 8

2.3

2.2.1

Deposition and Diffusion........................................................................ 9

2.2.2

Growth Classification ........................................................................... 10

Kinetic Monte Carlo Simulation Theory......................................................... 13 2.3.1

Nucleation and Growth of Quantum Dots ............................................ 13

2.3.2

Atomic Hopping Probability................................................................. 16

2.3.3

Kinetic Monte Carlo Algorithm............................................................ 19

III GROWTH PARAMETERS DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION………………………………………………... 23 3.1

Growth Parameters .......................................................................................... 23

3.2

Temperature T ................................................................................................. 24 vi

3.3

Flux Rate F...................................................................................................... 24

3.4

Surface Coverage c.......................................................................................... 25

3.5

Growth Interruption Time ti ............................................................................ 26

IV ELASTIC ANISOTROPY DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION……………………………………………….. 30 4.1

Strain Field in Quantum Dot Growth .............................................................. 30

4.2

Anisotropic Strain Energy Field Calculation .................................................. 32 4.2.1 Material Properties of Iso(001), GaAs(001), GaAs(111) and GaAs(113)......................................................................................................... 32 4.2.2 Strain Energy Field of Iso(001), GaAs(001), GaAs(111) and GaAs(113)......................................................................................................... 34

4.3

Effects of Elastic Strain Energy on Island Ordering ....................................... 35

4.4

Effect of Substrate Anisotropy and Growth Parameters on Island Ordering .. 37

4.5

Stacks of Quantum Dots.................................................................................. 41

V KINETIC MONTE CARLO SIMULATION OF THREE DIMENSIONAL QUANTUM DOT SELF-ASSEMBLY ……………………………………………54 5.1

Overview of Three Dimensional Quantum Dots Growth ............................... 54

5.2

Three Dimensional Model Description ........................................................... 56

5.3

Simulation Results........................................................................................... 59

5.4

Conclusion of 3D Simulation Model............................................................... 62

VI MOLECULAR DYNAMICS CALCULATION OF ELASTIC MODULI IN STRAINED SEMICONDUCOTR………………………………………………….70 6.1

Background of Elastic Moduli in Strained Semiconductor............................. 70

6.2

Elastic Stiffness Calculation............................................................................ 72 6.2.1

Bulk Modulus........................................................................................ 74

6.2.2

Modulus C11-C12 ................................................................................... 75

vii

6.2.3 Elastic Stiffness C44 .............................................................................. 76 6.3

6.4

6.5

Molecular Dynamics of Elastic Moduli Calculation ....................................... 77 6.3.1

Classical Interatomic Potential ............................................................. 78

6.3.2

MD Programming ................................................................................. 79

Elastic Stiffness Calculation of InAs............................................................... 80 6.4.1

Bulk Modulus........................................................................................ 81

6.4.2

C11 and C12 ............................................................................................ 81

6.4.3

C44 ......................................................................................................... 81

Elastic Stiffness Calculation of Silicon ........................................................... 82 6.5.1

Potential and Volumetric Strain............................................................ 83

6.5.2

C11, C12 and C44 ..................................................................................... 84

6.6

Effective Young’s modulus and Poisson’s ratio ............................................. 84

6.7

Conclusion of MD Simulation of Strained Semiconductor............................. 86

VII MOLECULAR DYNAMICS STUDY OF THE STRESS-STRAIN BEHAVIOR OF CARBON-NANOTUBE REINFORCED EPON862 ……………………………..96 7.1

Overview of Carbon Nanotube Composites.................................................... 96

7.2

Molecular Dynamics Simulation of CNT Composites.................................... 97 7.2.1

Model Building ……………………………………………………….98

7.2.2 Simulation of Stress-Strain Relationship ……………………………100 7.2.3 Rule-of-Mixture ……………………………………………………..102 7.3

Numerical Results ......................................................................................... 103

7.4

Concluding Remarks of MD Simulated CNT Reinforced Epon 862 ............ 105

VIII CONCLUSIONS ………………………………………………………………….112 REFERENCES ............................................................................................................. 114 APPENDICES ………………………………………………………………………...122 viii

APPENDIX A HALF-SPACE GREEN’S FUNCTIONS ...................................... 125 APPENDIX B DETAILED SIMULATION STEPS.............................................. 129 APPENDIX C VOLUME-CONSERVING OF ORTHORHOMBIC STRAIN..... 132 APPENDIX DVOLUME-CONSERVING OF MONOCLINIC STRAIN............. 133

ix

LIST OF FIGURES Figure

Page

1.1 Different laws of physics are required to describe properties and processes of fluids at different scales [E03]. ..................................................................................... 5 2.1 Possible ways of diffusion from A to B: a - desorption, condensation process; b surface diffusion; c - volume diffusion....................................................................... 21 2.2 Important growth modes in epitaxy: a) Volmer-Weber growth mode, b) FrankVan der Merwe growth mode, c) Stranski-Krastanov growth mode. ......................... 21 2.3 Strain relief by a) generation of mis-fit dislocations (green atoms) and b) StranskiKrastanov growth mode. The green line marks the end of the wetting layer. Atoms of the substrate are blue, deposited material is red. ........................................ 22 2.4 Schematic illustration of off, edge, and corner diffusions. ......................................... 22 2.5 Schematic illustration of nearest (horizontal and vertical directions in solid blue) and next nearest neighbor (diagonal directions in dashed brown) positions of each atom............................................................................................................................. 22 3.1 Atom island ordering under different temperatures (T is equal to 550K in (a), 600K in (b), 650K in (c), 700K in (d), 750K in (e), and 800K in (f)). Simulation parameters are F=1.0Ml/s, c=20% and ti=200s on a 200×200 grid. The strain field is updated every 2500 steps. ....................................................................................... 27 3.2 Island ordering of atoms under different flux rates. Simulation parameters are T=750K, c=20% and ti=200s on a 200×200 grid. Strain field is NOT included. Deposition stops after 0.2s in (a), 2s in (b), and 20s in (c). ........................................ 28 3.3 Atom island ordering with coverage c increasing from 10% to 50% in (a) to (e). Simulation parameters are T=700K, F=1.0Ml/s, ti=200s and on a 200×200 grid. The strain energy field is updated every 2500 steps. .................................................. 28

x

3.4 Atom island ordering for different interruption times. The ordering of (a) is the initial atom pattern when all the atoms just deposit to the surface and the interruption time ti=0. Simulation parameters are T=750K, F=1.0Ml/s and c=20% on a 200×200 grid. The strain energy is updated every 2500 steps............................ 29 4.1 Schematic illustration of crystal coordinates .............................................................. 44 4.2 Contours of normalized elastic strain energy distribution on the substrate surface of isotropic Iso(001) (a), and anisotropic GaAs(001) (b), GaAs(111) (c), and GaAs (113) (d). ........................................................................................................... 44 4.3 Island ordering on the surface of GaAs without elastic strain energy (a), with strain energy of Iso (001) (b), with anisotropic strain energy of GaAs (001) (c), GaAs (111) (d), and of GaAs (113) (e). Different island chain orderings (red dashed lines) corresponding to different growth orientations can be observed. Other simulation parameters are T=750K, F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid............................................................................................................... 45 4.4 Island density contours. The whole 200×200 grids are divided into 18×18 big grids, so there is 11×11 original grids in one big grid. The densities are plotted by 3 dimensional contours. .............................................................................................. 46 4.5 Surface island distributions vs. temperatures (550K, 650K, 750K) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations. ..................... 47 4.6 Number of islands and average island size vs. temperature for islands growth on substrates with different orientations for short temperature range (550K to 750K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s................... 48 4.7 Number of islands and average island size vs. temperature for islands growth on substrate GaAs (001) for long temperature range (550K to 950K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s. ............................ 48 4.8 Surface island distributions vs. coverage c (10%, 20% 30% and 50%) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations. ..................... 49 4.9 Number of islands and average island size vs. coverage for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and ti=200s................................................................................ 50 xi

4.10 Surface island distributions vs. interruption times ti (0s, 10s, and 200s) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and c=20%, on a 200×200 grid. Red dashed lines indicate the island chain orientations.................................................................................................................. 51 4.11 Number of islands and average island size vs. interruption time for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and c=20%...................................................................... 52 4.12 Schematic illustration of superlattice........................................................................ 52 4.13 Atom island ordering on the first three layers. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Thickness of each layer is 5 grids. Strain energy field is calculated at every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for growth interruption. ........................................................................................ 53 4.14 Atom island ordering for different layer thickness. Figures (b) to (h) are the simulated results on the second layer with thickness equal to 5 grids to 30 grids. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Strain energy field is updated every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for interruption. ...... 53 5.1 Schematic illustration of 3D QD self-assembly model. An atom on top of the substrate surface (x,y)-plane in (a), the relative locations of the atom grid on the (x,y)-plane in (b), the corresponding binding energy Es related to the in-plane locations in (c)............................................................................................................. 64 5.2 Illustration of an edge atom A’s jump up or down during 3D QD self-assembly. ..... 64 5.3 3D island distributions for different up/down ratios ρ with total coverage c=1.6Ml. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and interruption time ti=200s (The total simulation time = the deposition time (td=160s) plus the interruption time (ti=200s)). .......................................................................... 65 5.4 Average island height vs. up/down ratio ρ (for ρ from 3 to 50). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and interruption time ti=200s. ..................................................................... 65

xii

5.5 Island distributions for different interruption times (ti=0s, corresponding to deposition time td=160s, ti=100s, and ti=200s) and for different up/down ratios (ρ=10, 20 and 30). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and total coverage c=1.6Ml. .................................................................. 66 5.6 Island distributions during and after deposition. Deposition time td=50s with 0.5ML coverage in (a), td=100s with 1ML coverage in (b), td=160s with 1.6Ml coverage (i.e., at the beginning of interruption time) in (c) and interruption time ti=200s in (d). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and up/down ratio ρ=10.................................. 67 5.7 Average island diameter (of the bottom equivalent circle of the island) vs. number of islands during and after deposition. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml and up/down ratio ρ=10..... 68 5.8 Island distributions for different deposition processes of 0.5Ml, 1Ml and1.6 Ml with flux rates F=1Ml/s, 0.1Ml/s and 0.01Ml/s. Fixed growth parameters are temperature T=700K, total coverage c=1.6Ml and up/down ratio ρ=10. ................... 69 6.1 Illustration of the geometry of a triplet of atoms used in the definition of the threebody potentials. ........................................................................................................... 86 6.2 Flowchart of molecular dynamics programming process........................................... 87 6.3 Illustration of cubic crystalline under uniform strain. ................................................ 88 6.4 Effective bulk modulus of InAs is plotted as a function of volumetric strain. Solid line is the results obtained by our program. Dashed line with triangles is the results from the ref. [Ell02]. Dashdotted line is calculated by commercial software. ..................................................................................................................................... 88 6.5 Effective elastic constant C11-C12 of InAs is plotted as a function of volumetric strain. Solid line is the results obtained by our program and dashed line is from the ref. [Ell02]............................................................................................................. 89 6.6 Effective elastic constants C11 and C12 of InAs are plotted as a function of volumetric strain. Solid line is the results obtained by our program. Dashed line with triangles is from the ref. [Ell02]. Dashdotted line is calculated by commercial software....................................................................................................................... 90 6.7 Effective elastic constant C44 InAs is plotted as a function of volumetric strain. Solid line is the results obtained by our program. Dashed line is from the ref. [Ell02]. Dashdotted line is calculated by commercial software. ............................... 91 xiii

6.8 Lattice energy (per atom) vs single silicon cubic lattice length.................................. 91 6.9 Lattice energy (216 atom system) vs system volume (V-2/3). ..................................... 92 6.10 Bulk modulus of Silicon vs volumetric strain........................................................... 92 6.11 C11-C12 of Silicon vs volumetric strain. .................................................................... 93 6.12 C11 of silicon vs volumetric strain. ........................................................................... 93 6.13 C12 of silicon vs volumetric strain. ........................................................................... 94 6.14 C44 of Silicon vs volumetric strain............................................................................ 94 6.15 Young’s modulus plotted as a function of volumetric strain.................................... 95 6.16 Poisson’s ratio plotted as a function of volumetric strain......................................... 95 7.1 Schematic arrangement of Epon 862 matrix filled with long (a) and short (b) CNTs. ........................................................................................................................ 106 7.2 Epon 862 molecular structure (a) and computer constructed molecules of Epon 862 (b)....................................................................................................................... 107 7.3 Computer constructed Epon 862 matrix with the size of l1×l2×l3= 40.28×40.28×61.09 Å3. (The CNT is embedded along the longitudinal x3direction as shown in Figure4).................................................................................. 107 7.4 Long (a) and short (b) CNT (10,10)-reinforced Epon 862 matrix (40.28×40.28×61.09 Å3). Lengths of the long and short CNTs are 61.09Å and 29.32Å, respectively. ................................................................................................ 108 7.5 Equilibrium van der Waals separation distance between the CNT (10,10) and Epon 862 matrix........................................................................................................ 108 7.6 MD simulated stress-strain curves for CNT (10,10), long and short CNT composites, and Epon 862 matrix............................................................................. 109 7.7 MD simulated stress-strain curves for CNT-reinforced Epon 862 composites and Epon 862 matrix........................................................................................................ 109 7.8 Stress-strain curves of the long CNT composite: MD simulated vs. Rule-ofmixture. ..................................................................................................................... 110 xiv

7.9 Variation of MD simulated Young’s moduli of CNT (10,10) and Epon 862 matrix with increasing strain. ............................................................................................... 110 7.10 Variation of MD simulated Young’s moduli of long and short CNT-reinforced Epon 862 composites and Epon 862 matrix with increasing strain. ......................... 111 B1 Extended periodic condition...................................................................................... 130 B2 Atom neighbors. ........................................................................................................ 130 B3 Flow chart 2D KMC QDs Growth Model ………………………………………….130

xv

CHAPTER I

INTRODUCTION

1.1 Atomistic Simulation in Nanostructured Materials

Nanostructured materials have greatly attracted people’s attention in the past ten years due to two main reasons. First, they are attractive from a scientific point of view. Nanostructured semiconductors provide a means to create artificial potentials for carriers, electrons, and holes at length scales comparable to or smaller than the de Broglie wavelength. And, nanostructured composites promise to dramatically increase material strength and durability [Sie96] as well. Second, quantum mechanics becomes applicable not only in systems of academic interest, but also to systems of practical impact. In particular, semiconductor nanostructures have a large potential for applications in nanoand optoelectronics. Over the past few years, computational materials modeling has become an increasingly important subfield of materials research, with its growing importance a result of the enormous increase in computing power and the maturation of a number of computational methodologies. Modeling has, in many cases, moved from being a tool for finding explanations of materials properties to being a means for prediction of the 1

properties of new classes of materials. It is a field that encompasses a large variety of techniques applied over a correspondingly disparate range of length and time scales, describing the structure and properties of wide classes of materials systems. The choice of a technique is, of course, dependent on the system under consideration and the phenomena that are being modeled. At feature sizes smaller than 100 nm, conventional theoretical modeling and computer simulations, such as those based on continuum approaches, must be supplemented with atomistic models. With conventional approaches, such as the finiteelement method based on linear elasticity, modelers cannot include the effects of defects, surfaces, and other microstructures – all of which fundamentally affect the material’s performance. Atomistic simulations, in particular those based on molecular dynamics (MD) and Monte Carlo (MC), capture these microstructural effects as well as highly nonhomogeneous stresses at the atomic scale. Simulations based on state-of-the-art molecular dynamics can involve up to 100 million atoms, but many researchers believe that petaflop computers will be built by the end of the next decade (one petaflop is 1015 floating-point operations per second). Such computers will make it possible to simulate nanostructures that involve 1012 atoms – with a total system size of well over one micron. This in turn means that modelers can simulate an entire device structure atomistically.

1.2 Multiresolution in Atomistic Simulation

Multisolution modeling and computation is a rapidly evolving area of research that has a fundamental impact on computational science and applied mathematics and 2

influences the way we view the relation between mathematics and science. Currently, multiresolution problems are driven mainly by the use of mathematical models in the applied sciences: in particular, material science, chemistry, fluid dynamics, and biology. Though, modeling at the level of a single scale, such as molecular dynamics or continuum theory, is becoming relatively mature. There is an urgent need from science and technology - nano-science being a good example - for multiresolution modeling techniques. It is not an exaggeration to say that almost all problems have multiple scales. We organize our time according to days, months, and years, reflecting the multiple time scales in the dynamics of the solar system. In addition, different physical laws may be required to describe the system at different scales. Take the example of fluids (figure 1.1) [E03]. At the macroscale (meters or millimeters), fluids are accurately described by the density, velocity, and temperature fields, which obey the continuum Navier-Stokes equations. On the scale of the mean free path, it is necessary to use kinetic theory (Boltzmann’s equation) to get a more detailed description in terms of the one-particle phase-space distribution function. At the nanometer scale, molecular dynamics in the form of Newton’s law has to be used to give the actual position and velocity of each individual atom that makes up the fluid. A general framework for multiresolution methods should ideally unify existing methods, give guidelines on how to design new methods, improve existing ones and provide a mathematical theory for stability and accuracy of these methods [E03]. There is a long history in mathematics for the study of multiresolution problems. Fourier analysis has long been used as a way of representing functions according to their components at 3

different scales [Bak99]. On the computational side, several important classes of numerical methods have been developed which address explicitly the multiresolution nature of the solutions. These include multigrid methods [Wes92], domain decomposition methods [Tos05], fast multipole methods [Gum05], and adaptive mesh refinement techniques [Vay04]. Many other mathematical techniques have been developed to study multiresolution problems, including boundary-layer analysis [Kev96], semiclassical methods [Mas81], geometric theory of diffractions [Kel62], stochastic mode elimination [Maj00], and renormalization group methods [Cho03].

1.3 Research Purpose

In our research, combining continuum theory, we focus on using MC and MD to study quantum dots (QDs) self-assembly and carbon nanotube (CNT) composites. The growth of self-organized QDs in strained semiconductors in the Stranski-Krastanov growth mode has been studied using kinetic Monte Carlo (KMC) simulations method. The set of relevant growth parameters such as growth temperature, surface coverage, flux rate, and growth interruption time are investigated. The dependence of the self-organized quantum dots pattern on elastic strain energy of the substrate is studied as well. In our research, by incorporating the anisotropic strain energy field into a KMC algorithm for adatom diffusion, we show that the self-organized island pattern on the surface of an anisotropic substrate is closely correlated to the elastic energy distribution on the surface. Because of the misfit strain during the QD growth, we present a molecular dynamic prediction for the elastic stiffnesses C11, C12 and C44 in strained silicon and InAs as 4

functions of the volumetric strain level. Our approach combines basic continuum mechanics with the classical molecular dynamic approach, supplemented with the Stillinger–Weber potential. Using our approach, the bulk modulus, effective elastic stiffnesses C11, C12 and C44 of the strained silicon, including also the effective Young’s modulus and Poisson’s ratio, are all calculated and presented in terms of figures and formulae. Except for QD self-assembly, using MD method, we study CNT composites as well. Single-walled carbon nanotubes (CNTs) are used to reinforce epoxy Epon 862 matrix. Three periodic systems --- a long CNT reinforced Epon 862 composite, a short CNT reinforced Epon 862 composite, and the Epon 862 matrix itself ---, are studied using molecular dynamics. The stress-strain relations and the elastic Young’s moduli along the longitudinal direction (parallel to the CNT axis) are simulated and the results are compared with those from the rule-of-mixture.

Figure 1.1 Different laws of physics are required to describe properties and processes of fluids at different scales [E03].

5

CHAPTER II

MONTE CARLO SIMULATION OF SELF-ORGANIZED QUANTUM DOT GROWTH

2.1 Quantum Dots

Quantum dots (QD), coherent inclusions in a semiconductor matrix with truly zerodimensional electronic properties, present the utmost challenge and point of culmination of semiconductor physics. It was at the beginning of the 1990s that a modified StranskiKrastanow growth mechanism driven by self-organization phenomena at the surface of strongly strained heterostructures was realized for the fabrication of such dots. And this process presents a sound way to fabricate easily and fast large densities of quantum dots [Bim99]. A semiconductor quantum dot is the ultimate quantum confined structure. Its unique electronic properties rely on the delta-function-like energy dependence of the density of states due to the quantum confinement of carriers in all three dimensions. However, to exploit the electronic properties in novel quantum effect devices, the lateral dimensions of the structures have to be in the range of, or smaller than, the de Broglie wavelength of electrons (50 nm in GaAs) [Bim99]. 6

Moreover,

for

practicable

applications in semiconductor devices, millions of the quantum structures, densely packed and uniform on the atomic scale, are necessary to achieve the desired active volume. This requires more precise fabrication methods to provide better control of size and shape of large ensembles of nanostructures. If this can be done successfully, it will allow the realization of novel high-speed, quantum interference, optoelectronic and single-electron devices. Ultralow threshold currents with high temperature stability have been predicted for quantum dot lasers [Ara82] and high mobilities at room temperature in quantum-dot superlattices with appropriate choice of electronic coupling to suppress scattering by optical phonons. Advanced crystal growth techniques such as molecular beam epitaxy (MBE) and metalorganic vapour phase epitaxy (MOVPE) have made it possible to precisely fabricate two-dimensional layered semiconductors, i.e. heterojunctions, quantum wells and superlattices, on an atomic scale. Further reduction of the dimensionality of semiconductors in one-dimensional quantum wires and zerodimensional quantum dots, however, has turned out to be very difficult. In the past decade, the fabrication methods for nanometer-scale structures have been continuously improved and many fundamental aspects of low-dimensional semiconductors have been evaluated

[Kas90].

However,

the

reproducible

fabrication

of

device-quality

nanostructures is still a challenge, making the applicability of quantum wires and quantum dots somewhat haphazard. Originally, the fabrication methods of quantum wires and quantum dots were based on the lateral patterning of two-dimensional heterostructures by combining fineline lithography with wet and dry chemical etching. 7

Unfortunately,

however,

lithographic patterning always introduces irregularities in the shape of the nanostructures, and mechanical damage to the interfaces cannot be avoided. Therefore, the most promising quantum structures so far have been fabricated using techniques based on direct crystal growth. Along with the selective growth on prepatterned and masked substrates, the direct synthesis of nanostructures during the epitaxial growth process itself has become very important due to its potential for creating damage-free structures [Not93].

2.2 Self-organized Growth

Self-organized growth is a promising way to grow QD islands with uniform size and spatial distribution. However, it is still very challenging to self-organizationally grow high density and dislocation free quantum dot layers, even under the Stranski-Krastanov growth model [Dar97, Shc99b, Spr00]. Since experimentally this trial-and-error growth approach could be quite expensive and time-consuming, numerical simulation methods can be utilized to provide important information and guidance to the experimentation. As the self-organized QD growth is a competitive balance process between the binding and thermal energies, involving atom deposition and diffusion, certain kinetic Monte Carlo (KMC) models have been proposed to numerically simulate the self-organized growth [Mei03]. It was observed by Schöll and Bose [Sch98] and Elsholz et al. [Els03a] that a probability governed by the Arrhenius law including the sole atomic binding energy would result in the Ostwald ripening phenomenon in the KMC model. However, when the long-range strain energy is considered, cooperative growth where larger islands 8

lose some atoms to smaller ones could be achieved [Spr00].

2.2.1 Deposition and Diffusion The first step towards the formation of quantum dots is the deposition of a certain material on a given surface. Molecular beam epitaxy (MBE) belongs to the physical vapor deposition (PVD) methods. The substance to be deposited is vaporized, and the molecules formed into a beam. This beam is under ultra high vacuum conditions directed towards a surface, where the particles condensate. To allow for diffusive motion of the deposited atoms, the sample is heated. Since no other chemically active substances are used in this form of epitaxy, deposition and diffusion are characterized by only the most basic processes. For this reason the growth technique of MBE is very convenient for theoretical modeling of growth. On the other hand in chemical vapor deposition (CVD) the growth material is brought to the sample in form of a chemical carrier gas solution. A prominent example is metal organic CVD (MOCVD). Deposition occurs via chemical reactions at the sample surface. These reactions can be very complex and since the deposition rate depends sensitively on the concentrations of the various chemically active species the local growth kinetics can become quite involved. Furthermore, reactions usually occur both ways and consequently desorption is much more important in CVD than in MBE. Another widely used epitaxial method is the liquid phase epitaxy (LPE), where the substrate is submerged into an oversaturated liquid solution. Material precipitates from the solution and is deposited on the sample. Here, again, the relevant deposition processes are basic but diffusion is widely aided by the presence of a liquid phase. 9

While in MBE and low pressure CVD the interaction of the growing surface with the ambient atmosphere is of minor importance, the theory of growth from dense phases like LPE or growth from the melt can become most complex. Here, for instance, a simultaneous solution of the Navier-Stokes equation and mass- and energy-transport equations is required. Obviously, the technical details in various epitaxial setups can differ strongly. Nevertheless all methods used for epitaxy follow a simple scheme that consists of transporting material to the sample surface, depositing it and allowing for diffusion. The deposition might be as simple as scattering atoms on the surface (as in MBE or related sputtering techniques) or complex and dominated by chemical reactions (as in CVD) but ultimately any of the above mentioned techniques ensures a certain flux of particles to the surface. The flux is the relevant quantity in terms of growth kinetics and its physical effect is basically independent of the applied method of deposition. Once an atom is deposited on the surface it can travel from A to B on various paths (figure 2.1). Plain surface diffusion (figure 2.1b) is an important mechanism in all growth techniques and consists of consecutive hops from one lattice site to a neighboring one. Since it is the energetically most favorable way of diffusion, at least for short distances, it generally dominates all other ways of transportation.

2.2.2 Growth Classification If material A is deposited on material B it is not at all clear in which way the growth will occur. Additionally, for a given material system the mode of growth depends on external parameters like temperature and pressure. The simplest growth mode of 10

hetereogrowth, where one complete monolayer grows after the other is rather the exception than the rule. Under certain conditions rare gases grow layer-by-layer on graphite. Another example is the AlxGa1-x film growth on sapphire (001) as reported in [Wic94]. This growth mode is also called Frank-Van der Merwe [Fra49] growth mode (Figure 2.2b). Conversely, in growth experiments with lead on graphite the lead does not form monolayers. It rather forms little droplets very similar to water droplets on a freshly sealed car top. This growth mode is referred to as the Volmer-Weber [Vol26] growth mode (figure 2.2a). Indeed, this analogy between lead and water droplets is worth pursuing. By defining the surface tensions for a liquid droplet on a plane surface as σsv, σsl and σlv for the interfaces solid/vapor, solid/liquid and liquid/vapor, respectively, one can show that the droplet forms an angle θ with the solid given as:

cosθ =

σ sv − σ sl σ lv

(2.1)

For σsv Ecorner>Eedge.

2.3.2.2 Incorporation of the strain energy field As we have pointed out in the introduction, while the strain energy field is very important in the self-organized growth, it needs to be accurately calculated. In other words, for the strained semiconductor structure, one should use the strain energy corresponding to the point-eigenstrain, instead of that to the point-/line-force. We briefly present the correct formulation for the calculation of the strain energy below, leaving the details in Appendix A for the sake of reference. Within the framework of continuum elasticity theory, the elastic strain field induced by the islands can be obtained using the Green’s function solution for the anisotropic semiconductor substrate developed by refs [Pan02a, Pan02b]. Assume that there is a point misfit strain (or eigenstrain) γ ij* (due to the misfit lattice difference) at x, the induced elastic strain at y can be found as

1 2

k p * γ kp ( y ) = γ lm [σ ml , p ( x; y ) + σ ml , k ( x; y )] y

y

(2.6)

where the superscript k or p (=1,2,3) on the right side of equation (2.6) indicates the direction of the point-force, and the subscript prime followed by py or ky in the stress component σml denotes the derivative of the Green’s stress components. With the elastic strain equation (2.6), the strain energy at y due to an area A of island (with unit thickness in the out-of-plane z-direction) centered at x can be calculated as E str =

1 C ijkl γ ij γ kl A 2

(2.7)

18

where Cijkl is the elastic moduli. In the following calculation, we assume that the misfitstrain is hydrostatic, i.e., γ Ij* = γ *δ ij with γ*=0.07. The non-zero elastic coefficients for GaAs (001) [Pan02b] are taken to be C11=118, C12=54, and C44=59 (109N/m2).

2.3.3 Kinetic Monte Carlo Algorithm The simulation routine is based on a continuous time Monte Carlo scheme. The BKL algorithm [Bor75] named after Bortz, Kalos, and Lebowitz is very efficient for the problem at hand, since the independence of a particular time scale is of great advantage in simulating surface diffusion. In our program, atom diffusion processes are simulated one by one. The principal course of a simulated diffusion process is always the same and is sketched in the following. Each atom at most has four possible nearest and four possible next-nearest diffusion positions (figure 2.5), or four nearest and four next-nearest neighbors. Every possible diffusion step of a given atom is done by evaluating the probability p from equation (2.2). These eight possible diffusion probabilities to nearest and next nearest neighbor positions are stored and added to get the total diffusion probability for the atom, named pAtom. The total probability that anything might happen at all is then given by adding up all the pAtom to get pTotal. During the simulation, an atom is randomly chosen to move across the surface by hopping to nearest or next nearest lattice site. The key is to assign the proper weight to each atom to be chosen, considering that different atom contributes differently to the overall probability pTotal (see, e.g. equations. (2.2) & (2.5)). From the eight possible diffusion processes of the selected atom, one is chosen in accordance to its likelihood and executed. The corresponding time interval ∆t is 19

calculated and added to the elapsed simulation time. The movement of an atom also alters the diffusion barriers for the neighboring atoms in the previous neighborhood as well as in the new one. As such, the moving atom and all the atoms in its surrounding have to be recalculated to obtain the new diffusion probabilities. Strictly speaking, the strain energy field would have to be calculated at each step. However, since the strain evaluation is a lengthy procedure and the strain only slightly changes with the motion of a single atom, it turns out that it is sufficient to calculate the strain field, say, every 2500 jump steps. To speed up further the computation, the strain energy calculation does not extend over the whole system but only over a circular area of a given radius r around the point where the strain is to be evaluated, due to the rapid decay of this field. In the following examples, such a radius r is taken to be 30 lattice grids [Sch98]. The detailed simulation steps are listed in Appendix B.

20

B

a b A c

Figure 2.1 Possible ways of diffusion from A to B: a - desorption, condensation process; b -surface diffusion; c - volume diffusion

a)

b)

c)

Figure 2.2 Important growth modes in epitaxy: a) Volmer-Weber growth mode, b) FrankVan der Merwe growth mode, c) Stranski-Krastanov growth mode.

21

a)

b)

Figure 2.3 Strain relief by a) generation of mis-fit dislocations (green atoms) and b) Stranski- Krastanov growth mode. The green line marks the end of the wetting layer. Atoms of the substrate are blue, deposited material is red.

Eoff Eedge Ecorner

y

x Figure 2.4 Schematic illustration of off, edge, and corner diffusions.

y

x Figure 2.5 Schematic illustration of nearest (horizontal and vertical directions in solid blue) and next nearest neighbor (diagonal directions in dashed brown) positions of each atom. 22

CHAPTER III

GROWTH PARAMETERS DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION

3.1 Growth Parameters

For most electronic and optical applications that include active quantum dot layers, it is important to have as many and as equally shaped dots as possible. Unfortunately, the problem of generating self-organized QD patterns cannot be solved in a straightforward manner. Even though Stranski-Krastanov growth is known to be suitable for selforganized QD growth, there are many parameters that can considerably influence the growth result. Some of the relevant parameters are temperature T, flux rate F, coverage c, and interruption time ti. In what follows, we discuss the effect of these parameters on the island growth of material GaAs (001) on a 200×200 lattice size with periodic condition and with strain energy for all atoms being updated every 2500 hopping steps.

23

3.2 Temperature T

The growth temperature T affects the simulation result greatly. If T is chosen too low, the deposited atoms will just stick to the surface without having enough thermal energy to diffuse. In this case, of course, no self-organization is to be expected. If, on the other hand, the temperature is too high, then interatomic bonds are too easily overcome and one observes an ensemble of monomers and small polymers of atoms performing random walks over the surface. Under this situation, large islands are inherently unstable. Therefore, only for a distinct interval of temperatures self-organization is effective and accessible to production processes. The temperature dependence of the average island size is shown in figure 3.1 for fixed values of flux F=1Ml/s and coverage c=20%. Here, deposition stops after 0.2s of simulation time and the growth interruption time ti is 200s. It is observed from figure 3.1 that the higher the temperature, the larger the average island size, and that when the temperature is higher than 700K, e.g. from 700K to 800K, there is no big difference among the ultimate island size; however increasing temperature in this interval, we obtain more and more ordered island ordering.

3.3 Flux Rate F

The flux rate F is a measure of how fast the material is deposited on the surface. It does not include any equation concerning processes connected to self-organization, like the definition of diffusion barriers or the strain energy field. Thus, the flux rate 24

influences the kinetics of growth only indirectly by determining the density of atoms on the surface and the nucleation rate of islands. Since flux is only present during the time of deposition, one might assume that it is even less important during growth interruption. However, different values of flux can lead to significant differences in the surface morphology after the end of deposition. Under low flux rate, the system will have enough time during the deposition process to come close to an equilibrium size distribution with large-size islands. On the other hand, a small-size island distribution is associated with a high flux rate. In other words, in general, the lower the flux rate, the larger the equilibrium island size, as is also illustrated in figure 3.2 for different flux rates but with fixed temperature T=700K, coverage c=20%, and interruption time ti=200s. We point out that to highlight the effect of flux rate, strain energy is not included in this simulation.

3.4 Surface Coverage c

The surface coverage c describes how much material is deposited on the surface. This parameter certainly has an effect on size distributions, because one cannot expect islands of equilibrium size if there is not enough deposited material. On the other hand, if too much material is transported to the surface than necessary for optimal island sizes, too small distance between islands makes the formation of isolated islands impossible. For coverage below certain critical coverage cc, the growth is mainly kinetically controlled, so the atom size distribution is mainly controlled by a variation of coverage and one should expect an increase of the average island size with increasing coverage. Under this condition, every island collects atoms from its 25

immediate

vicinity

without

any

considerable exchange of material among islands, and its size is directly determined by the amount of the deposited material around the island. For temperature T=700K, flux F=1Ml/s, and interruption time ti=100s, the simulated results on the 200×200 grid with different coverage are plotted in figure 3.3. From the results, we can find that with coverage increase from 10% to 30%, the average island size becomes larger and larger. However, when the coverage is larger than 30%, isolated islands cannot be easily observed. In other words, the critical coverage is between 20% and 30%.

3.5 Growth Interruption Time ti

Another important factor to influence the growth result is the time between the end of deposition and the capping of the QD layer with another material, the so-called growth interruption. The time of this scale is called interruption time ti, with the time from the very beginning to the end of simulation is named total simulation time t (= ti + deposit time). It is known that the growth interruption affects the crystal surface dramatically and the results before and after interruption are different [Kam96]. Furthermore, during the growth interruption, the atoms can move to energetically favorable positions and approach thermodynamic equilibrium. Since there is a striking difference between kinetically controlled growth and a thermodynamically dominated size distribution, the effect of a growth interruption can be dramatic. With the increase of interruption time, the system exhibits ordered island size and ordering. For temperature T=750K, flux F=1Ml/s, and coverage c=20% on the grid of 26

200×200, the simulation results with different interruption times are shown in figure 3.4(a)-(f). It is noticed clearly that with increasing ti, the island size increases; however, when ti=100s, the size reaches the equilibrium, and we finally obtain perfect island ordering when ti=200s. In other words, with increasing interruption time, the system ultimately reaches the equilibrium.

(a) T=550K

(d) T=700K

(b) T=600K

(c) T=650K

(e) T=750K

(f) T=800K

Figure 3.1 Atom island ordering under different temperatures (T is equal to 550K in (a), 600K in (b), 650K in (c), 700K in (d), 750K in (e), and 800K in (f)). Simulation parameters are F=1.0Ml/s, c=20% and ti=200s on a 200×200 grid. The strain field is updated every 2500 steps.

27

(a) F=1Ml/s

(b) F=0.1Ml/s

(c) F=0.01Ml/s

Figure 3.2 Island ordering of atoms under different flux rates. Simulation parameters are T=750K, c=20% and ti=200s on a 200×200 grid. Strain field is NOT included. Deposition stops after 0.2s in (a), 2s in (b), and 20s in (c).

(a) c=10%

(b) c=20%

(d) c=40%

(c) c=30%

(e) c=50%

Figure 3.3 Atom island ordering with coverage c increasing from 10% to 50% in (a) to (e). Simulation parameters are T=700K, F=1.0Ml/s, ti=200s and on a 200×200 grid. The strain energy field is updated every 2500 steps. 28

(a) ti=0

(d) ti=100s

(b) ti=10s

(c) ti=50s

(e) ti=150s

(f) ti=200s

Figure 3.4 Atom island ordering for different interruption times. The ordering of (a) is the initial atom pattern when all the atoms just deposit to the surface and the interruption time ti=0. Simulation parameters are T=750K, F=1.0Ml/s and c=20% on a 200×200 grid. The strain energy is updated every 2500 steps.

29

CHAPTER IV

ELASTIC ANISOTROPY DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION

4.1 Strain Field in Quantum Dot Growth

In self-organized QD research, both theoretical [Ter96, Liu99] and experimental [Tei96, Dar97, Hol98] works have shown that the interaction among dots via their elastic strain field [Xie95, Sol96] may lead to a gradual improvement in size homogeneity, as well as to a more uniform lateral island spacing. Holy and coworkers have further shown that the vertical and lateral correlations in self-organized QD superlattices can be explained by taking into account the elastic anisotropy of the material [Hol99]. In these superlattices, dot correlations are induced by the interaction among dots via their elastic strain fields. Above the buried islands, the elastic energy distribution on the surface exhibits pronounced minima and maxima in the lateral directions. The continuum finite element analysis and many other numerical methods were proposed to study the effect of elastic anisotropy on the three dimensional (3D) strain field, as well as on the layer-tolayer dot correlations to the growth direction [Lee04, Que03, Spr02]. Although QDs ordering due to the elastic anisotropy has been studied [Spr02], the relationship 30

between elastic anisotropy and QD array lateral distribution, the very important character of QD self-organization, has not been investigated. Yet, depending on the elastic properties and growth orientation, a variety of lateral island ordering can be found, which is potentially correlated to the high elastic anisotropy [Sri91, Hol99]. The kinetic Monte Carlo (KMC) has been proposed to simulate QD self-organization by many researchers [Ter96, Joy04, Sri93, Mei03, Tha04]. It has been shown that the main parameters affecting QD island distribution pattern are the temperature T, flux rate F to the surface during deposition, surface coverage c, and growth interruption time ti [Mei03, Pan05]. The simulated island ordering and narrow size distribution in two dimensions were discussed and optimized growth parameters were identified in the former chapter. However, how the substrate anisotropy affects the island pattern has not been studied so far within the KMC simulation algorithm. Therefore, in this chapter, we use a kinetic Monte Carlo (KMC) self-organization model to show how the elastic anisotropy affects the QDs lateral ordering. In our simulation, the atom hopping rate is governed by the Arrenius law modified with the long-range strain energy field. The elastic strain field induced by the islands is obtained using the Green’s function solution for the anisotropic semiconductor substrate developed previously [Pan02a, Pan02b]. QDs array patterns without elastic strain energy, with elastic isotropic strain energy, and with anisotropic strain energy in three different growth orientations, i.e., GaAs (001), GaAs (111) and GaAs (113), are studied. Simulation examples with different growth parameters involving temperature T, coverage c, and interruption time ti are discussed, and our results show that without strain energy, isolated but random QDs nucleation can be

observed. With consideration of the strain 31

energy for isotropic and anisotropic cases, ordered and different island chains can be observed, which proves that strain energy affects the island ordering greatly. Results with different growth parameters show that temperature, coverage and interruption time affect the island ordering quality, but none of them affects island ordering tendency. Thus, the elastic strain energy seems to be the only factor that controls the self-organized island pattern. Furthermore, our simulation results also qualitatively agree with recent experimental results obtained via scanning tunneling microscopy (STM) [Bru95, Bru98].

4.2 Anisotropic Strain Energy Field Calculation

As we mentioned in chapter 2.3.2.2, the elastic strain field induced by the islands can be obtained using the Green’s function solution for the anisotropic semiconductor substrate developed previously and the non-zero elastic coefficients for GaAs (001) are calculated by taking C11=118.8×109N/m2, C12=53.8×109N/m2, and C44=59.4×109N/m2 [Pan02b]. Here, an isotropic model, called Iso (001) is also included for the purpose of comparison in which the elastic constants of GaAs (001) are replaced with its isotropic constants. In the isotropic Iso (001) model, the elastic constants C12 and C44 are the same as those of GaAs (001), i.e., C12=53.8×109N/m2 and C44=59.4×109N/m2, with C11=172.6×109N/m2 obtained by imposing the isotropic condition C11-C12-2C44=0 [Pan03]. 4.2.1 Material Properties of Iso(001), GaAs(001), GaAs(111) and GaAs(113) The elastic moduli of GaAs (111) and GaAs (113) are transformed from GaAs (001) 32

[Pan02a]. The QDs grow along z axis of substrate, which means, for 2D growth, QDs self-organize on xy plane. For GaAs (111) with the substrate coordinates of x axis along (11-2), the y axis along (-110) and the z axis along the (111), for GaAs (113) with the substrate coordinates of x axis along (33-2), the y axis along (-110) and the z axis along the (113), the global material properties of GaAs (111) and GaAs (113) are transformed by the well-known tensor transformation [Nye85] from GaAs (001). Each transformation is done by two steps. For GaAs (111), first, anticlockwise rotate xy surface along z axis

π 4 (from positive direction of z axis), then, anticlockwise rotate xz surface along y axis

(

π 2 − cos −1 2 3

)

(from positive direction of y axis). For GaAs (113), first,

anticlockwise rotate xy surface along z axis π 4 (from positive direction of z axis), then,

(

anticlockwise rotate xz surface along y axis tan −1 3

)

2 (from positive direction of y

axis). Figure 4.1 shows the coordinates of GaAs(001), (111) and (113) and the elastic constant matrixes of Iso (001), GaAs (001), (111) and (113) are shown in equations (4.14.4). Material properties of Iso(001): 0 0 0 ⎤ ⎡172.6 53.8 53.8 ⎢ 53.8 172.6 53.8 0 0 0 ⎥⎥ ⎢ ⎢ 53.8 53.8 172.6 0 0 0 ⎥ C = 10 9 × ⎢ ⎥ Pa 0 0 59.4 0 0 ⎥ ⎢ 0 ⎢ 0 0 0 0 59.4 0 ⎥ ⎢ ⎥ 0 0 0 0 59.4⎥⎦ ⎢⎣ 0 Material properties of GaAs (001):

33

(4.1)

0 0 0 ⎤ ⎡118.8 53.8 53.8 ⎢ 53.8 118.8 53.8 0 0 0 ⎥⎥ ⎢ ⎢ 53.8 53.8 118.8 0 0 0 ⎥ C = 10 9 × ⎢ ⎥ Pa 0 0 59.4 0 0 ⎥ ⎢ 0 ⎢ 0 0 0 0 59.4 0 ⎥ ⎢ ⎥ 0 0 0 0 59.4⎥⎦ ⎢⎣ 0

(4.2)

Material properties of GaAs (111): 45 36 0 12.73 0 ⎤ ⎡ 145 ⎢ 45 − 12.73 145 36 0 0 ⎥⎥ ⎢ ⎢ 36 36 154 0 0 0 ⎥ C = 10 9 × ⎢ ⎥ Pa − 12.73⎥ 0 0 41 0 ⎢ 0 ⎢12.73 − 12.73 0 0 41 0 ⎥ ⎢ ⎥ 0 0 − 12.73 0 50 ⎦⎥ ⎣⎢ 0

(4.3)

Material properties of GaAs (113): − 4.72 0 0 ⎤ ⎡152.81 31.79 41.79 ⎢ 31.79 145.7 − 10.38 48.91 0 0 ⎥⎥ ⎢ ⎢ 41.79 48.91 135.70 0 15.09 0 ⎥ C = 10 9 × ⎢ ⎥ Pa − 0 0 0 54 . 51 0 10 . 38 ⎢ ⎥ ⎢ − 4.72 − 10.38 15.09 0 47.39 0 ⎥ ⎢ ⎥ − 10.38 0 0 0 37.39 ⎦⎥ ⎣⎢ 0

(4.4)

4.2.2 Strain Energy Field of Iso(001), GaAs(001), GaAs(111) and GaAs(113) For the strain energy field calculation, due to its rapid decay behavior, the strain energy field is not evaluated over the whole system but only over a circular area of a given radius around the field point, which is chosen to be thirty grids [Pan05]. This decay feature can be also observed from figure 4.2 where the elastic strain energy field on the surface of substrates Iso (001), GaAs (001), GaAs (111) and GaAs (113) due to a buried point QD are plotted. On the other hand, figure 4.2 also shows clearly the direct 34

correlation between strain energy contour shapes and the growth orientations. The different strain energy patterns on the surface of Iso (001), GaAs (001), GaAs (111) and GaAs (113) will result in island distributions with different inclined island chain orientations. 4.3 Effects of Elastic Strain Energy on Island Ordering

As we have discussed in 2.3.2, during the formation of islands on the substrate, an elastic strain field is also generated around the islands. The nearest and next nearest neighbor binding energies contribute to the shape of the compact islands; however, the island-to-island interaction is very weak without strain energy. As the strain energy is calculated based on the point-eigenstrain Green’s function, it is expected that the absolute value of strain energy is the largest at the boundary area of an island, and with increasing distance away from the island it decays rapidly [Pan02a, Pan02b, Dar97]. Thus, the interaction among islands could lead to island lateral and vertical ordering due to the strain field. This is why that, without elastic strain energy, neither uniform size nor ordered spatial ordering can be observed and one would expect Ostwald ripening [Els03, Zhd99]. To demonstrate the importance of the strain energy and the effect of elastic anisotropy on QD arrays, we take the semiconductor material GaAs as an example. Island ordering of GaAs without strain energy is shown in figure 4.3e, and with strain energy on the surface of substrates Iso (001), GaAs (001), GaAs (111), and GaAs (113), respectively, in figure 4.3(a), (b), (c), and (d). The simulation parameters are temperature T=750K, flux rate F=1.0Ml/s, surface coverage c=20%, and interruption time ti=200s. To 35

show the island distribution more clearly, the island density contours are plotted in figure 4.4 by dividing the original 200×200 small grid size into 18×18 big grid size. So the big grid centers are used as the new coordinates and adatom number of each grid is used as the height. From figure 4.3, it can be observed clearly that without strain energy, isolated but randomly ordered islands can be formed (figure 4.3e) often associated with Ostwald ripening [Zhd99]; however, with the strain energy, isolated and ordered islands (for both size and distribution) can be formed (figure 4.3a,b,c,d). Furthermore, for different growth orientations GaAs (001), GaAs (111) and GaAs (113), different island arrays can be observed (figure 4.3b,c,d), which correspond to the different strain energy patterns shown in figure 4.2 due to material anisotropy with different growth orientations. These results indicate that the long-range strain energy could be one of the very important factors controlling ordered QD arrays and ordered island size. Actually, images of island inclinations on different anisotropic substrates from experiments have also shown the direct correlation between the QD island pattern and the substrate anisotropy (or growth orientation), as by Jacobi [Jac03] where InAs QDs were grown on GaAs (001) and GaAs (113) substrates and by Brune et al. [Bru95, Bur98] where the influence of elastic strain on Ag self-diffusion and nucleation was investigated experimentally. Besides the long-range strain energy field, other key growth parameters include the temperature, coverage, and interruption time, as have been discussed in the previous chapter. However, once we have selected the optimal parameters for temperature, coverage, and interruption time, the island ordering and QD pattern will be mostly controlled by the growth orientation, namely by the anisotropic elastic energy in 36

the substrate with material anisotropy. Therefore, in what follows, we will study these growth parameters by using different growth orientations for the semiconductor GaAs.

4.4 Effect of Substrate Anisotropy and Growth Parameters on Island Ordering

Numerous previous studies have shown that without strain energy field, one could obtain isolated but randomly ordered islands, which are caused by the nearest and next nearest binding energies among adatoms and substrate atoms; however, when the strain energy field is included, isolated and ordered island array (both in size and distribution) can be achieved. Brief results presented in the previous section have further shown that isotropic and anisotropic substrates with different growth orientations could result in different island chain orientations. Furthermore, in the previous chapter, we studied in detail the optimized growth parameters (T, F, c, and ti) for the substrate GaAs (001) and our simulation indicated that these growth parameters could greatly affect the island size and ordering. Here we would like to further study, for the anisotropic substrate case, how these parameters affect the island chain orientation. In other words, is anisotropic elastic strain energy the key parameter in controlling the island chain orientation? To answer this critical question, we have selected four different substrates in our simulation: isotropic substrate Iso (001) and three GaAs substrates but with different growth directions, namely, GaAs (001), GaAs (111) and GaAs (113).

4.4.1 Temperature T & Substrate Anisotropy Simulation

results

for

different

temperatures are shown in figure 4.5 with 37

fixed flux rate F=1Ml/s, coverage c=20%, and interruption time ti = 200s. We selected three different temperatures in the simulation, i.e., T=550K, T=650K, and T=750K, with the results on the substrate isotropic Iso (001) in figure 4.5a and on the anisotropic GaAs (001) in figure 4.5b, GaAs (111) in figure 4.5c and GaAs (113) in figure 4.5d. From figure 4.5, we observed that for all the substrate cases with low temperature (T=550K), no clear pattern can be formed. This is due to the fact that most adatoms do not have enough thermal energy to diffuse, and therefore, the interaction between islands is very weak and no island chains can be formed. With increasing temperature, the island size increases and the island chains become clear, which is also in agreement with experimental results [Sai99]; however, if the temperature is still below the critical value (around 750K as has been identified [Pan05]), the island distributions could be still irregular with no clear correlation to the growth directions (figure 4.5, the middle vertical column for T=650K). On the other hand, if the temperature is too high, the high thermal energy will make the big island unstable, as was discussed before [Pan05]. Therefore, only around the optimized temperature (i.e., T=750K), can one clearly observe ordering and regular islands distribution. More importantly, different island ordering patterns are closely correlated with different growth orientations. Therefore, temperature could affect the island size and island ordering, but it can not control the island chain orientation; the island pattern is directly dominated by the substrate anisotropy or the growth orientation once the optimal temperature is selected. Figure 4.6 and figure 4.7 show the variation of number of islands and average island size vs. different temperatures. Other fixed growth parameters are: flux rate F=1Ml/s, coverage c=20%, and interruption time ti =

200s. We observed from figure 4.6 that, 38

under the given growth conditions, the number of islands and average island size are insensitive to the substrate anisotropy (or growth orientation). For all the substrate, with increasing temperature to 750K, the island size increases and the number of islands decrease. In our previous discussion, we further identified by looking at the island size ordering and distribution that the optimal temperature was around 750K. The island distributions could be still irregular if temperature is too low (below this critical value) and be unstable if temperature is too high (above this critical value). Therefore, only around the optimized temperature (i.e., T=750K), can one clearly observe ordered and regular islands distribution. To further confirm our previous observation, we show here also the variation of the average island size and number of islands versus temperature for a large temperature range (figure 4.6 for GaAs (001) only). It is observed from figure 4.7 that the average size of the island reaches a maximum at T=750K where the number of islands reaches a minimum.

4.4.2 Surface Coverage c & Substrate Anisotropy Simulation results of island formation on substrates Iso (001), GaAs (001), GaAs (111) and GaAs (113) with fixed parameters of temperature T=750K, flux rate F=1Ml/s and interruption time ti = 200s are plotted in figure 4.8 for different coverage c=10%, c=20%, c=30%, and c=50%. It is observed from figure 4.8 that, if the coverage is relatively low, for example for coverage from c=10% to c=20%, the island size doesn’t increase much with increasing coverage, which is in agreement with the experiment results by Leonard et al [Leo94]. On the other hand, if the coverage is too high, for example at c=50%, chain-like islands, instead of isolated islands, are formed. The 39

most important feature is that, whether they are isolated or continuous chain-like islands, different chain orientations can be observed and are clearly correlated to the different growth directions. Consequently, coverage can only control whether the island is isolated or chain-like, whilst the substrate anisotropy controls the orientation of the island array or island-chain. Figure 4.9 shows the variation of the number of island and average island size on the surfaces of the substrates Iso (001), GaAs (001), GaAs (111) and GaAs (113) for different coverage c. The other fixed parameters are temperature T=750K, flux rate F=1Ml/s and interruption time ti = 200s. It is interesting that with increasing coverage, while the number of islands in the grid is nearly a constant value (around 50), the average island size increases, with GaAs (113) having the largest value followed by GaAs (111). The average island size on the surface of the substrate GaAs(001) has the smallest value.

4.4.3 Growth Interruption Time ti & Substrate Anisotropy We now investigate the effect of the growth interruption time ti on the island formation on the surface of four different substrates Iso (001), GaAs (001), GaAs (111), and GaAs (113) with fixed parameters T=750K, F=1.0Ml/s, and c=20% (figure 4.10). The growth interruption time ti is taken as 0s, 10s, and 200s. The results in figure 4.10 clearly show that the strain energy affects the QD arrays and the substrate growth orientation controls the island ordering pattern even from the earlier interruption time say at 10s. It is also observed that, with increasing interruption time, the island inclinations become clearer and clearer, and finally come to the final pattern for ti=200s (figure 4.10). Therefore, suitable amount of interruption time is required for the atom islands to 40

form ordered size and spatial distribution. Obviously, however, long interruption time can only make island size/distribution regular, while the island-chain inclination or island pattern is controlled by the substrate anisotropy or growth orientation. Figure 4.11 shows the effect of the growth interruption time ti on the number of islands and average island size. The fixed growth parameters are T=750K, F=1.0Ml/s, and c=20%. It is observed clearly from figure 4.11 that, with increasing interruption time, the number of islands decreases whilst the average island size increases. Furthermore, at large interruption time (ti=200s for our case), while the number of islands is nearly the same for different substrate orientations (or different material anisotropies), different substrate orientations predict slightly different island sizes with GaAs (111) having the smallest size among all the substrates.

4.5 Stacks of Quantum Dots

Having shown the effect of temperature, flux rate, surface coverage, and growthinterruption time on the island formation (where only a single layer over the substrate is concerned), we now study the interaction between atom layers. As before, the important contribution of the strain energy field will be emphasized. To increase the uniformity and ordering of the QD density over the surface, various approaches have been proposed, including growth of QDs over patterned substrates. Strained-layer heteroepitaxy has been studied in lateral as well as vertical direction [Moi94, Leo94, Not96]. Recently, it has been demonstrated that regular threedimensional QD superlattices with tunable lattice constants can be achieved utilizing 41

strain-mediated self-organization effects in Stranski-Krastanov growth mode [Spr00, Lee04]. There is evidence that stacks not only increase the overall QD density, but also improve the size distribution and the spatial ordering [Dar99]. In this section, we demonstrate the spatial correlation in the vertical direction by carrying out the simulation for a three-layer model with temperature T=700K, flux rate F=0.01Ml/s, coverage c=20%, interruption time ti=50s, and thickness for each layer H=5 grids (figure 4.12). The simulated results for the first three QD layers are shown in figure 4.13. In this model, we first deposit an initial atom pattern similar to figure 4.13a on the substrate GaAs (001). Then our KMC is executed and the simulated pattern for the first layer is shown in figure 4.13a. Next, on top of the first layer, we randomly deposit an initial pattern similar to that of the first layer with distance between the two layers being 5 grids. Simulation similar to the single layer case is then carried out. However, when we calculate the strain energy field at every 2500 hopping steps on the new surface, we consider not only the strain contribution from the surface as we did for figure 4.13a but also that from the substrate (first layer) due to the previous formed atomic pattern. The simulation again runs to the interruption time ti=50s and the result is plotted in figure 4.13b. Finally, following the same procedure, the third layer is deposited over the second layer, and the strain energy at each atom location in the third layer is evaluated by considering all the contributions of atoms in layer 1 to layer 3, with result being shown in figure 4.13c. From figure 4.13, we observed that, with increasing number of deposited layers the island size tends to become larger and the corresponding spatial distribution become more uniform [Spr00, Spr01, Ter96]. This phenomenon is understandable as the 42

additional strain energy from the buried layers enhances the mobility of the diffusing atoms and at the same time, controls the distance at which the atoms can hop. Another interesting factor that affects QD distribution is the thickness of the buffer (or spacer) layer. We have run our KMC program for different layer thickness, with results in figure 4.14b to figure 4.14f corresponding to layer thickness H equal to 5, 10, 15, 20, and 30 grids. It is observed that a thinner layer corresponds to a better island distribution [Ter96] since the thicker the layer, the smaller the strain energy could contribute from the buried layers.

43

z (001)

z (113)

z (111) y (010)

y (-110)

y (-110)

x (100)

x (33-2)

x (11-2)

GaAs (001)

GaAs (113)

GaAs (111)

Figure 4.1 Schematic illustration of crystal coordinates

5

(a)

5 4.4 4 3.2

1

2.8 2.4 2

-1

1.6

1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

3

3.6

x(100) Grid

x(100) Grid

3

(b)

1 -1

1.2

-3

0.8

-3

0.4

-5

-5

0

-5

-3

-1

1

3

5

-5

-3

y(010) Grid

(c)

5 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

x(11-2) Grid

3 1 -1 -3 -5 -5

-3

-1

1

y(-110) Grid

1

3

5

3

5

(d)

4 3.6

3 x(33-2) Grid

5

-1

y(010) Grid

3.2 2.8

1

2.4 2

-1

1.6 1.2

-3

0.8 0.4

-5

0

-5

-3

-1

1

3

5

y(-110) Grid

Figure 4.2 Contours of normalized elastic strain energy distribution on the substrate surface of isotropic Iso(001) (a), and anisotropic GaAs(001) (b), GaAs(111) (c), and GaAs (113) (d). 44

Figure 4.3 Island ordering on the surface of GaAs without elastic strain energy (a), with strain energy of Iso (001) (b), with anisotropic strain energy of GaAs (001) (c), GaAs (111) (d), and of GaAs (113) (e). Different island chain orderings (red dashed lines) corresponding to different growth orientations can be observed. Other simulation parameters are T=750K, F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid.

45

Figure 4.4 Island density contours. The whole 200×200 grids are divided into 18×18 big grids, so there is 11×11 original grids in one big grid. The densities are plotted by 3 dimensional contours.

46

(a) GaAs (Iso)

(b) GaAs (001)

(c) GaAs (111)

(d) GaAs (113)

T=550K

T=650K

T=750K

Figure 4.5 Surface island distributions vs. temperatures (550K, 650K, 750K) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations.

47

Figure 4.6 Number of islands and average island size vs. temperature for islands growth on substrates with different orientations for short temperature range (550K to 750K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s.

Figure 4.7 Number of islands and average island size vs. temperature for islands growth on substrate GaAs (001) for long temperature range (550K to 950K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s. 48

(a) GaAs (Iso)

(b) GaAs (001)

(c) GaAs (111)

(d) GaAs (113)

c=10%

c=20%

c=30%

c=50%

Figure 4.8 Surface island distributions vs. coverage c (10%, 20% 30% and 50%) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations.

49

Figure 4.9 Number of islands and average island size vs. coverage for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and ti=200s.

50

(a) GaAs (Iso)

(b) GaAs (001)

(c) GaAs (111)

(d) GaAs (113)

ti=0s

ti=10s

ti=200s

Figure 4.10 Surface island distributions vs. interruption times ti (0s, 10s, and 200s) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and c=20%, on a 200×200 grid. Red dashed lines indicate the island chain orientations.

51

Figure 4.11 Number of islands and average island size vs. interruption time for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and c=20%.

Third layer Second layer First layer Substrate

Figure 4.12 Schematic illustration of superlattice.

52

(a) First layer

(b) Second layer

(c) Third layer

Figure 4.13 Atom island ordering on the first three layers. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Thickness of each layer is 5 grids. Strain energy field is calculated at every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for growth interruption.

(a) First layer

(b) H=5 grids

(d) H=15 grids

(e) H=20 grids

(c) H=10

(f) H=30

Figure 4.14 Atom island ordering for different layer thickness. Figures (b) to (h) are the simulated results on the second layer with thickness equal to 5 grids to 30 grids. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Strain energy field is updated every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for interruption.

53

CHAPTER V

KINETIC MONTE CARLO SIMULATION OF THREE DIMENSIONAL QUANTUM DOT SELF-ASSEMBLY

5.1 Overview of Three Dimensional Quantum Dots Growth

Self-assembled quantum dots (QDs) are intensively studied in recent years due to their optical and electronic properties with potential applications in optoelectronics and semiconductor devices [Bre95, Bim98]. InAs/GaAs heterostructure is a typical example, which is characterized by a large lattice mismatch and undergoes a growth transition from a two-dimensional (2D) cluster to three-dimensional (3D) islands [Sta04, Joy04]. Their optical properties in the 2D cluster and 3D islands were also studied [Pol06]. Experimental data indicate that InAs/GaAs island array is governed by thermodynamics [Shc99, Pch97]. The physics behind the transition from 2D cluster to 3D island has been studied recently. For example, using the atomic force microscopy, Arciprete et al.[Arc06] studied how kinetics drives the 2D to 3D transition in InAs/GaAs(001) heterostructure. The transition from thermodynamically to kinetically controlled QDs self-assembly was also studied by Musikhin et al. [Mus05], both 54

experimentally and theoretically. It is argued that strains on the surface could change diffusivity and surface mobility [Zoe00, Sag00, Zan99] and therefore could be the reason for the enhanced evaporation [Sun00]. The epitaxial system can lower its free energy by transferring atoms from the island edge to the upper layer, because the transition leads to a decrease in the contact area between the substrate and a new layer [Arc06]. As such atomic transitions to the upper layers lead to the strain field relaxation. The tradeoff between the cost of additional surface energy and the gain of energy due to elastic relaxation is the very driving force for the transition from 2D cluster to 3D islands [Ker97, Cha00]. These arguments provide reason enough to consider situations when the atom hop probability to an upper layer can be higher than that to a lower layer. The kinetics of such a process could be described by adjusting the probability for atoms to hop up or hop down, which could subsequently lead to 2D cluster or 3D island self-assembly [Arc06, Pol06]. While some simple computational models for 2D to 3D self-assembly were discussed before (i.e., [Bru01]), no complete 3D QDs self-assembly model has been developed so far in which the growth from 2D cluster to 3D islands and the corresponding island size equalization can be clearly illustrated during growth process. We, therefore, propose a fast multiscale 3D QDs growth model. It is developed from our original (x,y)-plane growth model [Pan04, Pan06] by a fast algorithm for the long-range strain energy calculation and by introducing the up/down ratio for the lateral self-organization. With the proposed program, we have successfully simulated the transition of the QDs 2D cluster to 3D island process. Furthermore, depositions with different flux rates are studied to show how the flux rate affects the island equalization during deposition 55

process. Our model has also shown clearly the importance of the interruption time during the growth of QD islands. We believe that this model will provide an attractive means for producing ordered nanostructures.The 2D to 3D self-assembly was discussed by some computer modeling [Bru01]. But, so far, no complete 3D QDs self-assembly model has been developed to deal with 2D to 3D transition and island size equalization during growth process. Here, we developed a 3D QDs growth model from our original (x,y)plane growth model [Pan04, Pan06] to study QDs self-assembly in the real 3D growth situation. By adjusting the up/down ratio ρ and considering the lateral self-organization, we simulate the QDs 2D to 3D transition process and island size equalization. Furthermore, examples with different deposition flux rate are studied to show how it affects the island equalization during deposition process. We believe that this model will provide an attractive means for producing ordered nanostructures.

5.2 Three Dimensional Model Description

The 3D layer-by-layer growth model is developed from our 2D (x,y)-plane growth model [Pan04, Pan06]. As in our previous model, the important contribution from the long-range elastic strain energy is included using a fast algorithm based on the Green’s function method [Pan02]. Furthermore, to account for the out-of-plane (up and down) movement of the atoms, which is also called desorption and adsorption [Lam01, Lun05], an up/down ratio is introduced. First, the 2D hopping probability of an atom from one lattice site to a nearest or next nearest neighbor site in the (x,y)-plane is still governed by the Arrhenius law 56

enhanced by the long-range strain energy field [Nur01, Lar01, Pan04, Pan06]. ⎛ E + E n − Estr ( x, y ) ⎞ ⎟⎟ p = ν 0 exp⎜⎜ − s k BT ⎠ ⎝

(5.1)

where ν0 is the attempt frequency (=1013s-1), T the temperature, and kB the Boltzmann’s constant. Also in equation (1), Es and En are the binding energies to the surface and to the neighboring atoms, respectively. Finally, Estr(x,y), as a function of the plane coordinates (x,y), is the energy correction from the long-range strain field due to the lattice mismatch between the substrate and the deposited material. As this long-rang strain energy needs to be calculated repeatedly during the simulation, we have developed a fast algorithm by pre-calculating the energy along a unit circle interpreting the energy at any location afterwards [Pan03]. We calculate the binding energy to the neighboring atoms in the (x,y)-plane, En, using the following algorithm: We assume that the strength of a single nearest neighbor bond is Eb (=0.3eV), and it is reduced by a factor of α (= 1

2 ) for the next nearest

neighbors. To evaluate the diffusion barrier, the binding energy at the site P0, where the diffusing atom is located, is calculated to be EP0 = nEb + αmEb

(5.2)

with n≤4 and m≤4 being, respectively, the number of nearest and next nearest atoms. Similarly, for the site P1 where the atom is going to hop to, we have EP1 = g (n' Eb + αm' Eb )

(5.3)

where n′≤4 and m′≤4 are, respectively, the number of nearest and next nearest atoms at the new site P1, and g (=0.2) describes the coupling between adjacent lattice sites [Mei03]. 57

Therefore, the overall binding energy En caused by neighbor interactions for a diffusion process from site P0 to site P1 in the (x,y)-plane is given by the difference of the binding energies at the corresponding lattice sites E n = (n − g ⋅ n')Eb + (m − g ⋅ m')αEb

(5.4)

Second, the binding energy to the surface of the (x,y)-plane, i.e. Es, was assumed to be constant (actually, Es=1.3eV) in our previous 2D (x,y)-plane self-assembly model [Pan04, Pan06]. However, in 3D case (figure 5.1), the binding energy to the surface depends on the surface geometry of both the initial and final position of the atoms. Based on the recent molecular dynamics calculation [Mon01] and KMC simulation experience [Mei03, Pan04, Pan06] in 2D, we therefore propose the following simple equation for the 3D surface binding energy calculation. Es = (1 − g )Es0 + ( p − g ⋅ p')αEs0 + (q − g ⋅ q')α ⋅ αEs0

(5.5)

where Es0 is the contribution from the atom exactly under the selected adatom (figure 5.1). p, q are number of nearest and next nearest atoms in original positions (n'≤4, m'≤4). p', q' are number of nearest and next nearest atoms in new positions (n'≤4, m'≤4). The other parameters, such as g and α, are used to decrease the contribution from the atoms in the first and second squares (figure 5.1a, 5.1b). Assuming that the maximum surface binding energy be 1.3eV as before for the 2D growth simulation [Pan04, Pan06], which means that all the positions under this adatom are occupied by atoms, we can find Es0=0.28eV from equation (5.5) by back-calculation. Finally, as discussed in the introduction, the growth system always tries to decrease its free energy by moving atoms at the edge to upper layer to form 3D islands. The possibility for an atom to hop to upper or lower layers depends on the material 58

properties and growth condition [Arc06, Pol06]. To account for this up and down jump activities of the atoms, the parameter up/down ratio ρ is introduced, which is defined as the ratio of the probability for edge atoms hopping up to the probability hopping down. In other words, if we let Pup be the jumping up possibility and Pdown the jumping down possibility, then the up/down ratio ρ=Pup/Pdown. It is further remarked that the up/down ratio actually reflects the balance between surface energy increase and strain energy decrease. This physical activity is illustrated in figure 2 where layer 1 is the substrate and the atom “A” is on the edge of layer 2 (within the QD island). It is also understandable that the up or down jump possibility of atom “A” is controlled only by the up/down ratio instead of by the individual jumping up and down possibilities. Moreover, strictly speaking, since the strain energy changes along the island height direction, the up/down ratio is, in general, not constant in different island layers. However, in order to extract the important contribution of the up/down ratio, we assume that the up/down ratio ρ is constant. In other words, it is only affected by the geometry around it but not by the layer position. Furthermore, in order to form 3D islands, the number of atoms jumping up should be larger than those jumping down, which means that the up/down ratio ρ should be larger than 1. Otherwise, all the atoms will move to lower layer to form the layer-bylayer Frank-van der Merwe thin-film structure.

5.3 Simulation Results

Using the 3D QDs self-assembly approach described above, we can now simulate the epitaxial growth process. The growth model is on the 100×100 grids with 59

periodic boundary condition, and the material is InAs/GaAs (001).

5.3.1 Height Dependence of Up/Down Ratio (ρ) We first study the effect of the up/down ratio on the island height. Shown in figure 5.3 are the island distributions for different up/down ratio ρ. It is observed clearly from figure 5.3 that 1) the average island height increases with increasing up/down ratio ρ (for small up/down ratio, say, ρ30); 3) the inflection point within the sharp-change range, approximately at ρ=13, 60

physically corresponds to the critical transition point from 2D growth to 3D growth. This value is equivalent to an energy ratio between surface and bulk as has been demonstrated in [Bru01, Kra06] using a different approach.

5.3.2 Size Dependence of Up/Down Ratio (ρ) As we know, the island size equalization relies on the movement of atoms on the surface during epitaxial growth. The equalization starts to happen from the very beginning of deposition process and finally it will reach the equilibrium status if long enough interruption time after deposition [Pan04] is given. Here, three examples with different up/down ratios of 10, 20 and 30 have been run to show how the island size equalization happens during the deposition and interruption process (figure 5.5). For each example, the island distribution of exactly after deposition (ti=0), 100 seconds after deposition (ti=100) and 200 seconds after deposition (ti=200s) are shown in figure 5.5. It is easy to find that: 1) for ρ=10, 20 and 30, increasing the interruption time, small isolated islands assembles to big ones, and, with long interruption time, the island distribution is getting increasingly ordered; 2) by adjusting the up/down ratio, people can get island distribution with different shape, which could be a way for QDs growth experimentalists to get the optimized QDs distribution or shape by adjusting the up/down ratio. Taking the first example with up/down ratio of 10, in figure 5.5, we analyze the relationship between island diameter and number of island within this diameter (figure 5.6 and figure 5.7) using the simulated results of 0.5 ML, 1 ML, 1.6 ML deposition and 200s interruption after deposition. We find some remarks. 1) Island size equalization 61

happens starting from the beginning of deposition process for all three up/down ratio

ρ=10, 20 and 30, and, with long enough interruption time, finally we will get ordered island distribution with similar size; 2) In this example, low flux rate (F=0.01 ML/s) is used, so island equalization is obvious even just at the end of the deposition process; 3) We believe that if we decrease the flux rate, we will get the better island distribution during deposition process. To further confirm our conclusion in 3), more examples with same growth parameters except for different flux rate are run in figure 5.8. Finally, figure 5.8 shows the effect of the flux rate (F=1Ml/s, 0.1ML/s and 0.01ML/s) on the island distribution. For a lower flux rate (e.g. at F=0.01Ml/s), the deposited atoms will have more time to move to the equilibrium position during the deposition process and to assemble together. Furthermore, a low flux rate (say, F=0.01 ML/s) usually corresponds to large and ordered islands. Therefore, besides the up/down ratio, flux rate is also important in order to obtain ordered and uniform-size QD island distribution.

5.4 Conclusion of 3D Simulation Model

In this chapter, we propose a 3D QDs epitaxial growth model by updating our former (x,y)-plane growth model. Specifically, the balance between surface energy increase and strain energy decrease is demonstrated by introducing the parameter up/down ratio. Furthermore, island height dependence and island shape dependence of up/down ratio have been studied too. Combining up/down ratio and one of the growth parameters flux rate, island equalization during deposition process and after deposition is studied. The main conclusions from our 3D KMC model are as follows: 62

1.

The phenomena of self-assembly and shape transition of the island array can be simulated with changing parameter up/down ratio.

2.

For other fixed growth parameters, island height and size greatly depends on the up/down ratio. For the examples used in this paper, the average island height increases with increasing up/down ratio and island size and shape are very sensitive to the up/down ratio when it is less then 20. Furthermore, the critical transition point of up/down ratio in this example is around 13.

3.

It is the up/down ratio, the ratio characterizes the atom jump possibility to upper layer over that to lower layer, that determines the transition from 2D cluster to 3D islands, instead of the absolute jump possibility value.

4.

Island size equalization starts from the very beginning of deposition process. With enough interruption time, ordered island distribution with uniform size can be obtained.

5.

The lower the flux rate, the better the island distribution during deposition process.

63

(a)

(c)

(b)

Es y

y x x Figure 5.1 Schematic illustration of 3D QD self-assembly model. An atom on top of the substrate surface (x,y)-plane in (a), the relative locations of the atom grid on the (x,y)plane in (b), the corresponding binding energy Es related to the in-plane locations in (c).

4 3 2 1

A

Figure 5.2 Illustration of an edge atom A’s jump up or down during 3D QD self-assembly.

64

(a) ρ=3

(b) ρ=5

(c) ρ=7

(d) ρ=10

(e) ρ=20

(f) ρ=30

Figure 5.3 3D island distributions for different up/down ratios ρ with total coverage c=1.6Ml. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and interruption time ti=200s (The total simulation time = the deposition time (td=160s) plus the interruption time (ti=200s)).

Average Island Height (Grids)

20 16 12 8 4 0 0

10 20 30 40 Up/Down Ratio (ρ)

50

Figure 5.4 Average island height vs. up/down ratio ρ (for ρ from 3 to 50). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and interruption time ti=200s. 65

ti =0

ti =100s

ti =200s

ρ=10

ρ=20

ρ=30

Figure 5.5 Island distributions for different interruption times (ti=0s, corresponding to deposition time td=160s, ti=100s, and ti=200s) and for different up/down ratios (ρ=10, 20 and 30). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and total coverage c=1.6Ml.

66

(a) td=50s

(b) td=100s

(c) td=160s; ti=0

(d) ti=200s

Figure 5.6 Island distributions during and after deposition. Deposition time td=50s with 0.5Ml coverage in (a), td=100s with 1Ml coverage in (b), td=160s with 1.6Ml coverage (i.e., at the beginning of interruption time) in (c) and interruption time ti=200s in (d). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and up/down ratio ρ=10.

67

40

40

(a) td=50s

30 Num of Islands

Num of Islands

30

20

10

0

0

2

0

4 6 8 10 12 14 16 18 Island Diameter (Grids)

0

2

4 6 8 10 12 14 16 18 Island Diameter (Grids)

(d) ti=200s

16

(c) td=160s; ti=0

20

12 Num of Islands

Num of Islands

20

10

25

15 10

8

4

5 0

(b) td=100s

0

2

0

4 6 8 10 12 14 16 18 Island Diameter (Grids)

0

2

4 6 8 10 12 14 16 18 Island Diameter (Grids)

Figure 5.7 Average island diameter (of the bottom equivalent circle of the island) vs. number of islands during and after deposition. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml and up/down ratio ρ=10.

68

0.5 ML

1 ML

1.6 ML

1 ML/s

0.1 ML/s

0.01 ML/s

Figure 5.8 Island distributions for different deposition processes of 0.5Ml, 1Ml and1.6 Ml with flux rates F=1Ml/s, 0.1Ml/s and 0.01Ml/s. Fixed growth parameters are temperature T=700K, total coverage c=1.6Ml and up/down ratio ρ=10.

69

CHAPTER VI

MOLECULAR DYNAMICS CALCULATION OF ELASTIC MODULI IN STRAINED SEMICONDUCOTR

6.1 Background of Elastic Moduli in Strained Semiconductor

The quantum dots growth and the effect of anisotropic strain energy on island distribution have been studied in the former chapters. As we discussed, the anisotropic strain energy is caused by the misfit strain and basically is because of the unit crystal side length difference between the substrate and deposited material. Take InAs/GaAs system as an example [Ell02], usually we take the InAs QD and compress it uniformly by about 6.685% in each direction so that it has the same lattice constant as that of GaAs. It is then placed in a hole of exactly the right size, so that the InAs and GaAs crystals are coherent, and then the uniform compression is removed, which leads to straining of the surrounding matrix. Though it is not easy to accurately determine the strain distribution, there are some methods available [Dav99, Hol99, Pan01]. The strain effects of semiconductors have been widely investigated and lots of reports about influence on electronic and optical properties are available [Gle99, 70

Jai00, Jog98, Par99], but much less is known about their influence on the elastic constants. In the recent years, people begin to study the elastic constants of strained semiconductors. For example, Mehl [Meh93] studied moduli of Al-Li compounds under pressure, and Kanoun et al. studied the elastic properties of BN, AlN, GaN and InN under pressure effect by using ab initio calculations [Kan04a, Kan04b, Mer03]. And, recently, Ellaway and Faux [Ell02, Ell03] reported the effective elastic constants of InAs under uniform strain by using DLPOLY [Smi96] molecular dynamics (MD) package. We develop a MD program by using Stillinger-Weber potential. To check our program, using the parameters for InAs [Ell02, Ich96], elastic stiffnesses of InAs determined by our program match the results by Ellaway and Faux [Ell02] very well. Crystalline silicon is widely used in micro-electro-mechanical-systems (MEMS) and it is one of the most intensely investigated materials. Interest in silicon is mostly motivated by its great technological importance. Strained silicon has its own specific meaningful use. In 2002, Intel announced their employing a unique strained silicon transistor technology on its 90 nm logic process. This was the first process in the industry to implement strained silicon in production. This process improves transistor performance by 10-25 percent while increasing manufacturing costs by only two percent. This is because strained silicon increases transistor drive current which improves switching speed by making current flow more smoothly [Gha03]. So, to understand the mechanism of strained silicon well, to study the elastic moduli of strained silicon is very important. In this chapter, the effective moduli of silicon under hydrostatic strain are calculated as the functions of volumetric strain. While variation of the elastic moduli was studied previously when silicon was under pressure [Mcs64, Kar97], no literature is 71

available on the stress and strain behavior when silicon is under strain. Furthermore, due to the technical difficulty, it is hard to observe experimentally the elastic properties as functions of the applied strain. Here, we present a molecular dynamic prediction for the elastic stiffnesses in strained silicon as functions of the volumetric strain level. Our approach combines the basic continuum mechanics with the classical molecular dynamic approach, supplemented with the interatomic S-W potential [Sti85]. This potential has gained popularity due to its ability to describe fairly well diamond structures and has been used in numerous studies (e.g., [Klu86, Hum04, Men04]). In order to use our approach for the estimation of the three independent elastic stiffnesses C11, C12 and C44 in strained silicon, we apply three independent distortions as described by Mehl [Meh90] and utilized later by Ellaway and Faux [Ell02] for InAs. The elastic stiffnesses and effective Young’s modulus and Poisson’s ratio are all calculated and presented in terms of figures and formulas. In general, our simulation indicates that the bulk modulus, C11 and C12, increase with increasing volumetric strain whilst C44 is almost independent of volumetric strain. The difference between the strained moduli and the ones at zero strain can be very large, and therefore use of the standard free-strained moduli in the strained silicon case could result in wrong mechanical response of the silicon.

6.2 Elastic Stiffness Calculation

The III-V semiconductors are cubic crystals with three independent elastic stiffnesses C11, C12 and C44. An efficient way of calculating elastic constants is 72

important because these constants are directly employed in practical uses of materials. There are several ways to calculate elastic constants. The direct and traditional method is to apply a tension on the sample and calculate the corresponding strain and elastic constants from the tension-strain relationship [Ell02, Meh93, Ray88]. For a cubic crystal, the elastic moduli can be divided into two classes: the bulk modulus B=(C11+2C12)/3 and two “shear” moduli, C11-C12 and C44. S In the actual calculation, a small incremental distortion, which is usually less than 1%, is applied to the crystal with the latter being already subjected to a uniform volumetric strain (i.e., [Ell02]). Let’s assume the original cubic lengths and the distortions in x, y and z directions are εx, εy, εz and ∆x, ∆y, ∆z. After deformation, the new volume is: V = ∆x(1 + ε x )∆y (1 + ε y )∆z (1 + ε z )

= ∆x∆y∆z (1 + ε x + ε y + ε z + ε xε y + ε y ε z + ε z ε x + ε xε y ε z )

(6.1)

For small strains, the high order terms can be neglected and the volume simplifies as: V ≈ ∆x∆y∆z (1 + ε x + ε y + ε z ) = V0 (1 + ε x + ε y + ε z )

(6.2)

Thus, the volumetric strain becomes e=

V − V0 =εx +εy +εz V0

(6.3)

where V0 is the original volume of the unit cell, V is the volume after deformation and e is volumetric strain. Constitutive relation for cubic material is

73

⎡σ xx ⎤ ⎡C11 ⎢σ ⎥ ⎢ ⎢ yy ⎥ ⎢C12 ⎢σ zz ⎥ ⎢C12 ⎢ ⎥=⎢ ⎢ τ yz ⎥ ⎢ 0 ⎢ τ zx ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎣⎢ τ xy ⎦⎥ ⎣⎢ 0

C12

C12

0

0

C11

C12

0

0

C12

C11

0

0

0

0

C 44

0

0

0

0

C 44

0

0

0

0

0 ⎤ ⎡ε xx ⎤ ⎢ ⎥ 0 ⎥⎥ ⎢ε yy ⎥ 0 ⎥ ⎢ε zz ⎥ ⎥⎢ ⎥ 0 ⎥ ⎢γ yz ⎥ 0 ⎥ ⎢γ zx ⎥ ⎥⎢ ⎥ C 44 ⎦⎥ ⎣⎢γ xy ⎦⎥

(6.4)

And energy density expression for this system is 1 U = (σ xx ε xx + σ yy ε yy + σ zz ε zz + τ yz γ yz + τ zx γ zx + τ xy γ xy ) 2

(6.5)

In the following, we briefly list the relevant formulas needed to obtain the elastic constants of cubic crystal. The stress tensor and the elastic constants can be defined in different ways and it is useful to define our notation. We consider both strain with and without volume conservation, because volume-conserving strains alone do not provide enough constraints to determine all three elastic constants of cubic crystal. Only small lattice distortions are considered in order to remain within the elastic domain of the crystal.

6.2.1 Bulk Modulus Bulk modulus is related to the E(V ) curvature [Bir78],

B(V ) = −VP' (V ) = VE " (V )

(6.6)

where V is the volume of unit cell as mentioned before, E(V) is the energy/unit cell at volume V, E″(V) is the second derivative of E(V) with respect to the volume, and P(V) = – E′(V) is the pressure required to keep the cell at volume V. Since the calculations only provide a set of energies E(Vi) for a limited number of volumes Vi, the second derivative

E″(V) must be approximated [Meh93]. Birch [Bir78] proposed the following form 74

to calculate E(V) by fitting the cure. 2

3

N ⎡⎛ V0 ⎞ 2 / 3 ⎤ ⎡⎛ V0 ⎞ 2 / 3 ⎤ ⎡⎛ V0 ⎞ 2 / 3 ⎤ 9 9 ' E (V ) = E0 + B0V0 ⎢⎜ ⎟ − 1⎥ + B0V0 B0 − 4 ⎢⎜ ⎟ − 1⎥ + ∑ γ n ⎢⎜ ⎟ − 1⎥ 8 ⎥⎦ 16 ⎢⎣⎝ V ⎠ ⎥⎦ n=4 ⎣⎢⎝ V ⎠ ⎥⎦ ⎣⎢⎝ V ⎠

(

)

n

(6.7) where E0, V0, B0, and B0′ are, respectively, the equilibrium energy, volume, bulk modulus, and pressure derivative of the bulk modulus, whilst N is the curve fitting order. We remark that equation (6.7) is a special case of the following general expression N

E (V ) = ∑ anV −2 n / 3

(6.8)

n =0

where an are the fitting parameters and we find that a third polynomial is sufficient (N=3). In our simulation, the distortion used for calculating bulk modulus consists of a uniform compression (or expansion) of magnitude δ in each of x, y and z directions as in a conventional pressure experiment. The distortion can be represented by the matrix ⎛δ 0 0 ⎞ ⎜ ⎟ ⎜0 δ 0⎟ ⎜0 0 δ ⎟ ⎝ ⎠

(6.9)

6.2.2 Modulus C11-C12 In case of cubic of lattice, it is possible to choose the strain that keeps the volume of unit cell preserved. “Shear” modulus C11-C12 is evaluated by applying an orthorhombic volume-conserving (up to the first order of δ) strain, which can be represented in the strain tensor form as (proved in Appdendix C)

75

⎛ ⎜δ 0 ⎜ ⎜0 −δ ⎜⎜ 0 0 ⎝

⎞ ⎟ ⎟ ⎟ 2 δ ⎟ ⎟ 1−δ 2 ⎠ 0 0

(6.10)

Making use of the stress and strain relation, we found that the non-zero stresses are ⎡σ xx ⎤ ⎡C11 ⎢σ ⎥ ⎢ ⎢ yy ⎥ ⎢C12 ⎢σ zz ⎥ ⎢C12 ⎢ ⎥=⎢ ⎢ τ yz ⎥ ⎢ 0 ⎢ τ zx ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎢⎣ τ xy ⎥⎦ ⎣⎢ 0

C12

C12

0

0

C11

C12

0

0

C12 0

C11 0

0

0 0

0

0

C44 0

0

0

0

C44 0

0 ⎤⎡ δ ⎥ ⎢ 0 ⎥ −δ ⎢ 2 0 ⎥ ⎢δ / 1 − δ 2 ⎥⎢ 0 ⎥⎢ 0 0 ⎥⎢ 0 ⎥⎢ C44 ⎦⎥ ⎣⎢ 0

(

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦⎥

)

(6.11)

Then we have

δ2 1− δ 2 δ2 σ yy = (C12 −C 11)δ + C12 1− δ 2 δ2 σ zz = C11 1− δ 2 τ yz = τ zx = τ xy = 0 σ xx = (C11 −C 12 )δ + C12

(6.12)

So, the system energy as a function of δ is

E (δ ) = E (0 ) + (C11 − C12 )Vδ 2 + 0(δ 4 ) + ...

(6.13)

where terms of order δ4 and higher order can be safely neglected. So, the energy E(δ) may be found as a function of δ2 and a linear curve fitting yields C11-C12.

6.2.3 Elastic Stiffness C44 Elastic stiffness C44 is evaluated by a monoclinic volume-conserving strain (proved in Appendix D). The strain can be represented by the matrix 76

⎛ ⎜ 0 ⎜ ⎜1 ⎜2δ ⎜ ⎜⎜ 0 ⎝

1 δ 2

0

0

0

0

δ2 (4 − δ 2 )

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎟ ⎠

(6.14)

and, the stress and strain relation becomes: ⎡σ xx ⎤ ⎡C11 ⎢σ ⎥ ⎢ ⎢ yy ⎥ ⎢C12 ⎢σ zz ⎥ ⎢C12 ⎢ ⎥=⎢ ⎢τ yz ⎥ ⎢ 0 ⎢ τ zx ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎢⎣τ xy ⎥⎦ ⎣⎢ 0

C12

C12

0

0

C11

C12

0

0

C12

C11

0

0

0

0

C44

0

0

0

0

C44

0

0

0

0

0 ⎤⎡ 0 ⎥⎥ ⎢ ⎢ 0 ⎥ ⎢δ 2 / ⎥⎢ 0 ⎥⎢ 0 ⎥⎢ ⎥⎢ C44 ⎦⎥ ⎣⎢

⎤ ⎥ 0 ⎥ 4−δ 2 ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ δ ⎦⎥ 0

(

)

(6.15)

Then we have

σ xx = σ yy = σ zz = 0 τ yz = τ zx = 0

(6.16)

τ xy = δ So, application of this strain changes the total energy from its unstrained value to

( )

1 E (δ ) = E (0) + C 44Vδ 2 + 0 δ 4 + ... 2

(6.17)

Similarly, the energy E(δ) can be found as a function of δ2 and a linear curve fitting yields C44.

6.3 Molecular Dynamics of Elastic Moduli Calculation

Molecular dynamics is a technique to compute the equilibrium and transport properties of a classical many-body system.

77

6.3.1 Classical Interatomic Potential In order to use molecular dynamics, we have to define the rules that are governing interaction of atoms in the system. In classical and semi-calssical simulations, these rules are often expressed in terms of potential functions. The potential function U(ri) describes how the potential energy of a system of N atoms depends on the coordinates of the atoms. In the past, various empirical potentials have been proposed to describe the different phases of silicon. A number of these empirical potentials [Ter88, Sti85, Bas89, Bis85] have been proposed to describe the different phases of Si and the results have been tested and compared [Cow88, Coo93, Bal92]. Among these, the Stillinger and Weber (S-W) potential [Sti85] has gained popularity due to its ability to describe fairly well diamond structures and has been used in numerous studies [Klu86, Hum04, Men04] since it was proposed. The main advantage of this potential is its simplicity and fairly realistic description of crystal. The potential consists of a two-body term and an explicit angular interaction: V = ∑ν 2 (rij ) +

∑ν (r , r , r ) (r ) = εf (r σ ) (r , r , r ) = εf (r σ , r σ , r i< j

ν2 ν3

ij i

2

j

i< j

A Dissertation Presented to The Graduate Faculty of the University of Akron

In Partial Fulfillment of the Requirements for the Degree Doctor of Philosophy

Ronghua Zhu December, 2006

ATOMISTIC SIMULATION OF NANOSTRUCTURED MATERIALS

Ronghua Zhu Dissertation

Approved:

Accepted:

______________________________ Advisor Dr. Ernian Pan

______________________________ Department Chair Dr. Wieslaw K. Binienda

______________________________ Committee Member Dr. Wieslaw K. Binienda

______________________________ Dean of the College Dr. George K. Haritos

______________________________ Committee Member Dr. Ala R. Abbas

______________________________ Dean of the Graduate School Dr. George R. Newkome

______________________________ Committee Member Dr. Xiaosheng Gao

______________________________ Date

______________________________ Committee Member Dr. Alper Buldum ______________________________ Committee Member Dr. Peter W. Chung

ii

ABSTRACT

Atomistic based computer modeling and simulation of nanostructured materials has become an important subfield of materials research. Based on the multiresolution method, which combines the continuum mechanics, kinetic Monte Carlo method and molecular dynamics method, we study the nanostructured materials grown by quantum-dot selfassembly, mechanical properties of strained semiconductors, and mechanical properties of carbon nanotube reinforced composites. This thesis covers the following three main contributions. 1). Self-organization of semiconductors InAs/GaAs in Stranski-Krastanov growth mode is studied using kinetic Monte Carlo simulations method coupled with the Green’s function solution for the elastic strain energy distribution. The relevant growth parameters such as growth temperature, surface coverage, flux rate, and growth interruption time are investigated. It is shown clearly that when the long-range strain energy is included in the simulation, ordered uniform size distribution can be achieved. To address the effect of material anisotropy, the anisotropic substrates of GaAs with different growth orientations (001), (111), and (113) and an isotropic substrate Iso (001), reduced from cubic GaAs, are also investigated. Simulation results show that at selected growth

parameters

for

temperature, coverage, and growth interruption time, iii

strain energy field in the substrate is the key factor that controls the pattern of island distribution. Furthermore, layer-by-layer growth of quantum dots is also simulated briefly, and vertical alignment is observed that could lead to progressively uniform island sizes and spatial ordering. 2). Since the misfit strain will be induced during the quantum dots epitaxial growth, the mechanical property of the grown semiconductors will be influenced. In this thesis, utilizing the basic continuum mechanics, we present a molecular dynamic prediction for the elastic stiffness C11, C12 and C44 in strained silicon and InAs as functions of the volumetric (misfit) strain. The bulk modulus, effective elastic stiffness C11, C12 and C44 of the strained silicon, including also the effective Young’s modulus and Poisson’s ratio, are all calculated and presented. In general, our simulation indicates that the bulk moduli, C11 and C12 increase with increasing volumetric strain whilst C44 is almost independent of the volumetric strain. 3). Also using MD method, the carbon nanotube reinforced Epon 862 composite is studied. The stress-strain relations and the elastic Young’s moduli along the longitudinal direction (parallel to CNT) are simulated with the results being also compared with those from the rule-of-mixture. Our results show that, with increasing strain in the longitudinal direction, the Young’s modulus of CNT increases whilst that of the Epon 862 composite or matrix decreases. A long CNT can greatly improve the Young’s modulus of the Epon 862 composite (about 10 times stiffer), which is also consistent with the prediction based on the rule-of-mixture at low strain level. Furthermore, even a short CNT can also enhance the Young’s modulus of the Epon 862 composite, with an increment of 20% being observed as compared to that of the Epon 862 matrix. iv

ACKNOWLEDGEMENTS

This thesis is the end of my long journey in obtaining my doctor’s degree in Civil Engineering. I have not traveled in a vacuum in this journey. There are some people who made this journey easier with words of encouragement, with hands of help, with experience of guidance …… To my advisor Dr. Ernian Pan for his help, encouragement and suggestions during the past three and half years! Deeply indebted to Dr. Ernian Pan’s kindness and friendship. His diligence and unremittingness set me an excellent example, which will be helpful to me during my whole life. To my committee members: Dr. Wieslaw K. Binienda, Dr. Ala R. Abbas, Dr. Yu Qiao, Dr. Xiaosheng Gao, Dr. Alper Buldum and Dr. Peter W. Chung for their time, kind support and helpful suggestions. To my colleagues, Feng Han, Wael Alkasawn, Yan Zhang, Xu Wang, Mingkun Sun, Yuanguo Chen, James Ramsey and Kotchanu Jittava, for giving me the feeling of being at home at work. To my parents, my family and friends …… And, to Army Research Office for their support.

v

TABLE OF CONTENTS Page LIST OF FIGURES …….….……………………..…………………………………….. x CHAPTER I INTRODUCTION…….….……………………..……………………………………..... 1 1.1

Atomistic Simulation in Nanostructured Materials........................................... 1

1.2

Multiresolution in Atomistic Simulation........................................................... 2

1.3

Research Purpose .............................................................................................. 4

II MONTE CARLO SIMULATION OF SELF-ORGANIZED QUANTUM DOT GROWTH…………………………………………………………………………

6

2.1

Quantum Dots.................................................................................................... 6

2.2

Self-organized Growth ...................................................................................... 8

2.3

2.2.1

Deposition and Diffusion........................................................................ 9

2.2.2

Growth Classification ........................................................................... 10

Kinetic Monte Carlo Simulation Theory......................................................... 13 2.3.1

Nucleation and Growth of Quantum Dots ............................................ 13

2.3.2

Atomic Hopping Probability................................................................. 16

2.3.3

Kinetic Monte Carlo Algorithm............................................................ 19

III GROWTH PARAMETERS DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION………………………………………………... 23 3.1

Growth Parameters .......................................................................................... 23

3.2

Temperature T ................................................................................................. 24 vi

3.3

Flux Rate F...................................................................................................... 24

3.4

Surface Coverage c.......................................................................................... 25

3.5

Growth Interruption Time ti ............................................................................ 26

IV ELASTIC ANISOTROPY DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION……………………………………………….. 30 4.1

Strain Field in Quantum Dot Growth .............................................................. 30

4.2

Anisotropic Strain Energy Field Calculation .................................................. 32 4.2.1 Material Properties of Iso(001), GaAs(001), GaAs(111) and GaAs(113)......................................................................................................... 32 4.2.2 Strain Energy Field of Iso(001), GaAs(001), GaAs(111) and GaAs(113)......................................................................................................... 34

4.3

Effects of Elastic Strain Energy on Island Ordering ....................................... 35

4.4

Effect of Substrate Anisotropy and Growth Parameters on Island Ordering .. 37

4.5

Stacks of Quantum Dots.................................................................................. 41

V KINETIC MONTE CARLO SIMULATION OF THREE DIMENSIONAL QUANTUM DOT SELF-ASSEMBLY ……………………………………………54 5.1

Overview of Three Dimensional Quantum Dots Growth ............................... 54

5.2

Three Dimensional Model Description ........................................................... 56

5.3

Simulation Results........................................................................................... 59

5.4

Conclusion of 3D Simulation Model............................................................... 62

VI MOLECULAR DYNAMICS CALCULATION OF ELASTIC MODULI IN STRAINED SEMICONDUCOTR………………………………………………….70 6.1

Background of Elastic Moduli in Strained Semiconductor............................. 70

6.2

Elastic Stiffness Calculation............................................................................ 72 6.2.1

Bulk Modulus........................................................................................ 74

6.2.2

Modulus C11-C12 ................................................................................... 75

vii

6.2.3 Elastic Stiffness C44 .............................................................................. 76 6.3

6.4

6.5

Molecular Dynamics of Elastic Moduli Calculation ....................................... 77 6.3.1

Classical Interatomic Potential ............................................................. 78

6.3.2

MD Programming ................................................................................. 79

Elastic Stiffness Calculation of InAs............................................................... 80 6.4.1

Bulk Modulus........................................................................................ 81

6.4.2

C11 and C12 ............................................................................................ 81

6.4.3

C44 ......................................................................................................... 81

Elastic Stiffness Calculation of Silicon ........................................................... 82 6.5.1

Potential and Volumetric Strain............................................................ 83

6.5.2

C11, C12 and C44 ..................................................................................... 84

6.6

Effective Young’s modulus and Poisson’s ratio ............................................. 84

6.7

Conclusion of MD Simulation of Strained Semiconductor............................. 86

VII MOLECULAR DYNAMICS STUDY OF THE STRESS-STRAIN BEHAVIOR OF CARBON-NANOTUBE REINFORCED EPON862 ……………………………..96 7.1

Overview of Carbon Nanotube Composites.................................................... 96

7.2

Molecular Dynamics Simulation of CNT Composites.................................... 97 7.2.1

Model Building ……………………………………………………….98

7.2.2 Simulation of Stress-Strain Relationship ……………………………100 7.2.3 Rule-of-Mixture ……………………………………………………..102 7.3

Numerical Results ......................................................................................... 103

7.4

Concluding Remarks of MD Simulated CNT Reinforced Epon 862 ............ 105

VIII CONCLUSIONS ………………………………………………………………….112 REFERENCES ............................................................................................................. 114 APPENDICES ………………………………………………………………………...122 viii

APPENDIX A HALF-SPACE GREEN’S FUNCTIONS ...................................... 125 APPENDIX B DETAILED SIMULATION STEPS.............................................. 129 APPENDIX C VOLUME-CONSERVING OF ORTHORHOMBIC STRAIN..... 132 APPENDIX DVOLUME-CONSERVING OF MONOCLINIC STRAIN............. 133

ix

LIST OF FIGURES Figure

Page

1.1 Different laws of physics are required to describe properties and processes of fluids at different scales [E03]. ..................................................................................... 5 2.1 Possible ways of diffusion from A to B: a - desorption, condensation process; b surface diffusion; c - volume diffusion....................................................................... 21 2.2 Important growth modes in epitaxy: a) Volmer-Weber growth mode, b) FrankVan der Merwe growth mode, c) Stranski-Krastanov growth mode. ......................... 21 2.3 Strain relief by a) generation of mis-fit dislocations (green atoms) and b) StranskiKrastanov growth mode. The green line marks the end of the wetting layer. Atoms of the substrate are blue, deposited material is red. ........................................ 22 2.4 Schematic illustration of off, edge, and corner diffusions. ......................................... 22 2.5 Schematic illustration of nearest (horizontal and vertical directions in solid blue) and next nearest neighbor (diagonal directions in dashed brown) positions of each atom............................................................................................................................. 22 3.1 Atom island ordering under different temperatures (T is equal to 550K in (a), 600K in (b), 650K in (c), 700K in (d), 750K in (e), and 800K in (f)). Simulation parameters are F=1.0Ml/s, c=20% and ti=200s on a 200×200 grid. The strain field is updated every 2500 steps. ....................................................................................... 27 3.2 Island ordering of atoms under different flux rates. Simulation parameters are T=750K, c=20% and ti=200s on a 200×200 grid. Strain field is NOT included. Deposition stops after 0.2s in (a), 2s in (b), and 20s in (c). ........................................ 28 3.3 Atom island ordering with coverage c increasing from 10% to 50% in (a) to (e). Simulation parameters are T=700K, F=1.0Ml/s, ti=200s and on a 200×200 grid. The strain energy field is updated every 2500 steps. .................................................. 28

x

3.4 Atom island ordering for different interruption times. The ordering of (a) is the initial atom pattern when all the atoms just deposit to the surface and the interruption time ti=0. Simulation parameters are T=750K, F=1.0Ml/s and c=20% on a 200×200 grid. The strain energy is updated every 2500 steps............................ 29 4.1 Schematic illustration of crystal coordinates .............................................................. 44 4.2 Contours of normalized elastic strain energy distribution on the substrate surface of isotropic Iso(001) (a), and anisotropic GaAs(001) (b), GaAs(111) (c), and GaAs (113) (d). ........................................................................................................... 44 4.3 Island ordering on the surface of GaAs without elastic strain energy (a), with strain energy of Iso (001) (b), with anisotropic strain energy of GaAs (001) (c), GaAs (111) (d), and of GaAs (113) (e). Different island chain orderings (red dashed lines) corresponding to different growth orientations can be observed. Other simulation parameters are T=750K, F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid............................................................................................................... 45 4.4 Island density contours. The whole 200×200 grids are divided into 18×18 big grids, so there is 11×11 original grids in one big grid. The densities are plotted by 3 dimensional contours. .............................................................................................. 46 4.5 Surface island distributions vs. temperatures (550K, 650K, 750K) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations. ..................... 47 4.6 Number of islands and average island size vs. temperature for islands growth on substrates with different orientations for short temperature range (550K to 750K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s................... 48 4.7 Number of islands and average island size vs. temperature for islands growth on substrate GaAs (001) for long temperature range (550K to 950K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s. ............................ 48 4.8 Surface island distributions vs. coverage c (10%, 20% 30% and 50%) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations. ..................... 49 4.9 Number of islands and average island size vs. coverage for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and ti=200s................................................................................ 50 xi

4.10 Surface island distributions vs. interruption times ti (0s, 10s, and 200s) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and c=20%, on a 200×200 grid. Red dashed lines indicate the island chain orientations.................................................................................................................. 51 4.11 Number of islands and average island size vs. interruption time for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and c=20%...................................................................... 52 4.12 Schematic illustration of superlattice........................................................................ 52 4.13 Atom island ordering on the first three layers. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Thickness of each layer is 5 grids. Strain energy field is calculated at every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for growth interruption. ........................................................................................ 53 4.14 Atom island ordering for different layer thickness. Figures (b) to (h) are the simulated results on the second layer with thickness equal to 5 grids to 30 grids. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Strain energy field is updated every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for interruption. ...... 53 5.1 Schematic illustration of 3D QD self-assembly model. An atom on top of the substrate surface (x,y)-plane in (a), the relative locations of the atom grid on the (x,y)-plane in (b), the corresponding binding energy Es related to the in-plane locations in (c)............................................................................................................. 64 5.2 Illustration of an edge atom A’s jump up or down during 3D QD self-assembly. ..... 64 5.3 3D island distributions for different up/down ratios ρ with total coverage c=1.6Ml. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and interruption time ti=200s (The total simulation time = the deposition time (td=160s) plus the interruption time (ti=200s)). .......................................................................... 65 5.4 Average island height vs. up/down ratio ρ (for ρ from 3 to 50). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and interruption time ti=200s. ..................................................................... 65

xii

5.5 Island distributions for different interruption times (ti=0s, corresponding to deposition time td=160s, ti=100s, and ti=200s) and for different up/down ratios (ρ=10, 20 and 30). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and total coverage c=1.6Ml. .................................................................. 66 5.6 Island distributions during and after deposition. Deposition time td=50s with 0.5ML coverage in (a), td=100s with 1ML coverage in (b), td=160s with 1.6Ml coverage (i.e., at the beginning of interruption time) in (c) and interruption time ti=200s in (d). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and up/down ratio ρ=10.................................. 67 5.7 Average island diameter (of the bottom equivalent circle of the island) vs. number of islands during and after deposition. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml and up/down ratio ρ=10..... 68 5.8 Island distributions for different deposition processes of 0.5Ml, 1Ml and1.6 Ml with flux rates F=1Ml/s, 0.1Ml/s and 0.01Ml/s. Fixed growth parameters are temperature T=700K, total coverage c=1.6Ml and up/down ratio ρ=10. ................... 69 6.1 Illustration of the geometry of a triplet of atoms used in the definition of the threebody potentials. ........................................................................................................... 86 6.2 Flowchart of molecular dynamics programming process........................................... 87 6.3 Illustration of cubic crystalline under uniform strain. ................................................ 88 6.4 Effective bulk modulus of InAs is plotted as a function of volumetric strain. Solid line is the results obtained by our program. Dashed line with triangles is the results from the ref. [Ell02]. Dashdotted line is calculated by commercial software. ..................................................................................................................................... 88 6.5 Effective elastic constant C11-C12 of InAs is plotted as a function of volumetric strain. Solid line is the results obtained by our program and dashed line is from the ref. [Ell02]............................................................................................................. 89 6.6 Effective elastic constants C11 and C12 of InAs are plotted as a function of volumetric strain. Solid line is the results obtained by our program. Dashed line with triangles is from the ref. [Ell02]. Dashdotted line is calculated by commercial software....................................................................................................................... 90 6.7 Effective elastic constant C44 InAs is plotted as a function of volumetric strain. Solid line is the results obtained by our program. Dashed line is from the ref. [Ell02]. Dashdotted line is calculated by commercial software. ............................... 91 xiii

6.8 Lattice energy (per atom) vs single silicon cubic lattice length.................................. 91 6.9 Lattice energy (216 atom system) vs system volume (V-2/3). ..................................... 92 6.10 Bulk modulus of Silicon vs volumetric strain........................................................... 92 6.11 C11-C12 of Silicon vs volumetric strain. .................................................................... 93 6.12 C11 of silicon vs volumetric strain. ........................................................................... 93 6.13 C12 of silicon vs volumetric strain. ........................................................................... 94 6.14 C44 of Silicon vs volumetric strain............................................................................ 94 6.15 Young’s modulus plotted as a function of volumetric strain.................................... 95 6.16 Poisson’s ratio plotted as a function of volumetric strain......................................... 95 7.1 Schematic arrangement of Epon 862 matrix filled with long (a) and short (b) CNTs. ........................................................................................................................ 106 7.2 Epon 862 molecular structure (a) and computer constructed molecules of Epon 862 (b)....................................................................................................................... 107 7.3 Computer constructed Epon 862 matrix with the size of l1×l2×l3= 40.28×40.28×61.09 Å3. (The CNT is embedded along the longitudinal x3direction as shown in Figure4).................................................................................. 107 7.4 Long (a) and short (b) CNT (10,10)-reinforced Epon 862 matrix (40.28×40.28×61.09 Å3). Lengths of the long and short CNTs are 61.09Å and 29.32Å, respectively. ................................................................................................ 108 7.5 Equilibrium van der Waals separation distance between the CNT (10,10) and Epon 862 matrix........................................................................................................ 108 7.6 MD simulated stress-strain curves for CNT (10,10), long and short CNT composites, and Epon 862 matrix............................................................................. 109 7.7 MD simulated stress-strain curves for CNT-reinforced Epon 862 composites and Epon 862 matrix........................................................................................................ 109 7.8 Stress-strain curves of the long CNT composite: MD simulated vs. Rule-ofmixture. ..................................................................................................................... 110 xiv

7.9 Variation of MD simulated Young’s moduli of CNT (10,10) and Epon 862 matrix with increasing strain. ............................................................................................... 110 7.10 Variation of MD simulated Young’s moduli of long and short CNT-reinforced Epon 862 composites and Epon 862 matrix with increasing strain. ......................... 111 B1 Extended periodic condition...................................................................................... 130 B2 Atom neighbors. ........................................................................................................ 130 B3 Flow chart 2D KMC QDs Growth Model ………………………………………….130

xv

CHAPTER I

INTRODUCTION

1.1 Atomistic Simulation in Nanostructured Materials

Nanostructured materials have greatly attracted people’s attention in the past ten years due to two main reasons. First, they are attractive from a scientific point of view. Nanostructured semiconductors provide a means to create artificial potentials for carriers, electrons, and holes at length scales comparable to or smaller than the de Broglie wavelength. And, nanostructured composites promise to dramatically increase material strength and durability [Sie96] as well. Second, quantum mechanics becomes applicable not only in systems of academic interest, but also to systems of practical impact. In particular, semiconductor nanostructures have a large potential for applications in nanoand optoelectronics. Over the past few years, computational materials modeling has become an increasingly important subfield of materials research, with its growing importance a result of the enormous increase in computing power and the maturation of a number of computational methodologies. Modeling has, in many cases, moved from being a tool for finding explanations of materials properties to being a means for prediction of the 1

properties of new classes of materials. It is a field that encompasses a large variety of techniques applied over a correspondingly disparate range of length and time scales, describing the structure and properties of wide classes of materials systems. The choice of a technique is, of course, dependent on the system under consideration and the phenomena that are being modeled. At feature sizes smaller than 100 nm, conventional theoretical modeling and computer simulations, such as those based on continuum approaches, must be supplemented with atomistic models. With conventional approaches, such as the finiteelement method based on linear elasticity, modelers cannot include the effects of defects, surfaces, and other microstructures – all of which fundamentally affect the material’s performance. Atomistic simulations, in particular those based on molecular dynamics (MD) and Monte Carlo (MC), capture these microstructural effects as well as highly nonhomogeneous stresses at the atomic scale. Simulations based on state-of-the-art molecular dynamics can involve up to 100 million atoms, but many researchers believe that petaflop computers will be built by the end of the next decade (one petaflop is 1015 floating-point operations per second). Such computers will make it possible to simulate nanostructures that involve 1012 atoms – with a total system size of well over one micron. This in turn means that modelers can simulate an entire device structure atomistically.

1.2 Multiresolution in Atomistic Simulation

Multisolution modeling and computation is a rapidly evolving area of research that has a fundamental impact on computational science and applied mathematics and 2

influences the way we view the relation between mathematics and science. Currently, multiresolution problems are driven mainly by the use of mathematical models in the applied sciences: in particular, material science, chemistry, fluid dynamics, and biology. Though, modeling at the level of a single scale, such as molecular dynamics or continuum theory, is becoming relatively mature. There is an urgent need from science and technology - nano-science being a good example - for multiresolution modeling techniques. It is not an exaggeration to say that almost all problems have multiple scales. We organize our time according to days, months, and years, reflecting the multiple time scales in the dynamics of the solar system. In addition, different physical laws may be required to describe the system at different scales. Take the example of fluids (figure 1.1) [E03]. At the macroscale (meters or millimeters), fluids are accurately described by the density, velocity, and temperature fields, which obey the continuum Navier-Stokes equations. On the scale of the mean free path, it is necessary to use kinetic theory (Boltzmann’s equation) to get a more detailed description in terms of the one-particle phase-space distribution function. At the nanometer scale, molecular dynamics in the form of Newton’s law has to be used to give the actual position and velocity of each individual atom that makes up the fluid. A general framework for multiresolution methods should ideally unify existing methods, give guidelines on how to design new methods, improve existing ones and provide a mathematical theory for stability and accuracy of these methods [E03]. There is a long history in mathematics for the study of multiresolution problems. Fourier analysis has long been used as a way of representing functions according to their components at 3

different scales [Bak99]. On the computational side, several important classes of numerical methods have been developed which address explicitly the multiresolution nature of the solutions. These include multigrid methods [Wes92], domain decomposition methods [Tos05], fast multipole methods [Gum05], and adaptive mesh refinement techniques [Vay04]. Many other mathematical techniques have been developed to study multiresolution problems, including boundary-layer analysis [Kev96], semiclassical methods [Mas81], geometric theory of diffractions [Kel62], stochastic mode elimination [Maj00], and renormalization group methods [Cho03].

1.3 Research Purpose

In our research, combining continuum theory, we focus on using MC and MD to study quantum dots (QDs) self-assembly and carbon nanotube (CNT) composites. The growth of self-organized QDs in strained semiconductors in the Stranski-Krastanov growth mode has been studied using kinetic Monte Carlo (KMC) simulations method. The set of relevant growth parameters such as growth temperature, surface coverage, flux rate, and growth interruption time are investigated. The dependence of the self-organized quantum dots pattern on elastic strain energy of the substrate is studied as well. In our research, by incorporating the anisotropic strain energy field into a KMC algorithm for adatom diffusion, we show that the self-organized island pattern on the surface of an anisotropic substrate is closely correlated to the elastic energy distribution on the surface. Because of the misfit strain during the QD growth, we present a molecular dynamic prediction for the elastic stiffnesses C11, C12 and C44 in strained silicon and InAs as 4

functions of the volumetric strain level. Our approach combines basic continuum mechanics with the classical molecular dynamic approach, supplemented with the Stillinger–Weber potential. Using our approach, the bulk modulus, effective elastic stiffnesses C11, C12 and C44 of the strained silicon, including also the effective Young’s modulus and Poisson’s ratio, are all calculated and presented in terms of figures and formulae. Except for QD self-assembly, using MD method, we study CNT composites as well. Single-walled carbon nanotubes (CNTs) are used to reinforce epoxy Epon 862 matrix. Three periodic systems --- a long CNT reinforced Epon 862 composite, a short CNT reinforced Epon 862 composite, and the Epon 862 matrix itself ---, are studied using molecular dynamics. The stress-strain relations and the elastic Young’s moduli along the longitudinal direction (parallel to the CNT axis) are simulated and the results are compared with those from the rule-of-mixture.

Figure 1.1 Different laws of physics are required to describe properties and processes of fluids at different scales [E03].

5

CHAPTER II

MONTE CARLO SIMULATION OF SELF-ORGANIZED QUANTUM DOT GROWTH

2.1 Quantum Dots

Quantum dots (QD), coherent inclusions in a semiconductor matrix with truly zerodimensional electronic properties, present the utmost challenge and point of culmination of semiconductor physics. It was at the beginning of the 1990s that a modified StranskiKrastanow growth mechanism driven by self-organization phenomena at the surface of strongly strained heterostructures was realized for the fabrication of such dots. And this process presents a sound way to fabricate easily and fast large densities of quantum dots [Bim99]. A semiconductor quantum dot is the ultimate quantum confined structure. Its unique electronic properties rely on the delta-function-like energy dependence of the density of states due to the quantum confinement of carriers in all three dimensions. However, to exploit the electronic properties in novel quantum effect devices, the lateral dimensions of the structures have to be in the range of, or smaller than, the de Broglie wavelength of electrons (50 nm in GaAs) [Bim99]. 6

Moreover,

for

practicable

applications in semiconductor devices, millions of the quantum structures, densely packed and uniform on the atomic scale, are necessary to achieve the desired active volume. This requires more precise fabrication methods to provide better control of size and shape of large ensembles of nanostructures. If this can be done successfully, it will allow the realization of novel high-speed, quantum interference, optoelectronic and single-electron devices. Ultralow threshold currents with high temperature stability have been predicted for quantum dot lasers [Ara82] and high mobilities at room temperature in quantum-dot superlattices with appropriate choice of electronic coupling to suppress scattering by optical phonons. Advanced crystal growth techniques such as molecular beam epitaxy (MBE) and metalorganic vapour phase epitaxy (MOVPE) have made it possible to precisely fabricate two-dimensional layered semiconductors, i.e. heterojunctions, quantum wells and superlattices, on an atomic scale. Further reduction of the dimensionality of semiconductors in one-dimensional quantum wires and zerodimensional quantum dots, however, has turned out to be very difficult. In the past decade, the fabrication methods for nanometer-scale structures have been continuously improved and many fundamental aspects of low-dimensional semiconductors have been evaluated

[Kas90].

However,

the

reproducible

fabrication

of

device-quality

nanostructures is still a challenge, making the applicability of quantum wires and quantum dots somewhat haphazard. Originally, the fabrication methods of quantum wires and quantum dots were based on the lateral patterning of two-dimensional heterostructures by combining fineline lithography with wet and dry chemical etching. 7

Unfortunately,

however,

lithographic patterning always introduces irregularities in the shape of the nanostructures, and mechanical damage to the interfaces cannot be avoided. Therefore, the most promising quantum structures so far have been fabricated using techniques based on direct crystal growth. Along with the selective growth on prepatterned and masked substrates, the direct synthesis of nanostructures during the epitaxial growth process itself has become very important due to its potential for creating damage-free structures [Not93].

2.2 Self-organized Growth

Self-organized growth is a promising way to grow QD islands with uniform size and spatial distribution. However, it is still very challenging to self-organizationally grow high density and dislocation free quantum dot layers, even under the Stranski-Krastanov growth model [Dar97, Shc99b, Spr00]. Since experimentally this trial-and-error growth approach could be quite expensive and time-consuming, numerical simulation methods can be utilized to provide important information and guidance to the experimentation. As the self-organized QD growth is a competitive balance process between the binding and thermal energies, involving atom deposition and diffusion, certain kinetic Monte Carlo (KMC) models have been proposed to numerically simulate the self-organized growth [Mei03]. It was observed by Schöll and Bose [Sch98] and Elsholz et al. [Els03a] that a probability governed by the Arrhenius law including the sole atomic binding energy would result in the Ostwald ripening phenomenon in the KMC model. However, when the long-range strain energy is considered, cooperative growth where larger islands 8

lose some atoms to smaller ones could be achieved [Spr00].

2.2.1 Deposition and Diffusion The first step towards the formation of quantum dots is the deposition of a certain material on a given surface. Molecular beam epitaxy (MBE) belongs to the physical vapor deposition (PVD) methods. The substance to be deposited is vaporized, and the molecules formed into a beam. This beam is under ultra high vacuum conditions directed towards a surface, where the particles condensate. To allow for diffusive motion of the deposited atoms, the sample is heated. Since no other chemically active substances are used in this form of epitaxy, deposition and diffusion are characterized by only the most basic processes. For this reason the growth technique of MBE is very convenient for theoretical modeling of growth. On the other hand in chemical vapor deposition (CVD) the growth material is brought to the sample in form of a chemical carrier gas solution. A prominent example is metal organic CVD (MOCVD). Deposition occurs via chemical reactions at the sample surface. These reactions can be very complex and since the deposition rate depends sensitively on the concentrations of the various chemically active species the local growth kinetics can become quite involved. Furthermore, reactions usually occur both ways and consequently desorption is much more important in CVD than in MBE. Another widely used epitaxial method is the liquid phase epitaxy (LPE), where the substrate is submerged into an oversaturated liquid solution. Material precipitates from the solution and is deposited on the sample. Here, again, the relevant deposition processes are basic but diffusion is widely aided by the presence of a liquid phase. 9

While in MBE and low pressure CVD the interaction of the growing surface with the ambient atmosphere is of minor importance, the theory of growth from dense phases like LPE or growth from the melt can become most complex. Here, for instance, a simultaneous solution of the Navier-Stokes equation and mass- and energy-transport equations is required. Obviously, the technical details in various epitaxial setups can differ strongly. Nevertheless all methods used for epitaxy follow a simple scheme that consists of transporting material to the sample surface, depositing it and allowing for diffusion. The deposition might be as simple as scattering atoms on the surface (as in MBE or related sputtering techniques) or complex and dominated by chemical reactions (as in CVD) but ultimately any of the above mentioned techniques ensures a certain flux of particles to the surface. The flux is the relevant quantity in terms of growth kinetics and its physical effect is basically independent of the applied method of deposition. Once an atom is deposited on the surface it can travel from A to B on various paths (figure 2.1). Plain surface diffusion (figure 2.1b) is an important mechanism in all growth techniques and consists of consecutive hops from one lattice site to a neighboring one. Since it is the energetically most favorable way of diffusion, at least for short distances, it generally dominates all other ways of transportation.

2.2.2 Growth Classification If material A is deposited on material B it is not at all clear in which way the growth will occur. Additionally, for a given material system the mode of growth depends on external parameters like temperature and pressure. The simplest growth mode of 10

hetereogrowth, where one complete monolayer grows after the other is rather the exception than the rule. Under certain conditions rare gases grow layer-by-layer on graphite. Another example is the AlxGa1-x film growth on sapphire (001) as reported in [Wic94]. This growth mode is also called Frank-Van der Merwe [Fra49] growth mode (Figure 2.2b). Conversely, in growth experiments with lead on graphite the lead does not form monolayers. It rather forms little droplets very similar to water droplets on a freshly sealed car top. This growth mode is referred to as the Volmer-Weber [Vol26] growth mode (figure 2.2a). Indeed, this analogy between lead and water droplets is worth pursuing. By defining the surface tensions for a liquid droplet on a plane surface as σsv, σsl and σlv for the interfaces solid/vapor, solid/liquid and liquid/vapor, respectively, one can show that the droplet forms an angle θ with the solid given as:

cosθ =

σ sv − σ sl σ lv

(2.1)

For σsv Ecorner>Eedge.

2.3.2.2 Incorporation of the strain energy field As we have pointed out in the introduction, while the strain energy field is very important in the self-organized growth, it needs to be accurately calculated. In other words, for the strained semiconductor structure, one should use the strain energy corresponding to the point-eigenstrain, instead of that to the point-/line-force. We briefly present the correct formulation for the calculation of the strain energy below, leaving the details in Appendix A for the sake of reference. Within the framework of continuum elasticity theory, the elastic strain field induced by the islands can be obtained using the Green’s function solution for the anisotropic semiconductor substrate developed by refs [Pan02a, Pan02b]. Assume that there is a point misfit strain (or eigenstrain) γ ij* (due to the misfit lattice difference) at x, the induced elastic strain at y can be found as

1 2

k p * γ kp ( y ) = γ lm [σ ml , p ( x; y ) + σ ml , k ( x; y )] y

y

(2.6)

where the superscript k or p (=1,2,3) on the right side of equation (2.6) indicates the direction of the point-force, and the subscript prime followed by py or ky in the stress component σml denotes the derivative of the Green’s stress components. With the elastic strain equation (2.6), the strain energy at y due to an area A of island (with unit thickness in the out-of-plane z-direction) centered at x can be calculated as E str =

1 C ijkl γ ij γ kl A 2

(2.7)

18

where Cijkl is the elastic moduli. In the following calculation, we assume that the misfitstrain is hydrostatic, i.e., γ Ij* = γ *δ ij with γ*=0.07. The non-zero elastic coefficients for GaAs (001) [Pan02b] are taken to be C11=118, C12=54, and C44=59 (109N/m2).

2.3.3 Kinetic Monte Carlo Algorithm The simulation routine is based on a continuous time Monte Carlo scheme. The BKL algorithm [Bor75] named after Bortz, Kalos, and Lebowitz is very efficient for the problem at hand, since the independence of a particular time scale is of great advantage in simulating surface diffusion. In our program, atom diffusion processes are simulated one by one. The principal course of a simulated diffusion process is always the same and is sketched in the following. Each atom at most has four possible nearest and four possible next-nearest diffusion positions (figure 2.5), or four nearest and four next-nearest neighbors. Every possible diffusion step of a given atom is done by evaluating the probability p from equation (2.2). These eight possible diffusion probabilities to nearest and next nearest neighbor positions are stored and added to get the total diffusion probability for the atom, named pAtom. The total probability that anything might happen at all is then given by adding up all the pAtom to get pTotal. During the simulation, an atom is randomly chosen to move across the surface by hopping to nearest or next nearest lattice site. The key is to assign the proper weight to each atom to be chosen, considering that different atom contributes differently to the overall probability pTotal (see, e.g. equations. (2.2) & (2.5)). From the eight possible diffusion processes of the selected atom, one is chosen in accordance to its likelihood and executed. The corresponding time interval ∆t is 19

calculated and added to the elapsed simulation time. The movement of an atom also alters the diffusion barriers for the neighboring atoms in the previous neighborhood as well as in the new one. As such, the moving atom and all the atoms in its surrounding have to be recalculated to obtain the new diffusion probabilities. Strictly speaking, the strain energy field would have to be calculated at each step. However, since the strain evaluation is a lengthy procedure and the strain only slightly changes with the motion of a single atom, it turns out that it is sufficient to calculate the strain field, say, every 2500 jump steps. To speed up further the computation, the strain energy calculation does not extend over the whole system but only over a circular area of a given radius r around the point where the strain is to be evaluated, due to the rapid decay of this field. In the following examples, such a radius r is taken to be 30 lattice grids [Sch98]. The detailed simulation steps are listed in Appendix B.

20

B

a b A c

Figure 2.1 Possible ways of diffusion from A to B: a - desorption, condensation process; b -surface diffusion; c - volume diffusion

a)

b)

c)

Figure 2.2 Important growth modes in epitaxy: a) Volmer-Weber growth mode, b) FrankVan der Merwe growth mode, c) Stranski-Krastanov growth mode.

21

a)

b)

Figure 2.3 Strain relief by a) generation of mis-fit dislocations (green atoms) and b) Stranski- Krastanov growth mode. The green line marks the end of the wetting layer. Atoms of the substrate are blue, deposited material is red.

Eoff Eedge Ecorner

y

x Figure 2.4 Schematic illustration of off, edge, and corner diffusions.

y

x Figure 2.5 Schematic illustration of nearest (horizontal and vertical directions in solid blue) and next nearest neighbor (diagonal directions in dashed brown) positions of each atom. 22

CHAPTER III

GROWTH PARAMETERS DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION

3.1 Growth Parameters

For most electronic and optical applications that include active quantum dot layers, it is important to have as many and as equally shaped dots as possible. Unfortunately, the problem of generating self-organized QD patterns cannot be solved in a straightforward manner. Even though Stranski-Krastanov growth is known to be suitable for selforganized QD growth, there are many parameters that can considerably influence the growth result. Some of the relevant parameters are temperature T, flux rate F, coverage c, and interruption time ti. In what follows, we discuss the effect of these parameters on the island growth of material GaAs (001) on a 200×200 lattice size with periodic condition and with strain energy for all atoms being updated every 2500 hopping steps.

23

3.2 Temperature T

The growth temperature T affects the simulation result greatly. If T is chosen too low, the deposited atoms will just stick to the surface without having enough thermal energy to diffuse. In this case, of course, no self-organization is to be expected. If, on the other hand, the temperature is too high, then interatomic bonds are too easily overcome and one observes an ensemble of monomers and small polymers of atoms performing random walks over the surface. Under this situation, large islands are inherently unstable. Therefore, only for a distinct interval of temperatures self-organization is effective and accessible to production processes. The temperature dependence of the average island size is shown in figure 3.1 for fixed values of flux F=1Ml/s and coverage c=20%. Here, deposition stops after 0.2s of simulation time and the growth interruption time ti is 200s. It is observed from figure 3.1 that the higher the temperature, the larger the average island size, and that when the temperature is higher than 700K, e.g. from 700K to 800K, there is no big difference among the ultimate island size; however increasing temperature in this interval, we obtain more and more ordered island ordering.

3.3 Flux Rate F

The flux rate F is a measure of how fast the material is deposited on the surface. It does not include any equation concerning processes connected to self-organization, like the definition of diffusion barriers or the strain energy field. Thus, the flux rate 24

influences the kinetics of growth only indirectly by determining the density of atoms on the surface and the nucleation rate of islands. Since flux is only present during the time of deposition, one might assume that it is even less important during growth interruption. However, different values of flux can lead to significant differences in the surface morphology after the end of deposition. Under low flux rate, the system will have enough time during the deposition process to come close to an equilibrium size distribution with large-size islands. On the other hand, a small-size island distribution is associated with a high flux rate. In other words, in general, the lower the flux rate, the larger the equilibrium island size, as is also illustrated in figure 3.2 for different flux rates but with fixed temperature T=700K, coverage c=20%, and interruption time ti=200s. We point out that to highlight the effect of flux rate, strain energy is not included in this simulation.

3.4 Surface Coverage c

The surface coverage c describes how much material is deposited on the surface. This parameter certainly has an effect on size distributions, because one cannot expect islands of equilibrium size if there is not enough deposited material. On the other hand, if too much material is transported to the surface than necessary for optimal island sizes, too small distance between islands makes the formation of isolated islands impossible. For coverage below certain critical coverage cc, the growth is mainly kinetically controlled, so the atom size distribution is mainly controlled by a variation of coverage and one should expect an increase of the average island size with increasing coverage. Under this condition, every island collects atoms from its 25

immediate

vicinity

without

any

considerable exchange of material among islands, and its size is directly determined by the amount of the deposited material around the island. For temperature T=700K, flux F=1Ml/s, and interruption time ti=100s, the simulated results on the 200×200 grid with different coverage are plotted in figure 3.3. From the results, we can find that with coverage increase from 10% to 30%, the average island size becomes larger and larger. However, when the coverage is larger than 30%, isolated islands cannot be easily observed. In other words, the critical coverage is between 20% and 30%.

3.5 Growth Interruption Time ti

Another important factor to influence the growth result is the time between the end of deposition and the capping of the QD layer with another material, the so-called growth interruption. The time of this scale is called interruption time ti, with the time from the very beginning to the end of simulation is named total simulation time t (= ti + deposit time). It is known that the growth interruption affects the crystal surface dramatically and the results before and after interruption are different [Kam96]. Furthermore, during the growth interruption, the atoms can move to energetically favorable positions and approach thermodynamic equilibrium. Since there is a striking difference between kinetically controlled growth and a thermodynamically dominated size distribution, the effect of a growth interruption can be dramatic. With the increase of interruption time, the system exhibits ordered island size and ordering. For temperature T=750K, flux F=1Ml/s, and coverage c=20% on the grid of 26

200×200, the simulation results with different interruption times are shown in figure 3.4(a)-(f). It is noticed clearly that with increasing ti, the island size increases; however, when ti=100s, the size reaches the equilibrium, and we finally obtain perfect island ordering when ti=200s. In other words, with increasing interruption time, the system ultimately reaches the equilibrium.

(a) T=550K

(d) T=700K

(b) T=600K

(c) T=650K

(e) T=750K

(f) T=800K

Figure 3.1 Atom island ordering under different temperatures (T is equal to 550K in (a), 600K in (b), 650K in (c), 700K in (d), 750K in (e), and 800K in (f)). Simulation parameters are F=1.0Ml/s, c=20% and ti=200s on a 200×200 grid. The strain field is updated every 2500 steps.

27

(a) F=1Ml/s

(b) F=0.1Ml/s

(c) F=0.01Ml/s

Figure 3.2 Island ordering of atoms under different flux rates. Simulation parameters are T=750K, c=20% and ti=200s on a 200×200 grid. Strain field is NOT included. Deposition stops after 0.2s in (a), 2s in (b), and 20s in (c).

(a) c=10%

(b) c=20%

(d) c=40%

(c) c=30%

(e) c=50%

Figure 3.3 Atom island ordering with coverage c increasing from 10% to 50% in (a) to (e). Simulation parameters are T=700K, F=1.0Ml/s, ti=200s and on a 200×200 grid. The strain energy field is updated every 2500 steps. 28

(a) ti=0

(d) ti=100s

(b) ti=10s

(c) ti=50s

(e) ti=150s

(f) ti=200s

Figure 3.4 Atom island ordering for different interruption times. The ordering of (a) is the initial atom pattern when all the atoms just deposit to the surface and the interruption time ti=0. Simulation parameters are T=750K, F=1.0Ml/s and c=20% on a 200×200 grid. The strain energy is updated every 2500 steps.

29

CHAPTER IV

ELASTIC ANISOTROPY DEPENDENCE OF QUANTUM DOTS TWO DIMENSIONAL DISTRIBUTION

4.1 Strain Field in Quantum Dot Growth

In self-organized QD research, both theoretical [Ter96, Liu99] and experimental [Tei96, Dar97, Hol98] works have shown that the interaction among dots via their elastic strain field [Xie95, Sol96] may lead to a gradual improvement in size homogeneity, as well as to a more uniform lateral island spacing. Holy and coworkers have further shown that the vertical and lateral correlations in self-organized QD superlattices can be explained by taking into account the elastic anisotropy of the material [Hol99]. In these superlattices, dot correlations are induced by the interaction among dots via their elastic strain fields. Above the buried islands, the elastic energy distribution on the surface exhibits pronounced minima and maxima in the lateral directions. The continuum finite element analysis and many other numerical methods were proposed to study the effect of elastic anisotropy on the three dimensional (3D) strain field, as well as on the layer-tolayer dot correlations to the growth direction [Lee04, Que03, Spr02]. Although QDs ordering due to the elastic anisotropy has been studied [Spr02], the relationship 30

between elastic anisotropy and QD array lateral distribution, the very important character of QD self-organization, has not been investigated. Yet, depending on the elastic properties and growth orientation, a variety of lateral island ordering can be found, which is potentially correlated to the high elastic anisotropy [Sri91, Hol99]. The kinetic Monte Carlo (KMC) has been proposed to simulate QD self-organization by many researchers [Ter96, Joy04, Sri93, Mei03, Tha04]. It has been shown that the main parameters affecting QD island distribution pattern are the temperature T, flux rate F to the surface during deposition, surface coverage c, and growth interruption time ti [Mei03, Pan05]. The simulated island ordering and narrow size distribution in two dimensions were discussed and optimized growth parameters were identified in the former chapter. However, how the substrate anisotropy affects the island pattern has not been studied so far within the KMC simulation algorithm. Therefore, in this chapter, we use a kinetic Monte Carlo (KMC) self-organization model to show how the elastic anisotropy affects the QDs lateral ordering. In our simulation, the atom hopping rate is governed by the Arrenius law modified with the long-range strain energy field. The elastic strain field induced by the islands is obtained using the Green’s function solution for the anisotropic semiconductor substrate developed previously [Pan02a, Pan02b]. QDs array patterns without elastic strain energy, with elastic isotropic strain energy, and with anisotropic strain energy in three different growth orientations, i.e., GaAs (001), GaAs (111) and GaAs (113), are studied. Simulation examples with different growth parameters involving temperature T, coverage c, and interruption time ti are discussed, and our results show that without strain energy, isolated but random QDs nucleation can be

observed. With consideration of the strain 31

energy for isotropic and anisotropic cases, ordered and different island chains can be observed, which proves that strain energy affects the island ordering greatly. Results with different growth parameters show that temperature, coverage and interruption time affect the island ordering quality, but none of them affects island ordering tendency. Thus, the elastic strain energy seems to be the only factor that controls the self-organized island pattern. Furthermore, our simulation results also qualitatively agree with recent experimental results obtained via scanning tunneling microscopy (STM) [Bru95, Bru98].

4.2 Anisotropic Strain Energy Field Calculation

As we mentioned in chapter 2.3.2.2, the elastic strain field induced by the islands can be obtained using the Green’s function solution for the anisotropic semiconductor substrate developed previously and the non-zero elastic coefficients for GaAs (001) are calculated by taking C11=118.8×109N/m2, C12=53.8×109N/m2, and C44=59.4×109N/m2 [Pan02b]. Here, an isotropic model, called Iso (001) is also included for the purpose of comparison in which the elastic constants of GaAs (001) are replaced with its isotropic constants. In the isotropic Iso (001) model, the elastic constants C12 and C44 are the same as those of GaAs (001), i.e., C12=53.8×109N/m2 and C44=59.4×109N/m2, with C11=172.6×109N/m2 obtained by imposing the isotropic condition C11-C12-2C44=0 [Pan03]. 4.2.1 Material Properties of Iso(001), GaAs(001), GaAs(111) and GaAs(113) The elastic moduli of GaAs (111) and GaAs (113) are transformed from GaAs (001) 32

[Pan02a]. The QDs grow along z axis of substrate, which means, for 2D growth, QDs self-organize on xy plane. For GaAs (111) with the substrate coordinates of x axis along (11-2), the y axis along (-110) and the z axis along the (111), for GaAs (113) with the substrate coordinates of x axis along (33-2), the y axis along (-110) and the z axis along the (113), the global material properties of GaAs (111) and GaAs (113) are transformed by the well-known tensor transformation [Nye85] from GaAs (001). Each transformation is done by two steps. For GaAs (111), first, anticlockwise rotate xy surface along z axis

π 4 (from positive direction of z axis), then, anticlockwise rotate xz surface along y axis

(

π 2 − cos −1 2 3

)

(from positive direction of y axis). For GaAs (113), first,

anticlockwise rotate xy surface along z axis π 4 (from positive direction of z axis), then,

(

anticlockwise rotate xz surface along y axis tan −1 3

)

2 (from positive direction of y

axis). Figure 4.1 shows the coordinates of GaAs(001), (111) and (113) and the elastic constant matrixes of Iso (001), GaAs (001), (111) and (113) are shown in equations (4.14.4). Material properties of Iso(001): 0 0 0 ⎤ ⎡172.6 53.8 53.8 ⎢ 53.8 172.6 53.8 0 0 0 ⎥⎥ ⎢ ⎢ 53.8 53.8 172.6 0 0 0 ⎥ C = 10 9 × ⎢ ⎥ Pa 0 0 59.4 0 0 ⎥ ⎢ 0 ⎢ 0 0 0 0 59.4 0 ⎥ ⎢ ⎥ 0 0 0 0 59.4⎥⎦ ⎢⎣ 0 Material properties of GaAs (001):

33

(4.1)

0 0 0 ⎤ ⎡118.8 53.8 53.8 ⎢ 53.8 118.8 53.8 0 0 0 ⎥⎥ ⎢ ⎢ 53.8 53.8 118.8 0 0 0 ⎥ C = 10 9 × ⎢ ⎥ Pa 0 0 59.4 0 0 ⎥ ⎢ 0 ⎢ 0 0 0 0 59.4 0 ⎥ ⎢ ⎥ 0 0 0 0 59.4⎥⎦ ⎢⎣ 0

(4.2)

Material properties of GaAs (111): 45 36 0 12.73 0 ⎤ ⎡ 145 ⎢ 45 − 12.73 145 36 0 0 ⎥⎥ ⎢ ⎢ 36 36 154 0 0 0 ⎥ C = 10 9 × ⎢ ⎥ Pa − 12.73⎥ 0 0 41 0 ⎢ 0 ⎢12.73 − 12.73 0 0 41 0 ⎥ ⎢ ⎥ 0 0 − 12.73 0 50 ⎦⎥ ⎣⎢ 0

(4.3)

Material properties of GaAs (113): − 4.72 0 0 ⎤ ⎡152.81 31.79 41.79 ⎢ 31.79 145.7 − 10.38 48.91 0 0 ⎥⎥ ⎢ ⎢ 41.79 48.91 135.70 0 15.09 0 ⎥ C = 10 9 × ⎢ ⎥ Pa − 0 0 0 54 . 51 0 10 . 38 ⎢ ⎥ ⎢ − 4.72 − 10.38 15.09 0 47.39 0 ⎥ ⎢ ⎥ − 10.38 0 0 0 37.39 ⎦⎥ ⎣⎢ 0

(4.4)

4.2.2 Strain Energy Field of Iso(001), GaAs(001), GaAs(111) and GaAs(113) For the strain energy field calculation, due to its rapid decay behavior, the strain energy field is not evaluated over the whole system but only over a circular area of a given radius around the field point, which is chosen to be thirty grids [Pan05]. This decay feature can be also observed from figure 4.2 where the elastic strain energy field on the surface of substrates Iso (001), GaAs (001), GaAs (111) and GaAs (113) due to a buried point QD are plotted. On the other hand, figure 4.2 also shows clearly the direct 34

correlation between strain energy contour shapes and the growth orientations. The different strain energy patterns on the surface of Iso (001), GaAs (001), GaAs (111) and GaAs (113) will result in island distributions with different inclined island chain orientations. 4.3 Effects of Elastic Strain Energy on Island Ordering

As we have discussed in 2.3.2, during the formation of islands on the substrate, an elastic strain field is also generated around the islands. The nearest and next nearest neighbor binding energies contribute to the shape of the compact islands; however, the island-to-island interaction is very weak without strain energy. As the strain energy is calculated based on the point-eigenstrain Green’s function, it is expected that the absolute value of strain energy is the largest at the boundary area of an island, and with increasing distance away from the island it decays rapidly [Pan02a, Pan02b, Dar97]. Thus, the interaction among islands could lead to island lateral and vertical ordering due to the strain field. This is why that, without elastic strain energy, neither uniform size nor ordered spatial ordering can be observed and one would expect Ostwald ripening [Els03, Zhd99]. To demonstrate the importance of the strain energy and the effect of elastic anisotropy on QD arrays, we take the semiconductor material GaAs as an example. Island ordering of GaAs without strain energy is shown in figure 4.3e, and with strain energy on the surface of substrates Iso (001), GaAs (001), GaAs (111), and GaAs (113), respectively, in figure 4.3(a), (b), (c), and (d). The simulation parameters are temperature T=750K, flux rate F=1.0Ml/s, surface coverage c=20%, and interruption time ti=200s. To 35

show the island distribution more clearly, the island density contours are plotted in figure 4.4 by dividing the original 200×200 small grid size into 18×18 big grid size. So the big grid centers are used as the new coordinates and adatom number of each grid is used as the height. From figure 4.3, it can be observed clearly that without strain energy, isolated but randomly ordered islands can be formed (figure 4.3e) often associated with Ostwald ripening [Zhd99]; however, with the strain energy, isolated and ordered islands (for both size and distribution) can be formed (figure 4.3a,b,c,d). Furthermore, for different growth orientations GaAs (001), GaAs (111) and GaAs (113), different island arrays can be observed (figure 4.3b,c,d), which correspond to the different strain energy patterns shown in figure 4.2 due to material anisotropy with different growth orientations. These results indicate that the long-range strain energy could be one of the very important factors controlling ordered QD arrays and ordered island size. Actually, images of island inclinations on different anisotropic substrates from experiments have also shown the direct correlation between the QD island pattern and the substrate anisotropy (or growth orientation), as by Jacobi [Jac03] where InAs QDs were grown on GaAs (001) and GaAs (113) substrates and by Brune et al. [Bru95, Bur98] where the influence of elastic strain on Ag self-diffusion and nucleation was investigated experimentally. Besides the long-range strain energy field, other key growth parameters include the temperature, coverage, and interruption time, as have been discussed in the previous chapter. However, once we have selected the optimal parameters for temperature, coverage, and interruption time, the island ordering and QD pattern will be mostly controlled by the growth orientation, namely by the anisotropic elastic energy in 36

the substrate with material anisotropy. Therefore, in what follows, we will study these growth parameters by using different growth orientations for the semiconductor GaAs.

4.4 Effect of Substrate Anisotropy and Growth Parameters on Island Ordering

Numerous previous studies have shown that without strain energy field, one could obtain isolated but randomly ordered islands, which are caused by the nearest and next nearest binding energies among adatoms and substrate atoms; however, when the strain energy field is included, isolated and ordered island array (both in size and distribution) can be achieved. Brief results presented in the previous section have further shown that isotropic and anisotropic substrates with different growth orientations could result in different island chain orientations. Furthermore, in the previous chapter, we studied in detail the optimized growth parameters (T, F, c, and ti) for the substrate GaAs (001) and our simulation indicated that these growth parameters could greatly affect the island size and ordering. Here we would like to further study, for the anisotropic substrate case, how these parameters affect the island chain orientation. In other words, is anisotropic elastic strain energy the key parameter in controlling the island chain orientation? To answer this critical question, we have selected four different substrates in our simulation: isotropic substrate Iso (001) and three GaAs substrates but with different growth directions, namely, GaAs (001), GaAs (111) and GaAs (113).

4.4.1 Temperature T & Substrate Anisotropy Simulation

results

for

different

temperatures are shown in figure 4.5 with 37

fixed flux rate F=1Ml/s, coverage c=20%, and interruption time ti = 200s. We selected three different temperatures in the simulation, i.e., T=550K, T=650K, and T=750K, with the results on the substrate isotropic Iso (001) in figure 4.5a and on the anisotropic GaAs (001) in figure 4.5b, GaAs (111) in figure 4.5c and GaAs (113) in figure 4.5d. From figure 4.5, we observed that for all the substrate cases with low temperature (T=550K), no clear pattern can be formed. This is due to the fact that most adatoms do not have enough thermal energy to diffuse, and therefore, the interaction between islands is very weak and no island chains can be formed. With increasing temperature, the island size increases and the island chains become clear, which is also in agreement with experimental results [Sai99]; however, if the temperature is still below the critical value (around 750K as has been identified [Pan05]), the island distributions could be still irregular with no clear correlation to the growth directions (figure 4.5, the middle vertical column for T=650K). On the other hand, if the temperature is too high, the high thermal energy will make the big island unstable, as was discussed before [Pan05]. Therefore, only around the optimized temperature (i.e., T=750K), can one clearly observe ordering and regular islands distribution. More importantly, different island ordering patterns are closely correlated with different growth orientations. Therefore, temperature could affect the island size and island ordering, but it can not control the island chain orientation; the island pattern is directly dominated by the substrate anisotropy or the growth orientation once the optimal temperature is selected. Figure 4.6 and figure 4.7 show the variation of number of islands and average island size vs. different temperatures. Other fixed growth parameters are: flux rate F=1Ml/s, coverage c=20%, and interruption time ti =

200s. We observed from figure 4.6 that, 38

under the given growth conditions, the number of islands and average island size are insensitive to the substrate anisotropy (or growth orientation). For all the substrate, with increasing temperature to 750K, the island size increases and the number of islands decrease. In our previous discussion, we further identified by looking at the island size ordering and distribution that the optimal temperature was around 750K. The island distributions could be still irregular if temperature is too low (below this critical value) and be unstable if temperature is too high (above this critical value). Therefore, only around the optimized temperature (i.e., T=750K), can one clearly observe ordered and regular islands distribution. To further confirm our previous observation, we show here also the variation of the average island size and number of islands versus temperature for a large temperature range (figure 4.6 for GaAs (001) only). It is observed from figure 4.7 that the average size of the island reaches a maximum at T=750K where the number of islands reaches a minimum.

4.4.2 Surface Coverage c & Substrate Anisotropy Simulation results of island formation on substrates Iso (001), GaAs (001), GaAs (111) and GaAs (113) with fixed parameters of temperature T=750K, flux rate F=1Ml/s and interruption time ti = 200s are plotted in figure 4.8 for different coverage c=10%, c=20%, c=30%, and c=50%. It is observed from figure 4.8 that, if the coverage is relatively low, for example for coverage from c=10% to c=20%, the island size doesn’t increase much with increasing coverage, which is in agreement with the experiment results by Leonard et al [Leo94]. On the other hand, if the coverage is too high, for example at c=50%, chain-like islands, instead of isolated islands, are formed. The 39

most important feature is that, whether they are isolated or continuous chain-like islands, different chain orientations can be observed and are clearly correlated to the different growth directions. Consequently, coverage can only control whether the island is isolated or chain-like, whilst the substrate anisotropy controls the orientation of the island array or island-chain. Figure 4.9 shows the variation of the number of island and average island size on the surfaces of the substrates Iso (001), GaAs (001), GaAs (111) and GaAs (113) for different coverage c. The other fixed parameters are temperature T=750K, flux rate F=1Ml/s and interruption time ti = 200s. It is interesting that with increasing coverage, while the number of islands in the grid is nearly a constant value (around 50), the average island size increases, with GaAs (113) having the largest value followed by GaAs (111). The average island size on the surface of the substrate GaAs(001) has the smallest value.

4.4.3 Growth Interruption Time ti & Substrate Anisotropy We now investigate the effect of the growth interruption time ti on the island formation on the surface of four different substrates Iso (001), GaAs (001), GaAs (111), and GaAs (113) with fixed parameters T=750K, F=1.0Ml/s, and c=20% (figure 4.10). The growth interruption time ti is taken as 0s, 10s, and 200s. The results in figure 4.10 clearly show that the strain energy affects the QD arrays and the substrate growth orientation controls the island ordering pattern even from the earlier interruption time say at 10s. It is also observed that, with increasing interruption time, the island inclinations become clearer and clearer, and finally come to the final pattern for ti=200s (figure 4.10). Therefore, suitable amount of interruption time is required for the atom islands to 40

form ordered size and spatial distribution. Obviously, however, long interruption time can only make island size/distribution regular, while the island-chain inclination or island pattern is controlled by the substrate anisotropy or growth orientation. Figure 4.11 shows the effect of the growth interruption time ti on the number of islands and average island size. The fixed growth parameters are T=750K, F=1.0Ml/s, and c=20%. It is observed clearly from figure 4.11 that, with increasing interruption time, the number of islands decreases whilst the average island size increases. Furthermore, at large interruption time (ti=200s for our case), while the number of islands is nearly the same for different substrate orientations (or different material anisotropies), different substrate orientations predict slightly different island sizes with GaAs (111) having the smallest size among all the substrates.

4.5 Stacks of Quantum Dots

Having shown the effect of temperature, flux rate, surface coverage, and growthinterruption time on the island formation (where only a single layer over the substrate is concerned), we now study the interaction between atom layers. As before, the important contribution of the strain energy field will be emphasized. To increase the uniformity and ordering of the QD density over the surface, various approaches have been proposed, including growth of QDs over patterned substrates. Strained-layer heteroepitaxy has been studied in lateral as well as vertical direction [Moi94, Leo94, Not96]. Recently, it has been demonstrated that regular threedimensional QD superlattices with tunable lattice constants can be achieved utilizing 41

strain-mediated self-organization effects in Stranski-Krastanov growth mode [Spr00, Lee04]. There is evidence that stacks not only increase the overall QD density, but also improve the size distribution and the spatial ordering [Dar99]. In this section, we demonstrate the spatial correlation in the vertical direction by carrying out the simulation for a three-layer model with temperature T=700K, flux rate F=0.01Ml/s, coverage c=20%, interruption time ti=50s, and thickness for each layer H=5 grids (figure 4.12). The simulated results for the first three QD layers are shown in figure 4.13. In this model, we first deposit an initial atom pattern similar to figure 4.13a on the substrate GaAs (001). Then our KMC is executed and the simulated pattern for the first layer is shown in figure 4.13a. Next, on top of the first layer, we randomly deposit an initial pattern similar to that of the first layer with distance between the two layers being 5 grids. Simulation similar to the single layer case is then carried out. However, when we calculate the strain energy field at every 2500 hopping steps on the new surface, we consider not only the strain contribution from the surface as we did for figure 4.13a but also that from the substrate (first layer) due to the previous formed atomic pattern. The simulation again runs to the interruption time ti=50s and the result is plotted in figure 4.13b. Finally, following the same procedure, the third layer is deposited over the second layer, and the strain energy at each atom location in the third layer is evaluated by considering all the contributions of atoms in layer 1 to layer 3, with result being shown in figure 4.13c. From figure 4.13, we observed that, with increasing number of deposited layers the island size tends to become larger and the corresponding spatial distribution become more uniform [Spr00, Spr01, Ter96]. This phenomenon is understandable as the 42

additional strain energy from the buried layers enhances the mobility of the diffusing atoms and at the same time, controls the distance at which the atoms can hop. Another interesting factor that affects QD distribution is the thickness of the buffer (or spacer) layer. We have run our KMC program for different layer thickness, with results in figure 4.14b to figure 4.14f corresponding to layer thickness H equal to 5, 10, 15, 20, and 30 grids. It is observed that a thinner layer corresponds to a better island distribution [Ter96] since the thicker the layer, the smaller the strain energy could contribute from the buried layers.

43

z (001)

z (113)

z (111) y (010)

y (-110)

y (-110)

x (100)

x (33-2)

x (11-2)

GaAs (001)

GaAs (113)

GaAs (111)

Figure 4.1 Schematic illustration of crystal coordinates

5

(a)

5 4.4 4 3.2

1

2.8 2.4 2

-1

1.6

1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

3

3.6

x(100) Grid

x(100) Grid

3

(b)

1 -1

1.2

-3

0.8

-3

0.4

-5

-5

0

-5

-3

-1

1

3

5

-5

-3

y(010) Grid

(c)

5 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

x(11-2) Grid

3 1 -1 -3 -5 -5

-3

-1

1

y(-110) Grid

1

3

5

3

5

(d)

4 3.6

3 x(33-2) Grid

5

-1

y(010) Grid

3.2 2.8

1

2.4 2

-1

1.6 1.2

-3

0.8 0.4

-5

0

-5

-3

-1

1

3

5

y(-110) Grid

Figure 4.2 Contours of normalized elastic strain energy distribution on the substrate surface of isotropic Iso(001) (a), and anisotropic GaAs(001) (b), GaAs(111) (c), and GaAs (113) (d). 44

Figure 4.3 Island ordering on the surface of GaAs without elastic strain energy (a), with strain energy of Iso (001) (b), with anisotropic strain energy of GaAs (001) (c), GaAs (111) (d), and of GaAs (113) (e). Different island chain orderings (red dashed lines) corresponding to different growth orientations can be observed. Other simulation parameters are T=750K, F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid.

45

Figure 4.4 Island density contours. The whole 200×200 grids are divided into 18×18 big grids, so there is 11×11 original grids in one big grid. The densities are plotted by 3 dimensional contours.

46

(a) GaAs (Iso)

(b) GaAs (001)

(c) GaAs (111)

(d) GaAs (113)

T=550K

T=650K

T=750K

Figure 4.5 Surface island distributions vs. temperatures (550K, 650K, 750K) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are F=1.0Ml/s, c=20%, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations.

47

Figure 4.6 Number of islands and average island size vs. temperature for islands growth on substrates with different orientations for short temperature range (550K to 750K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s.

Figure 4.7 Number of islands and average island size vs. temperature for islands growth on substrate GaAs (001) for long temperature range (550K to 950K). Other simulation parameters are fixed: F=1.0Ml/s, c=20%, and ti=200s. 48

(a) GaAs (Iso)

(b) GaAs (001)

(c) GaAs (111)

(d) GaAs (113)

c=10%

c=20%

c=30%

c=50%

Figure 4.8 Surface island distributions vs. coverage c (10%, 20% 30% and 50%) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and ti=200s, on a 200×200 grid. Red dashed lines indicate the island chain orientations.

49

Figure 4.9 Number of islands and average island size vs. coverage for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and ti=200s.

50

(a) GaAs (Iso)

(b) GaAs (001)

(c) GaAs (111)

(d) GaAs (113)

ti=0s

ti=10s

ti=200s

Figure 4.10 Surface island distributions vs. interruption times ti (0s, 10s, and 200s) on different substrates: Iso (001) in (a), GaAs (001) in (b), GaAs (111) in (c), and GaAs (113) in (d). Other simulation parameters are T=750K, F=1.0Ml/s, and c=20%, on a 200×200 grid. Red dashed lines indicate the island chain orientations.

51

Figure 4.11 Number of islands and average island size vs. interruption time for islands growth on substrates with different orientations. Other simulation parameters are fixed: T=750K, F=1.0Ml/s, and c=20%.

Third layer Second layer First layer Substrate

Figure 4.12 Schematic illustration of superlattice.

52

(a) First layer

(b) Second layer

(c) Third layer

Figure 4.13 Atom island ordering on the first three layers. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Thickness of each layer is 5 grids. Strain energy field is calculated at every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for growth interruption.

(a) First layer

(b) H=5 grids

(d) H=15 grids

(e) H=20 grids

(c) H=10

(f) H=30

Figure 4.14 Atom island ordering for different layer thickness. Figures (b) to (h) are the simulated results on the second layer with thickness equal to 5 grids to 30 grids. Simulation parameters are T=700K, F=0.01Ml/s, c=20% and t=70s for each layer on a 200×200 grid. Strain energy field is updated every 2500 steps for each layer. It takes each layer 20s to deposit atoms on the surface and 50s for interruption.

53

CHAPTER V

KINETIC MONTE CARLO SIMULATION OF THREE DIMENSIONAL QUANTUM DOT SELF-ASSEMBLY

5.1 Overview of Three Dimensional Quantum Dots Growth

Self-assembled quantum dots (QDs) are intensively studied in recent years due to their optical and electronic properties with potential applications in optoelectronics and semiconductor devices [Bre95, Bim98]. InAs/GaAs heterostructure is a typical example, which is characterized by a large lattice mismatch and undergoes a growth transition from a two-dimensional (2D) cluster to three-dimensional (3D) islands [Sta04, Joy04]. Their optical properties in the 2D cluster and 3D islands were also studied [Pol06]. Experimental data indicate that InAs/GaAs island array is governed by thermodynamics [Shc99, Pch97]. The physics behind the transition from 2D cluster to 3D island has been studied recently. For example, using the atomic force microscopy, Arciprete et al.[Arc06] studied how kinetics drives the 2D to 3D transition in InAs/GaAs(001) heterostructure. The transition from thermodynamically to kinetically controlled QDs self-assembly was also studied by Musikhin et al. [Mus05], both 54

experimentally and theoretically. It is argued that strains on the surface could change diffusivity and surface mobility [Zoe00, Sag00, Zan99] and therefore could be the reason for the enhanced evaporation [Sun00]. The epitaxial system can lower its free energy by transferring atoms from the island edge to the upper layer, because the transition leads to a decrease in the contact area between the substrate and a new layer [Arc06]. As such atomic transitions to the upper layers lead to the strain field relaxation. The tradeoff between the cost of additional surface energy and the gain of energy due to elastic relaxation is the very driving force for the transition from 2D cluster to 3D islands [Ker97, Cha00]. These arguments provide reason enough to consider situations when the atom hop probability to an upper layer can be higher than that to a lower layer. The kinetics of such a process could be described by adjusting the probability for atoms to hop up or hop down, which could subsequently lead to 2D cluster or 3D island self-assembly [Arc06, Pol06]. While some simple computational models for 2D to 3D self-assembly were discussed before (i.e., [Bru01]), no complete 3D QDs self-assembly model has been developed so far in which the growth from 2D cluster to 3D islands and the corresponding island size equalization can be clearly illustrated during growth process. We, therefore, propose a fast multiscale 3D QDs growth model. It is developed from our original (x,y)-plane growth model [Pan04, Pan06] by a fast algorithm for the long-range strain energy calculation and by introducing the up/down ratio for the lateral self-organization. With the proposed program, we have successfully simulated the transition of the QDs 2D cluster to 3D island process. Furthermore, depositions with different flux rates are studied to show how the flux rate affects the island equalization during deposition 55

process. Our model has also shown clearly the importance of the interruption time during the growth of QD islands. We believe that this model will provide an attractive means for producing ordered nanostructures.The 2D to 3D self-assembly was discussed by some computer modeling [Bru01]. But, so far, no complete 3D QDs self-assembly model has been developed to deal with 2D to 3D transition and island size equalization during growth process. Here, we developed a 3D QDs growth model from our original (x,y)plane growth model [Pan04, Pan06] to study QDs self-assembly in the real 3D growth situation. By adjusting the up/down ratio ρ and considering the lateral self-organization, we simulate the QDs 2D to 3D transition process and island size equalization. Furthermore, examples with different deposition flux rate are studied to show how it affects the island equalization during deposition process. We believe that this model will provide an attractive means for producing ordered nanostructures.

5.2 Three Dimensional Model Description

The 3D layer-by-layer growth model is developed from our 2D (x,y)-plane growth model [Pan04, Pan06]. As in our previous model, the important contribution from the long-range elastic strain energy is included using a fast algorithm based on the Green’s function method [Pan02]. Furthermore, to account for the out-of-plane (up and down) movement of the atoms, which is also called desorption and adsorption [Lam01, Lun05], an up/down ratio is introduced. First, the 2D hopping probability of an atom from one lattice site to a nearest or next nearest neighbor site in the (x,y)-plane is still governed by the Arrhenius law 56

enhanced by the long-range strain energy field [Nur01, Lar01, Pan04, Pan06]. ⎛ E + E n − Estr ( x, y ) ⎞ ⎟⎟ p = ν 0 exp⎜⎜ − s k BT ⎠ ⎝

(5.1)

where ν0 is the attempt frequency (=1013s-1), T the temperature, and kB the Boltzmann’s constant. Also in equation (1), Es and En are the binding energies to the surface and to the neighboring atoms, respectively. Finally, Estr(x,y), as a function of the plane coordinates (x,y), is the energy correction from the long-range strain field due to the lattice mismatch between the substrate and the deposited material. As this long-rang strain energy needs to be calculated repeatedly during the simulation, we have developed a fast algorithm by pre-calculating the energy along a unit circle interpreting the energy at any location afterwards [Pan03]. We calculate the binding energy to the neighboring atoms in the (x,y)-plane, En, using the following algorithm: We assume that the strength of a single nearest neighbor bond is Eb (=0.3eV), and it is reduced by a factor of α (= 1

2 ) for the next nearest

neighbors. To evaluate the diffusion barrier, the binding energy at the site P0, where the diffusing atom is located, is calculated to be EP0 = nEb + αmEb

(5.2)

with n≤4 and m≤4 being, respectively, the number of nearest and next nearest atoms. Similarly, for the site P1 where the atom is going to hop to, we have EP1 = g (n' Eb + αm' Eb )

(5.3)

where n′≤4 and m′≤4 are, respectively, the number of nearest and next nearest atoms at the new site P1, and g (=0.2) describes the coupling between adjacent lattice sites [Mei03]. 57

Therefore, the overall binding energy En caused by neighbor interactions for a diffusion process from site P0 to site P1 in the (x,y)-plane is given by the difference of the binding energies at the corresponding lattice sites E n = (n − g ⋅ n')Eb + (m − g ⋅ m')αEb

(5.4)

Second, the binding energy to the surface of the (x,y)-plane, i.e. Es, was assumed to be constant (actually, Es=1.3eV) in our previous 2D (x,y)-plane self-assembly model [Pan04, Pan06]. However, in 3D case (figure 5.1), the binding energy to the surface depends on the surface geometry of both the initial and final position of the atoms. Based on the recent molecular dynamics calculation [Mon01] and KMC simulation experience [Mei03, Pan04, Pan06] in 2D, we therefore propose the following simple equation for the 3D surface binding energy calculation. Es = (1 − g )Es0 + ( p − g ⋅ p')αEs0 + (q − g ⋅ q')α ⋅ αEs0

(5.5)

where Es0 is the contribution from the atom exactly under the selected adatom (figure 5.1). p, q are number of nearest and next nearest atoms in original positions (n'≤4, m'≤4). p', q' are number of nearest and next nearest atoms in new positions (n'≤4, m'≤4). The other parameters, such as g and α, are used to decrease the contribution from the atoms in the first and second squares (figure 5.1a, 5.1b). Assuming that the maximum surface binding energy be 1.3eV as before for the 2D growth simulation [Pan04, Pan06], which means that all the positions under this adatom are occupied by atoms, we can find Es0=0.28eV from equation (5.5) by back-calculation. Finally, as discussed in the introduction, the growth system always tries to decrease its free energy by moving atoms at the edge to upper layer to form 3D islands. The possibility for an atom to hop to upper or lower layers depends on the material 58

properties and growth condition [Arc06, Pol06]. To account for this up and down jump activities of the atoms, the parameter up/down ratio ρ is introduced, which is defined as the ratio of the probability for edge atoms hopping up to the probability hopping down. In other words, if we let Pup be the jumping up possibility and Pdown the jumping down possibility, then the up/down ratio ρ=Pup/Pdown. It is further remarked that the up/down ratio actually reflects the balance between surface energy increase and strain energy decrease. This physical activity is illustrated in figure 2 where layer 1 is the substrate and the atom “A” is on the edge of layer 2 (within the QD island). It is also understandable that the up or down jump possibility of atom “A” is controlled only by the up/down ratio instead of by the individual jumping up and down possibilities. Moreover, strictly speaking, since the strain energy changes along the island height direction, the up/down ratio is, in general, not constant in different island layers. However, in order to extract the important contribution of the up/down ratio, we assume that the up/down ratio ρ is constant. In other words, it is only affected by the geometry around it but not by the layer position. Furthermore, in order to form 3D islands, the number of atoms jumping up should be larger than those jumping down, which means that the up/down ratio ρ should be larger than 1. Otherwise, all the atoms will move to lower layer to form the layer-bylayer Frank-van der Merwe thin-film structure.

5.3 Simulation Results

Using the 3D QDs self-assembly approach described above, we can now simulate the epitaxial growth process. The growth model is on the 100×100 grids with 59

periodic boundary condition, and the material is InAs/GaAs (001).

5.3.1 Height Dependence of Up/Down Ratio (ρ) We first study the effect of the up/down ratio on the island height. Shown in figure 5.3 are the island distributions for different up/down ratio ρ. It is observed clearly from figure 5.3 that 1) the average island height increases with increasing up/down ratio ρ (for small up/down ratio, say, ρ30); 3) the inflection point within the sharp-change range, approximately at ρ=13, 60

physically corresponds to the critical transition point from 2D growth to 3D growth. This value is equivalent to an energy ratio between surface and bulk as has been demonstrated in [Bru01, Kra06] using a different approach.

5.3.2 Size Dependence of Up/Down Ratio (ρ) As we know, the island size equalization relies on the movement of atoms on the surface during epitaxial growth. The equalization starts to happen from the very beginning of deposition process and finally it will reach the equilibrium status if long enough interruption time after deposition [Pan04] is given. Here, three examples with different up/down ratios of 10, 20 and 30 have been run to show how the island size equalization happens during the deposition and interruption process (figure 5.5). For each example, the island distribution of exactly after deposition (ti=0), 100 seconds after deposition (ti=100) and 200 seconds after deposition (ti=200s) are shown in figure 5.5. It is easy to find that: 1) for ρ=10, 20 and 30, increasing the interruption time, small isolated islands assembles to big ones, and, with long interruption time, the island distribution is getting increasingly ordered; 2) by adjusting the up/down ratio, people can get island distribution with different shape, which could be a way for QDs growth experimentalists to get the optimized QDs distribution or shape by adjusting the up/down ratio. Taking the first example with up/down ratio of 10, in figure 5.5, we analyze the relationship between island diameter and number of island within this diameter (figure 5.6 and figure 5.7) using the simulated results of 0.5 ML, 1 ML, 1.6 ML deposition and 200s interruption after deposition. We find some remarks. 1) Island size equalization 61

happens starting from the beginning of deposition process for all three up/down ratio

ρ=10, 20 and 30, and, with long enough interruption time, finally we will get ordered island distribution with similar size; 2) In this example, low flux rate (F=0.01 ML/s) is used, so island equalization is obvious even just at the end of the deposition process; 3) We believe that if we decrease the flux rate, we will get the better island distribution during deposition process. To further confirm our conclusion in 3), more examples with same growth parameters except for different flux rate are run in figure 5.8. Finally, figure 5.8 shows the effect of the flux rate (F=1Ml/s, 0.1ML/s and 0.01ML/s) on the island distribution. For a lower flux rate (e.g. at F=0.01Ml/s), the deposited atoms will have more time to move to the equilibrium position during the deposition process and to assemble together. Furthermore, a low flux rate (say, F=0.01 ML/s) usually corresponds to large and ordered islands. Therefore, besides the up/down ratio, flux rate is also important in order to obtain ordered and uniform-size QD island distribution.

5.4 Conclusion of 3D Simulation Model

In this chapter, we propose a 3D QDs epitaxial growth model by updating our former (x,y)-plane growth model. Specifically, the balance between surface energy increase and strain energy decrease is demonstrated by introducing the parameter up/down ratio. Furthermore, island height dependence and island shape dependence of up/down ratio have been studied too. Combining up/down ratio and one of the growth parameters flux rate, island equalization during deposition process and after deposition is studied. The main conclusions from our 3D KMC model are as follows: 62

1.

The phenomena of self-assembly and shape transition of the island array can be simulated with changing parameter up/down ratio.

2.

For other fixed growth parameters, island height and size greatly depends on the up/down ratio. For the examples used in this paper, the average island height increases with increasing up/down ratio and island size and shape are very sensitive to the up/down ratio when it is less then 20. Furthermore, the critical transition point of up/down ratio in this example is around 13.

3.

It is the up/down ratio, the ratio characterizes the atom jump possibility to upper layer over that to lower layer, that determines the transition from 2D cluster to 3D islands, instead of the absolute jump possibility value.

4.

Island size equalization starts from the very beginning of deposition process. With enough interruption time, ordered island distribution with uniform size can be obtained.

5.

The lower the flux rate, the better the island distribution during deposition process.

63

(a)

(c)

(b)

Es y

y x x Figure 5.1 Schematic illustration of 3D QD self-assembly model. An atom on top of the substrate surface (x,y)-plane in (a), the relative locations of the atom grid on the (x,y)plane in (b), the corresponding binding energy Es related to the in-plane locations in (c).

4 3 2 1

A

Figure 5.2 Illustration of an edge atom A’s jump up or down during 3D QD self-assembly.

64

(a) ρ=3

(b) ρ=5

(c) ρ=7

(d) ρ=10

(e) ρ=20

(f) ρ=30

Figure 5.3 3D island distributions for different up/down ratios ρ with total coverage c=1.6Ml. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and interruption time ti=200s (The total simulation time = the deposition time (td=160s) plus the interruption time (ti=200s)).

Average Island Height (Grids)

20 16 12 8 4 0 0

10 20 30 40 Up/Down Ratio (ρ)

50

Figure 5.4 Average island height vs. up/down ratio ρ (for ρ from 3 to 50). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and interruption time ti=200s. 65

ti =0

ti =100s

ti =200s

ρ=10

ρ=20

ρ=30

Figure 5.5 Island distributions for different interruption times (ti=0s, corresponding to deposition time td=160s, ti=100s, and ti=200s) and for different up/down ratios (ρ=10, 20 and 30). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, and total coverage c=1.6Ml.

66

(a) td=50s

(b) td=100s

(c) td=160s; ti=0

(d) ti=200s

Figure 5.6 Island distributions during and after deposition. Deposition time td=50s with 0.5Ml coverage in (a), td=100s with 1Ml coverage in (b), td=160s with 1.6Ml coverage (i.e., at the beginning of interruption time) in (c) and interruption time ti=200s in (d). Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml, and up/down ratio ρ=10.

67

40

40

(a) td=50s

30 Num of Islands

Num of Islands

30

20

10

0

0

2

0

4 6 8 10 12 14 16 18 Island Diameter (Grids)

0

2

4 6 8 10 12 14 16 18 Island Diameter (Grids)

(d) ti=200s

16

(c) td=160s; ti=0

20

12 Num of Islands

Num of Islands

20

10

25

15 10

8

4

5 0

(b) td=100s

0

2

0

4 6 8 10 12 14 16 18 Island Diameter (Grids)

0

2

4 6 8 10 12 14 16 18 Island Diameter (Grids)

Figure 5.7 Average island diameter (of the bottom equivalent circle of the island) vs. number of islands during and after deposition. Fixed growth parameters are temperature T=700K, flux rate F=0.01Ml/s, total coverage c=1.6Ml and up/down ratio ρ=10.

68

0.5 ML

1 ML

1.6 ML

1 ML/s

0.1 ML/s

0.01 ML/s

Figure 5.8 Island distributions for different deposition processes of 0.5Ml, 1Ml and1.6 Ml with flux rates F=1Ml/s, 0.1Ml/s and 0.01Ml/s. Fixed growth parameters are temperature T=700K, total coverage c=1.6Ml and up/down ratio ρ=10.

69

CHAPTER VI

MOLECULAR DYNAMICS CALCULATION OF ELASTIC MODULI IN STRAINED SEMICONDUCOTR

6.1 Background of Elastic Moduli in Strained Semiconductor

The quantum dots growth and the effect of anisotropic strain energy on island distribution have been studied in the former chapters. As we discussed, the anisotropic strain energy is caused by the misfit strain and basically is because of the unit crystal side length difference between the substrate and deposited material. Take InAs/GaAs system as an example [Ell02], usually we take the InAs QD and compress it uniformly by about 6.685% in each direction so that it has the same lattice constant as that of GaAs. It is then placed in a hole of exactly the right size, so that the InAs and GaAs crystals are coherent, and then the uniform compression is removed, which leads to straining of the surrounding matrix. Though it is not easy to accurately determine the strain distribution, there are some methods available [Dav99, Hol99, Pan01]. The strain effects of semiconductors have been widely investigated and lots of reports about influence on electronic and optical properties are available [Gle99, 70

Jai00, Jog98, Par99], but much less is known about their influence on the elastic constants. In the recent years, people begin to study the elastic constants of strained semiconductors. For example, Mehl [Meh93] studied moduli of Al-Li compounds under pressure, and Kanoun et al. studied the elastic properties of BN, AlN, GaN and InN under pressure effect by using ab initio calculations [Kan04a, Kan04b, Mer03]. And, recently, Ellaway and Faux [Ell02, Ell03] reported the effective elastic constants of InAs under uniform strain by using DLPOLY [Smi96] molecular dynamics (MD) package. We develop a MD program by using Stillinger-Weber potential. To check our program, using the parameters for InAs [Ell02, Ich96], elastic stiffnesses of InAs determined by our program match the results by Ellaway and Faux [Ell02] very well. Crystalline silicon is widely used in micro-electro-mechanical-systems (MEMS) and it is one of the most intensely investigated materials. Interest in silicon is mostly motivated by its great technological importance. Strained silicon has its own specific meaningful use. In 2002, Intel announced their employing a unique strained silicon transistor technology on its 90 nm logic process. This was the first process in the industry to implement strained silicon in production. This process improves transistor performance by 10-25 percent while increasing manufacturing costs by only two percent. This is because strained silicon increases transistor drive current which improves switching speed by making current flow more smoothly [Gha03]. So, to understand the mechanism of strained silicon well, to study the elastic moduli of strained silicon is very important. In this chapter, the effective moduli of silicon under hydrostatic strain are calculated as the functions of volumetric strain. While variation of the elastic moduli was studied previously when silicon was under pressure [Mcs64, Kar97], no literature is 71

available on the stress and strain behavior when silicon is under strain. Furthermore, due to the technical difficulty, it is hard to observe experimentally the elastic properties as functions of the applied strain. Here, we present a molecular dynamic prediction for the elastic stiffnesses in strained silicon as functions of the volumetric strain level. Our approach combines the basic continuum mechanics with the classical molecular dynamic approach, supplemented with the interatomic S-W potential [Sti85]. This potential has gained popularity due to its ability to describe fairly well diamond structures and has been used in numerous studies (e.g., [Klu86, Hum04, Men04]). In order to use our approach for the estimation of the three independent elastic stiffnesses C11, C12 and C44 in strained silicon, we apply three independent distortions as described by Mehl [Meh90] and utilized later by Ellaway and Faux [Ell02] for InAs. The elastic stiffnesses and effective Young’s modulus and Poisson’s ratio are all calculated and presented in terms of figures and formulas. In general, our simulation indicates that the bulk modulus, C11 and C12, increase with increasing volumetric strain whilst C44 is almost independent of volumetric strain. The difference between the strained moduli and the ones at zero strain can be very large, and therefore use of the standard free-strained moduli in the strained silicon case could result in wrong mechanical response of the silicon.

6.2 Elastic Stiffness Calculation

The III-V semiconductors are cubic crystals with three independent elastic stiffnesses C11, C12 and C44. An efficient way of calculating elastic constants is 72

important because these constants are directly employed in practical uses of materials. There are several ways to calculate elastic constants. The direct and traditional method is to apply a tension on the sample and calculate the corresponding strain and elastic constants from the tension-strain relationship [Ell02, Meh93, Ray88]. For a cubic crystal, the elastic moduli can be divided into two classes: the bulk modulus B=(C11+2C12)/3 and two “shear” moduli, C11-C12 and C44. S In the actual calculation, a small incremental distortion, which is usually less than 1%, is applied to the crystal with the latter being already subjected to a uniform volumetric strain (i.e., [Ell02]). Let’s assume the original cubic lengths and the distortions in x, y and z directions are εx, εy, εz and ∆x, ∆y, ∆z. After deformation, the new volume is: V = ∆x(1 + ε x )∆y (1 + ε y )∆z (1 + ε z )

= ∆x∆y∆z (1 + ε x + ε y + ε z + ε xε y + ε y ε z + ε z ε x + ε xε y ε z )

(6.1)

For small strains, the high order terms can be neglected and the volume simplifies as: V ≈ ∆x∆y∆z (1 + ε x + ε y + ε z ) = V0 (1 + ε x + ε y + ε z )

(6.2)

Thus, the volumetric strain becomes e=

V − V0 =εx +εy +εz V0

(6.3)

where V0 is the original volume of the unit cell, V is the volume after deformation and e is volumetric strain. Constitutive relation for cubic material is

73

⎡σ xx ⎤ ⎡C11 ⎢σ ⎥ ⎢ ⎢ yy ⎥ ⎢C12 ⎢σ zz ⎥ ⎢C12 ⎢ ⎥=⎢ ⎢ τ yz ⎥ ⎢ 0 ⎢ τ zx ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎣⎢ τ xy ⎦⎥ ⎣⎢ 0

C12

C12

0

0

C11

C12

0

0

C12

C11

0

0

0

0

C 44

0

0

0

0

C 44

0

0

0

0

0 ⎤ ⎡ε xx ⎤ ⎢ ⎥ 0 ⎥⎥ ⎢ε yy ⎥ 0 ⎥ ⎢ε zz ⎥ ⎥⎢ ⎥ 0 ⎥ ⎢γ yz ⎥ 0 ⎥ ⎢γ zx ⎥ ⎥⎢ ⎥ C 44 ⎦⎥ ⎣⎢γ xy ⎦⎥

(6.4)

And energy density expression for this system is 1 U = (σ xx ε xx + σ yy ε yy + σ zz ε zz + τ yz γ yz + τ zx γ zx + τ xy γ xy ) 2

(6.5)

In the following, we briefly list the relevant formulas needed to obtain the elastic constants of cubic crystal. The stress tensor and the elastic constants can be defined in different ways and it is useful to define our notation. We consider both strain with and without volume conservation, because volume-conserving strains alone do not provide enough constraints to determine all three elastic constants of cubic crystal. Only small lattice distortions are considered in order to remain within the elastic domain of the crystal.

6.2.1 Bulk Modulus Bulk modulus is related to the E(V ) curvature [Bir78],

B(V ) = −VP' (V ) = VE " (V )

(6.6)

where V is the volume of unit cell as mentioned before, E(V) is the energy/unit cell at volume V, E″(V) is the second derivative of E(V) with respect to the volume, and P(V) = – E′(V) is the pressure required to keep the cell at volume V. Since the calculations only provide a set of energies E(Vi) for a limited number of volumes Vi, the second derivative

E″(V) must be approximated [Meh93]. Birch [Bir78] proposed the following form 74

to calculate E(V) by fitting the cure. 2

3

N ⎡⎛ V0 ⎞ 2 / 3 ⎤ ⎡⎛ V0 ⎞ 2 / 3 ⎤ ⎡⎛ V0 ⎞ 2 / 3 ⎤ 9 9 ' E (V ) = E0 + B0V0 ⎢⎜ ⎟ − 1⎥ + B0V0 B0 − 4 ⎢⎜ ⎟ − 1⎥ + ∑ γ n ⎢⎜ ⎟ − 1⎥ 8 ⎥⎦ 16 ⎢⎣⎝ V ⎠ ⎥⎦ n=4 ⎣⎢⎝ V ⎠ ⎥⎦ ⎣⎢⎝ V ⎠

(

)

n

(6.7) where E0, V0, B0, and B0′ are, respectively, the equilibrium energy, volume, bulk modulus, and pressure derivative of the bulk modulus, whilst N is the curve fitting order. We remark that equation (6.7) is a special case of the following general expression N

E (V ) = ∑ anV −2 n / 3

(6.8)

n =0

where an are the fitting parameters and we find that a third polynomial is sufficient (N=3). In our simulation, the distortion used for calculating bulk modulus consists of a uniform compression (or expansion) of magnitude δ in each of x, y and z directions as in a conventional pressure experiment. The distortion can be represented by the matrix ⎛δ 0 0 ⎞ ⎜ ⎟ ⎜0 δ 0⎟ ⎜0 0 δ ⎟ ⎝ ⎠

(6.9)

6.2.2 Modulus C11-C12 In case of cubic of lattice, it is possible to choose the strain that keeps the volume of unit cell preserved. “Shear” modulus C11-C12 is evaluated by applying an orthorhombic volume-conserving (up to the first order of δ) strain, which can be represented in the strain tensor form as (proved in Appdendix C)

75

⎛ ⎜δ 0 ⎜ ⎜0 −δ ⎜⎜ 0 0 ⎝

⎞ ⎟ ⎟ ⎟ 2 δ ⎟ ⎟ 1−δ 2 ⎠ 0 0

(6.10)

Making use of the stress and strain relation, we found that the non-zero stresses are ⎡σ xx ⎤ ⎡C11 ⎢σ ⎥ ⎢ ⎢ yy ⎥ ⎢C12 ⎢σ zz ⎥ ⎢C12 ⎢ ⎥=⎢ ⎢ τ yz ⎥ ⎢ 0 ⎢ τ zx ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎢⎣ τ xy ⎥⎦ ⎣⎢ 0

C12

C12

0

0

C11

C12

0

0

C12 0

C11 0

0

0 0

0

0

C44 0

0

0

0

C44 0

0 ⎤⎡ δ ⎥ ⎢ 0 ⎥ −δ ⎢ 2 0 ⎥ ⎢δ / 1 − δ 2 ⎥⎢ 0 ⎥⎢ 0 0 ⎥⎢ 0 ⎥⎢ C44 ⎦⎥ ⎣⎢ 0

(

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦⎥

)

(6.11)

Then we have

δ2 1− δ 2 δ2 σ yy = (C12 −C 11)δ + C12 1− δ 2 δ2 σ zz = C11 1− δ 2 τ yz = τ zx = τ xy = 0 σ xx = (C11 −C 12 )δ + C12

(6.12)

So, the system energy as a function of δ is

E (δ ) = E (0 ) + (C11 − C12 )Vδ 2 + 0(δ 4 ) + ...

(6.13)

where terms of order δ4 and higher order can be safely neglected. So, the energy E(δ) may be found as a function of δ2 and a linear curve fitting yields C11-C12.

6.2.3 Elastic Stiffness C44 Elastic stiffness C44 is evaluated by a monoclinic volume-conserving strain (proved in Appendix D). The strain can be represented by the matrix 76

⎛ ⎜ 0 ⎜ ⎜1 ⎜2δ ⎜ ⎜⎜ 0 ⎝

1 δ 2

0

0

0

0

δ2 (4 − δ 2 )

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎟ ⎠

(6.14)

and, the stress and strain relation becomes: ⎡σ xx ⎤ ⎡C11 ⎢σ ⎥ ⎢ ⎢ yy ⎥ ⎢C12 ⎢σ zz ⎥ ⎢C12 ⎢ ⎥=⎢ ⎢τ yz ⎥ ⎢ 0 ⎢ τ zx ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎢⎣τ xy ⎥⎦ ⎣⎢ 0

C12

C12

0

0

C11

C12

0

0

C12

C11

0

0

0

0

C44

0

0

0

0

C44

0

0

0

0

0 ⎤⎡ 0 ⎥⎥ ⎢ ⎢ 0 ⎥ ⎢δ 2 / ⎥⎢ 0 ⎥⎢ 0 ⎥⎢ ⎥⎢ C44 ⎦⎥ ⎣⎢

⎤ ⎥ 0 ⎥ 4−δ 2 ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ δ ⎦⎥ 0

(

)

(6.15)

Then we have

σ xx = σ yy = σ zz = 0 τ yz = τ zx = 0

(6.16)

τ xy = δ So, application of this strain changes the total energy from its unstrained value to

( )

1 E (δ ) = E (0) + C 44Vδ 2 + 0 δ 4 + ... 2

(6.17)

Similarly, the energy E(δ) can be found as a function of δ2 and a linear curve fitting yields C44.

6.3 Molecular Dynamics of Elastic Moduli Calculation

Molecular dynamics is a technique to compute the equilibrium and transport properties of a classical many-body system.

77

6.3.1 Classical Interatomic Potential In order to use molecular dynamics, we have to define the rules that are governing interaction of atoms in the system. In classical and semi-calssical simulations, these rules are often expressed in terms of potential functions. The potential function U(ri) describes how the potential energy of a system of N atoms depends on the coordinates of the atoms. In the past, various empirical potentials have been proposed to describe the different phases of silicon. A number of these empirical potentials [Ter88, Sti85, Bas89, Bis85] have been proposed to describe the different phases of Si and the results have been tested and compared [Cow88, Coo93, Bal92]. Among these, the Stillinger and Weber (S-W) potential [Sti85] has gained popularity due to its ability to describe fairly well diamond structures and has been used in numerous studies [Klu86, Hum04, Men04] since it was proposed. The main advantage of this potential is its simplicity and fairly realistic description of crystal. The potential consists of a two-body term and an explicit angular interaction: V = ∑ν 2 (rij ) +

∑ν (r , r , r ) (r ) = εf (r σ ) (r , r , r ) = εf (r σ , r σ , r i< j

ν2 ν3

ij i

2

j

i< j