1 STRUCTURE AND DYNAMICS OF IN TERFA CE S IN ORGANIC AND INORGANIC MATERIALS USING ATOMIC LEVEL SIMULATION By DONGHWA LEE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010
2 2010 D onghwa L ee
3 To my f amily with love
4 ACKNOWLEDGMENTS I want to show my gratitude to my advisor Prof. Simon R. Phillpot for his kind instruction and all effort s to gui de me during Ph. D study. His continuous support s, useful suggestions and helpful discussions have made me pass through all difficulties during my research. I also would like to thank Prof Susan B. Sinnott for her advices for resolving lots of issues for the study I also wish to express my apprec iation to committee members Prof. Aravind Asthagiri, Prof. David Nor ton and Prof. Jacob L Jones. I cannot resist mention ing a ll computational science focus group member s for their kindness, friendship and supports. From the beginning of my study, they are always my mentors, friends and resources for all kinds. Without their help, I cannot get to this stage. I personally express my appreciation to my family My parents have been always mental supporters in my s ide I deeply thanks to m y wife, Minjung Kim for all she has made. I never forget how much efforts she spent for me during my graduate study periods.
5 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................. 4 LIST OF TABLES ............................................................................................................ 8 LIST OF FIGURES .......................................................................................................... 9 ABSTRACT ................................................................................................................... 13 CHAPTER 1 GENERAL INTRODUCTION .................................................................................. 16 1.1 Imperfections I n S olids ...................................................................................... 16 1.1.1 P oint D efects (0 D) .................................................................................. 17 1.1.2 L ine D efects (1 D) ................................................................................... 18 1.1.3 Interfacial or P lanar D efects (2 D) ........................................................... 19 1.1.4 Bulk or V olume D efects (3 D) .................................................................. 20 1.2 Interface S cience .............................................................................................. 20 2 SIMULATION METHODS ....................................................................................... 22 2.1 Simulation Methodologies ................................................................................. 22 2 2 First Principles Calculations .............................................................................. 23 2.3 Density Functional Theory (DFT) ...................................................................... 26 2 3.1 Hohenberg Kohn T heorem s .................................................................. 26 2.3 2 KohnSham Self Consistent Field Methodology .................................... 27 2 3 3 ExchangeCorrelation Functional .......................................................... 30 2.3. 4 Pseudopotentials ................................................................................... 31 2 3 5 Planew ave Implementation .................................................................. 32 2.3. 6 Linear Scaling ....................................................................................... 33 2.3. 7 Application to LiNbO3 S tudy .................................................................. 35 2. 3. 8 Application to Molecular B eam D eposition ............................................ 36 2.3. 9 Charge analysis M ethods ...................................................................... 37 2 .3. 10 Limitations on DFT ................................................................................ 39 2 4 Molecular Dynamics (MD) Simulation ............................................................... 40 2.4.1 Introduction .............................................................................................. 40 2.4.2 Interatomi c Interactions ........................................................................... 41 2. 4 .3 Shell Model for Describing P olarizability of E lectrons .............................. 44 2. 4 .4 Threebody I nteractions ........................................................................... 45 2. 4 .5 Potential for LiNbO3 study ................................ ....................................... 45 2. 4 .6 Integration Algorithms .............................................................................. 46 2. 4 .7 Geo metry O ptimization ............................................................................ 49 2. 4 .8 Thermodynamic Conditions ( Ensemble s) ................................................ 52
6 3 STRUCTURAL ANALYSIS ..................................................................................... 54 3.1 Trigonal Ferroelectrics ...................................................................................... 54 3.1.1 Introduction on LiNbO3 ............................................................................ 54 3.1.2 Crystallography of LiNbO3 ....................................................................... 55 3.1.3 Temperature E ffects on LiNbO3............................................................... 63 3.1.4 Domain W alls in LiNbO3 .......................................................................... 65 3.2 Conductive P olymers ........................................................................................ 67 3.2.1 Introduction .............................................................................................. 67 3.2.2 Structure of P olythiophene ...................................................................... 69 3.2.3 Conduction M echanism in P olythiophene ................................................ 71 3.2.4 Growth of P olythiophene ......................................................................... 73 4 DOMAIN WALLS IN FER ROELECTRIC LITHIUM NIOBATE ................................. 76 4.1 Introduction ....................................................................................................... 76 4. 2 Domain W all P reparation .................................................................................. 77 4.3 Domain W all E nergy ......................................................................................... 78 4.4 Atomic D isplacement around D omain W alls ..................................................... 80 4.4.1 Y walls, P arallel to c glide P la nes ........................................................... 81 4.4.2 X walls, P erpendicular to c glide P lane ................................................... 85 4.5 Polarization at Domain Walls ............................................................................ 89 4.5.1 Y wall s ..................................................................................................... 90 4.5.2 X walls ..................................................................................................... 92 4.6 Domain W all W idth ........................................................................................... 94 4.7 Threshold F ield for D omain W all M otion ........................................................... 96 4.8 Conclusions and Discussion ............................................................................. 98 5 DOMAIN WALL DEFECT INTE RACTIONS IN FERROELECTRIC LITHIUM NIOBATE .............................................................................................................. 100 5 .1 Introduction ..................................................................................................... 100 5.2 Methodology and Validation ............................................................................ 100 5.2.1 Empirical S tudy with T wo D ifferent M ethods ......................................... 100 5.2.2 DFT with T hermodynamic F rame W orks ............................................... 102 5.2.3 Comparison of DFEs of N eutral D efect C omplexes .............................. 103 5.3 Configurations of T wo D ifferent I ntrinsic D efect P airs ..................................... 105 5.3.1 Li Frenkel P air (Lii + VLi ) ........................................................................ 105 5.3.2 Nb antisite P air(NbLi +4VLi ) ................................................................... 106 5.4 Domain W all intrinsic D efect I nteractions ....................................................... 109 5.4.1 Y wall with I ndividual P oint D efects ....................................................... 110 5.4.2 Y wall with Li Frenkel P air ..................................................................... 111 5.4.3 Y wall with AntisiteI (NbLi + 4VLi ) .......................................................... 112 5.4.4 X walls with P oint D efects ..................................................................... 113 5.4.5 X walls with Li Frenkel P air ................................................................... 114 5.4.6 Two X walls with AntisiteI P air ............................................................... 115 5.5 Extrinsic D efect, Mg ........................................................................................ 116
7 5.5.1 Reaction S chemes with Mg D opants ..................................................... 116 5.5.2 Domain W all Mg D opant I nteractions .................................................... 119 5.6 Conclusion ...................................................................................................... 124 6 DOMAIN WALLS IN LITHIUM TANTALATE ......................................................... 126 6.1 Introduction ..................................................................................................... 126 6.2 Potential D evelopment .................................................................................... 128 6.3 Validation of LiTaO3 P otential ......................................................................... 128 6.4 Energetics of D omain W alls in Lithium Tantalate ............................................ 132 6.5 Domain W all S tructure at T ransition S tate ...................................................... 136 6.6 Defects in Lithium Tantalate ........................................................................... 138 6 .7 C onclusion ...................................................................................................... 140 7 SURFACE DEPOSITION ON POLYTHIOPHENE MOLECULES ......................... 142 7.1 Introduction ..................................................................................................... 142 7.2 SPIAD ............................................................................................................. 143 7.3 C2H+ ion deposition on 3T polymer with PbS nanocrystalline ......................... 151 7.4 Conclusions .................................................................................................... 152 8 ADATIVE KINETIC MONTE CARLO (AKMC) ...................................................... 154 8.1 Introduction ..................................................................................................... 154 8.2 Mechanism ..................................................................................................... 155 8.2.1 Kinetic Monte Carlo (KMC) .................................................................... 155 8.2.2 Adaptive K inetic Monte Carlo (AKMC) ................................................... 157 8.3 D imer M ethod with T wo D imensional P otential ............................................... 158 8.4 Diffusion P rocess on (111) A luminum(Al) S urface .......................................... 159 8.5 Conclusion ...................................................................................................... 163 9 GENERAL CONCLUSIONS ................................................................................. 164 9.1 Domain W alls .................................................................................................. 164 9.2 Organic Gas I nterface ..................................................................................... 167 9.3 Metal gas I nterface ......................................................................................... 168 LIST OF REFERENCES ............................................................................................. 170 BIOGRA PHICAL SKETCH .......................................................................................... 187
8 LIST OF TABLES Table page 2 1 The potential parameters for LiNbO3. ................................................................. 46 3 1 Wyckoff positions for R 3 c ................................................................................... 57 3 2 Wyckoff positions for R3c ................................................................................... 60 3 3 Comparison of positional paramet ers of LiNbO3 between empirical calculations, DFT calculations and experimental data. ....................................... 62 3 4 The comparison on conductivity, stability and Processing attributes of representative conducting polymers ................................................................... 69 4 1 Comparison o f constants .................................................................................... 97 5 1 DFEs calculated by various methods. The interstitials are places at the center of a n empty oxygen cage. ...................................................................... 102 5 2 Lattice energy (eV) per unit formula of oxides for empirical study and DFT ..... 104 5 3 Lattice energy of oxides .................................................................................... 117 5 4 DFE of isolated defects: comparison of bulk DFEs calculated with a supercell method and the Mott Littleton method. ............................................................. 118 5 5 Energies for five reactions, calculated from the DFEs of isolated defects and the bulk lattice energies in Tables 52, 5 3 and 54. ......................................... 118 6 1 The comparison o f potential parameters between LiNbO3 and LiTaO3 ............ 128 6 2 Comparison of structural parameters and elastic properties of LiTaO3 between experiments and simulations .............................................................. 130 6 3 Intrinsic defect formation energies. ................................................................... 138 6 4 Frenkel, Schottky and pseudoSchottky defect formation energies .................. 139 6 5 Lattice energies of different oxide states .......................................................... 139 6 6 Energies of possible reactions that lead to Li2O deficiency .............................. 140
9 LIST OF FIGURES Figure pag e 1 1 Various classes of imperfections over the size range ......................................... 16 1 2 Geometrical production of dislocations ............................................................... 18 2 1 Flow chart of the Kohn Sham SCF procedure .................................................... 29 2 2 2D Schematic of the charge neutralization scheme used in the direct sum method ............................................................................................................... 43 2 3 Illustation of potential energy surface of a arbitrary function ............................... 50 2 4 Schematic of steepest descent method. ............................................................. 51 3 1 Crystallographic direction of rhombohedral axes and hexagonal axes in a unit cell of LiNbO3 ............................................................................................... 55 3 2 Schematic of the paraelectric and ferroel ectric phase ........................................ 58 3 3  projection of the rhombohedral structure ................................................... 61 3 4 Schematic of Y1, Y2 and Y3 axis along with three oxygen ions, O1, O2 and O3 ....................................................................................................................... 63 3 5 Temperature dependence of polarization of LiNbO3 from experiment and MD calculation .................................................................................................... 64 3 6 Schemat ic of surface charge associated with spontaneous polarization, and formation of 180 domains to minimize electrostatic energy ............................... 66 3 7 Schematic view parallel to the (0001) plane showing two str ucturally distinct domain wall structures ........................................................................................ 66 3 8 The comparison of calculated thiophene structure ............................................. 69 3 9 The structure of terthi ophene .............................................................................. 70 3 10 The variation of polythiophene polymer chain with ptype doping ...................... 72 3 11 The variation of electron band during oxi dation .................................................. 73 4 1 Structure of LiNbO3 ............................................................................................ 77 4 2 Energies of domain walls in equilibrium and in transition states. ........................ 80 4 3 Displacement pattern in z direction of oxygen ions and cations around Y wall .. 82
10 4 4 Atomic displacements of O1, O2 and O3 ions .................................................... 84 4 5 Crystallographies of the XIwall and XIIwall. ...................................................... 86 4 6 The Z displacements patterns around XIwall and XIIwall .................................. 87 4 7 Normal displacements of three different oxygen ions ......................................... 88 4 8 Comparison of z directional polarization (Pz) around Y walls, determined by DFT a nd empirical simulation methods .............................................................. 91 4 9 Schematic representation of three different polarization, Pz, Pt and Pn ............. 92 4 10 Comparison on Z directional polarization around X walls between DFT and empirical simulation methods ............................................................................. 93 4 1 1 A fit of Pz to hyper tangent function around Y wal l ............................................. 95 4 1 2 Fits of Pz of XIwall and XIIwall ........................................................................... 96 5 1 DFE comparison of stoichiometric and nonstoichiometric defect pairs in LiNbO3 using DFT and empirical method ......................................................... 105 5 2 Two different configurations of Li Frenkel pair .................................................. 106 5 3 Different structural arrangem ents of Nb antisite defect pair .............................. 107 5 4 The comparison of DFE and association energy .............................................. 108 5 5 The change in DFE of s ingle point defects near Y wall .................................... 111 5 6 Comparison on DEFs of Li Frenkel pair in two different configurations ............ 111 5 7 Formation energies of AntisiteI clusters as a distance from a Y domain w all. .. 112 5 8 Variation of DFE of point defects near XIwall and XIIwall ............................... 114 5 9 DFE of Li Frenkel defect pair near XIwall and XIIwall ..................................... 115 5 10 DFE of Nb antisite defect pair near XIwall and XIIwall .................................... 116 5 11 The variation of DFE of Mg dopant occupying cation sites near Y wall ............ 120 5 12 Possible arrangement for reaction 1 ................................................................. 121 5 13 The variation of DFE of defect pair using reaction 1 near Y wall ...................... 121 5 14 Possible arrangement of Mg dopants in LiNbO3 for reaction 2 ......................... 122
11 5 15 The variation of DFE of defect pai r with reaction 2 near Y wall ........................ 123 5 16 Possible arrangement of Mg dopants in LiNbO3 ............................................... 123 5 17 The variation of DFE of defect pair with reaction 5 near Y wall ........................ 124 6 1 Different shapes of domains in LiNbO3 and LiTaO3 .......................................... 127 6 2 Crystal structure of LiTaO3 for +Ps and Ps in z direction. ................................. 129 6 4 The comparison of normalized polarization between simulation and experiment ........................................................................................................ 131 6 5 Domain wal l energy comparison in LiTaO3 ....................................................... 133 6 6 Domain wall energy comparison between LiNbO3 and LiTaO3. ........................ 134 6 7 The schematic view of two possible shapes of X walls .................................... 136 6 8 Structural arrangement of cations on the Y wall plane and the X wall plane at the transition state. ........................................................................................... 137 7 1 Comparison between hybridized HF calculation and linear scaled DFT calculation. ........................................................................................................ 144 7 2 The schematics of t wo layered structure with four terthiophenes (3T) and f ive layered st ructure of thiophene oligomers .......................................................... 145 7 3 The angle between incident thiophene ring and substrate 3T ring is determined in order to study angle dependency of reaction ............................. 146 7 4 P olymerization between two 3Ts and an incident thiophene. ........................... 148 7 5 C rosslinking between two 3Ts. ......................................................................... 149 7 6 P olymerization between a 3T and an incident thiophene. ................................ 149 7 7 B roken 3T chain and an incident thiophene. .................................................... 150 7 8 The distribution of molecular weight for four different trajectories which lead successful reaction with 3Ts in the substrate .................................................. 150 7 9 The iso surface of electron density ................................................................... 152 8 1 Schematic illustration of the procedure for choosing the reaction pathway to advance the system to the next stat e in the standard KMC algorithm .............. 156 8 2 Flow chart illustrating the basic steps in the KMC algorithm. ............................ 156
12 8 3 The schematic view and flow chart of dimer saddle point searching mechanism ....................................................................................................... 158 8 4 The comparison on searching behavior of dimer method ................................. 159 8 5 The schematic view of six layered intial structure on (111) Al surface .............. 161 8 6 Two possible adsorption sites for an Al adatom. .............................................. 161 8 7 Four different diffusion processes of an Al adatom on (111) Al surface ........... 162 8 8 The evolution of time for AKMC simulation ....................................................... 163
13 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy S TRUCTURE AND DYNAMICS OF INTERF A CES IN ORGANIC AND INORGANIC MATERIALS USING ATOMIC LEVEL SIMULATION By D onghwa L ee A ugust 2010 Chair: Simon R. Phillpot Major: Materials Science and Engineering Interfaces in materials play a key role for industrial applications The structures and dynamics at various interfaces including ferroelectric domain walls, gas organic interface, organic semiconductor interface and metal gas interface are investigated with different atomic levels of simulation approaches. Ferroelec tricity: Due to their unique ferroelectric and nonlinear optical properties, trigonal ferroelectrics such as LiNbO3 and LiTaO3, are of wide interest for their potential applicati ons in optoelectronics and nonlinear optics. The properties of these materials are heavily influenced by the shape of ferroelectric domains and domain walls. Therefore, investigation of the local structure and energetics of the ferroelectric domain walls a nd their interaction with defects on atomic scales, which is not clearly understood is extremely important The structure and energetics of ferroelectric domain walls in LiNbO3 are examined using density functional theory (DFT) and molecular dynamics (M D) methods. The energetically favorable structures of 180 domain walls and the activation energy for domain wall motion are determined by atomic level simulations. The variation of
14 polarization due to the presence of domain walls is also discussed. Defects can be pinned by domain walls V arious defects domain w alls interactions and the effects on domain wall motion are described using atomic level simulation methods. Although the structure of LiTaO3 is very similar with LiNbO3, it has been said experimentally that the shapes of domain walls are different with the presence of particular defects. Using both DFT and a newly developed interatomic potential for LiTaO3, the differences in domain wall structure are understood in terms of the difference in energetics of domain walls between two materials Polymerization: Surface polymerization by ion assisted deposition (SPIAD) enables the control of thin film chemistry and morphology on the nanoscale during growth of conductive polymer thin films. This method allows fine tuning of optical band gaps and other optoelectronic properties of a polymer film by controlling the structure and kinetic energy of the depositing ions and neutrals. Thus, a comprehensive understanding of various mechanisms on the atomic level will contribute to optimizing growth conditions during SPIAD SPIAD simulations are performed to study polymerization and crosslinking behavior of polythiophene molecules at the gas organic interfaces using DFT MD method. The growth processes for polythiophene molecules are studied by depositing thiophene molecules with 25 eV kinetic energy on terthiophene surface. T he mechanism and various processes for polymerization and crosslinking of polythiophenes will be discussed The change s in bond chemistry at the polythiophene molecules and at a PbS nanocrystalline quantum dot (organic semiconductor interface) after a collision of C2H+ molecules with the substrate are also addressed.
15 Surface diffusion: Surface diffusion is a key concept for understanding catalytic behavior at the surface. We develop a new code implementing adaptive kinetic Monte Carlo ( AKMC ) method with the dimer transition searching mechanism. The code is developed with a simple LennardJones (LJ) potential A test of dimer method is performed by using 2dimensional testing potential. R esults of surface di ffusion process es of an Al adatom on Al (111) surface using AKMC method are presented.
16 CHAPTER 1 GENERAL INTRODUCTION 1.1 Imperfections I n S olids Colin Humphreys said Crystals are like people, it is the defects in them which tend to make them interesting !. As we can see from thermodynamic s, r eal crystals are never perfect : the y a lways contain defects. Defects can exist in vari eties of typ es inside materials and can have a pro f ound impact on the macroscopic properties of material s. Although some defects act as the or i gin of beautiful colors in many minerals or the conductivity in ionic conductors many natural defects degrade the properties of materials and make them difficulty for industrial usage. However, as we can control these defects, we can capture the desired properties of materials for specific applications Therefore, understanding imperfect ion s is an important problem for material scientists and engineers. Defects can be categorized by the dimension ality also ; 0 dimension (0D), 1 dimension (1D) 2 dimension (2D) and 3 dimension (3D). The size range over different classes of defects is shown in figure 1 1. Figure 11 Various classes of imperfections over the size range
17 1.1.1 P oint D efects (0 D) A crystal has a l ocalized defective region of atomic dimensions, the regio n is termed a p oint defect Different types of point defect s exist, including vacant lattice sites, and extra or missing atoms or electrons, and foreign atoms. Vacancies : A vacancy is a lattice position which is vacant because of thermal ly induced diffusi on of atoms. It can be also described to off stoichiometry. The number of vacancies formed by thermal agitation can be understood as ) exp( KT Q N NV A V (1 1) Where NA is the total number of atoms in a solid, QV is the energy to form a vacancy, K is the Boltzmann constant, and T is the temperature. Vacancies are important because they can determine the diffusion rate in the matrix. Interstitials : An interstitial atom occupies a site other than a normal lattice position. An interstitial atom can be either the same type of atom (self interstitial) or foreign atom (impurity). Impurities : An impurity atom is a non matrix atom which is not present in pure crystalline materials. Impurities atoms can be found at either vacancy or interstitial sites. B ecause they are electronically different, they would cause a difference in property, even mothered lattice structure. While vacancies and self interstitials are considered as intrinsic point defects, impurities are regarded as extrinsic point defects. P oin t defects apply strain to neighbor ing atoms cause the distortion of structure, and affect the properties of materials
18 1.1.2 L ine D efects (1 D) Line defects are defective regions of the crystal along a line. The l ine does not need to be straight; it can b e curved or even form a closed loop. The line defects are often referr ed to as dislocations ; disclinations are another type of line defects, however they are only important in liquid crystals. Because dislocation lines show abrupt changes in the regular or dering of atoms, t hey cause a change in density and their motion lead to change s in mechanical properties of material. They can be categorized into three different types depending on how the atom s displace along the dislocation line ; edge dislocation, screw dislocation and mixed dislocation. Edge dislocation: When atoms are displaced in a direction perpendicular to the dislocation line, the dislocation is termed an edge dislocation, see figure 1 2(a). Edge dislocations occur with a n extra atomic plane. Scre w dislocation : If the atom displacement is parallel to dislocation line, the dislocation is called as a screw dislocation, see figure 1 2(b). Screw dislocations result from shear displacement of atoms. Mixed dislocation: When atoms are displaced at some angle to the dislocation line it can have both edge and screw components. Figure 12 Geometrical production of dis locations  (a) (b)
19 1.1.3 Interfac ial or P lanar D efects (2 D) Interfacial or planar defects in solids can be categ orized into three different types; interfaces, domain walls a nd grain boundaries. I nterface s : Interfaces betw een solid and gas are called free surface s. Because t he environment of an atom at a surface differs from that of an atom in the bulk, the structure near the surface has a slightly larger lattice parameter tha n in the interior. The unsatisfied atomic bonds at the densest atomic density plane usually shows the lowest surface energy. The surface energy is normally in the range of a few J/m2. Ferroelectric Domain walls : A change in electronic structure, chemical composition or atomic arrangement across the interface leads to interphase boundaries In particular the inter face leading a change in the atomic arrangement but no change in the symmetry of the crystal structure, is a ferroelectric domain wall More details on ferroe lectric domain walls are given in chapter 3. Grain boundaries : Grain boundaries are planar defects which separate regions of different crystallographic orientations. A grain is a region of crystal in which the principal crystallographic directions do not change. Because of the different crystallogr a phical orientation of each grain the grain boundaries lead to atomic mismatch at the interfaces and create lower density structures. Grain boundaries are typically very reactive, so that point defects are usual ly segregate there. The interaction between interfaces and point defects affects properties of materials.
20 1.1.4 Bulk or V olume D efects (3 D) Bulk defects are generally created during production of materials such as casting, forging, and rolling. Three com mon types of bulk defects are inclusions, cracks and pores. Inclusion: An i nclusion is an unwanted secondphase particle which is pr oduced during material proces sing Inclusions can be found in a variety of shapes and sizes in crystals. Cracks : A crack is the region where material has fractures or discontinuities Cracks are normally created by either mechanical or thermal stresses dur ing forging or welding Since c racks can propagate with the external load, they affect the mechanical properties. Pores : Por es are generated during the solidification process. The void created by the volume change from liquid to solid remains as pores in a solid. Pores can affect optical, thermal and mechanical properties. 1.2 Interface S cience Interface science is the study o f physical and chemical phenomena that occur at the interface of two phases such as domain walls, organic gas interfaces, organic semiconductor interfaces and solid gas interfaces As advanced materials become increasingly multifunctional, there is a concerted focus on the influen ce of interfaces. Interfaces not only change mechanical pr operties of materials, such as elasticity, impact strength, and wear resistance but also affect electrical and optical properties, such as electrical conductivity, polarizability and refractive index. To tailor the interfacial properties of materials, it is therefore necessary to have a solid fundamental understanding of atomistic phenomena at the interface. For ferroelectric materials, two
21 oppositely polarized ferroelectric s will have internal interfaces known as ferroelectric domain walls. Therefore, the first system to be addressed in this dissertation is domain wall structures and properties in LiNbO3, which is a widely used material for microactuator s or optical filter s. This understanding will help for its potential usage for ferroelectric random access memory (FeRAM) or optical switch. We also perform a study on LiTaO3 which has very similar structure but different material properties. The difference in material proper ties of LiNbO3 and LiTaO3 are understood by atomic level simulation. Defects d omain wall interactions can change the properties of ferroelectric materials. Thus, the next study in this thesis I address the interaction between point defects ( intrinsic and extrinsic ) and domain walls in LiNbO3. The next interfacial study is the deposition kinetics on a gas polymer interface. The deposition of polymer particles into the conductive polymer chain is performed to understand polymerization and crosslinking proc ess of poly thiophene polymers. Because p olythiophene is one of the most widely studied conductive polymers for diverse electronic application areas, such as organic light emitting diode and photovoltaic devices, the deposition of thiophene molecules on the terthiophene (3T) substrate is studied. The variation of chemical boding at the gas semiconductor interface is also studied in order to understand the electronic conduction mechanism between conductive thiophene polymer chains and nanocrystalline quantum dots. The third and last area which this theses address is the diffusion process of adsorbed particle on the FCC (111) surface. For this study, a n effective tool for studying surface migration is developed. The diffusion kinetic s of an a luminum (Al) adatom on Al (111) surface is studied.
22 CHAPTER 2 SIMULATION METHODS 2 .1 Simulation Methodologies Various simulation methodologies have been developed t o elucidate materials behaviors ov er a wide length scale and time scale. Atomic simulation is a powerful me thod to investigate structure and energetics of materials O n the scale of atoms, t wo types of simulation methods are widely performed; consider ing interactions between essentially structureless atoms with empirical potential s ( atomi sti c methods) or consid er ing the interactions between atom s and electrons explicitly (electronic structure method) For the atomistic methods, the interactions between atoms are simplified by using empirical potential s. The atomistic method can be performed quickly without high computational load. However, the choice of empirical potential is critical for the accuracy. The electronic method can improve the accuracy by considering electronelectron interactions. However, the computational demands for calculating electronelectron interactions are extremely high, thus limiting its application to small system sizes for prototype study Thus, density functional theory (DFT) which regards the electronelectron interactions as interactions between electron and elect r onic density is developed. DFT improve s the efficiency of electronic methods by reducing significant amount of time and resource required, and make s it possible to perform the electronic structure calculations on several hundred atomic system. In this study, empirical and D FT methods are used to study the structure and energetics of materials. M olecular dynamics ( MD ) simulations with DFT rather than empirical methods enables a different level of dynamic studies Typically MD simulations are performed in small range of time scale s (~109~1012 s), because they predict the atomic behavior
23 from reproducing the vibration of each atom species ; the characteristic vibrational period being ~1013 s. In this study, MD simulations combined with density functional theory (DFT) are per formed in order to study dynamic process in materials including chemical reactions, polymerization, crosslinking and charge transfer at the interface. R are kinetic events happen in g at longer time scale (103 ~ 106 s) can be also explored either with temp erature accelerated dynamics (TAD) or with Adaptive Kinetic Monte Carlo ( AKMC ) methods E lectronic and atomistic scale modeling is used to understand structure, energetics and bonding chemistry of material properties in detail ; AKMC and TAD are used to pre dict the kinetic process and unknown rare events In this study, the newly implemented AKMC method into empirical potential is used in order to interpret the surface diffusion process of a particle on a metal surface. 2 2 First Principles Calculations F irs t principles calculations refer to calculation s made by considering the int eractions of electrons and nucleus in a system by solving Schrdingers equation within a set of approximations. The time dependent Schrdinger e quation for a general quantum system can be written as:  ^H t i ( 2 1) where ^H is the Hamiltonian operator, is th e wave function and t i is the ener gy operator. For a single particle in a potential V equation ( 2 1) can be rewritten as: ) ( 22 2 i ir V m t i ( 2 2)
24 where m is the mass of the particle. The first term and second term o n the right side of the equation (2 2) represent the kinetic en ergy and potential energy of the system respectively For a few specific case s such as a hydrogen atom e quation ( 2 2) can be solved. However, for more complex systems, it is impossible to obtain the analytic solution for equation ( 2 2 ) If we assume the w ave function separated from a time dependent eigen function, we can change the wave function in equation ( 2 2) to the timeindependent eigen function. The Born Oppenheimer approximation[4 ] allows the wave function of molecule to be separated into electronic and nucleus components the e quation ( 2 2 ) can be e xpanded to : E r r U r V mN i i j j i N i i N i i 1 1 1 2 2, ) ( 2 ( 2 3 ) where, m is the mass of electron and N is the total number of electrons in the system. The individual terms in the Hamiltonian in e quation ( 2 3) are, in order, the kinetic energy of each electron, the electronnuclei interactions, and the electronic electron interactions respectively. We can separate the electrons into two categories : core and valence electrons The v alence electrons represent for the electrons at the outmost shell and the core electrons term the electrons at the inner shell. Intuitively, the core electrions of an atom are only weakly perturbed by chemical reactions, geometry and bonding change s. Thus, several approximations are developed in order to reduce the computational expense by treating them only approximately The simplest level of approximation is the frozencore approximation which treats the electron occupa ncy of core states as fixed. Alternatively we can use an e ffective potential, combining opposite
25 sign charges between core electron and nucleus a so called pseudopotential. We give more detail on pseudopotential s in Sec 2 3 4 The electronelectron interaction is the most challenging and computationally expensive calculation It requires various assumptions for the solution, since the interaction of each wave function with all other electrons has to be considered. The total wave function can be denoted by a function of the spatial coordinates of the N electrons : Nr r r ..., ,......... ,2 1 ( 2 4 ) With the Hartree product the equation ( 2 4) can be approximated as a product of individual wave functions : r r rN ........2 1 ( 2 5 ) W hile t he Hartree product approximation is fairly convenient, it fails to describe the antisymmetric condition of the wavefunction of fermions. The Hartree Fock method was developed to overcome this problem. In HartreeFock calculations, the oneelectron wave functions are approximated by linear combination of atomic orbital s. Commonly, the atomic orbitals are assumed to be composed of a linear combination of Gaussian type orbital s. Various types of sets of Gaussian functions, known as basis sets, have been developed and incorporated in configuration interaction (CI), coupled cluster (CC), Moller Plesset perturbation theory (MP), and the quadratic configuration interaction (QCI) approaches  A lthough different basis sets have been suggested molecular Hartre e Fock calculations are not us ed in most modern simulations due to the high numerical costs.
26 An alternate formalism has been developed based on the probability of finding an electron at a particular system coordinate The electron density at a particular position r in space can be given as: i i ir r r n*2 ( 2 6 ) where, the asterisk indicates the complex conjugate of the wave f unction. The prefactor of 2 in e quation ( 2 6) comes from the electrons of two different spin states which can occup y each orbital. T he use of the electron density reduces the many body problems of N electrons with 3N spatial coor dinates to 3 spatial coordinates. This idea is the foundation of density f unctional theory (DFT). 2 3 Density Functional Theory (DFT) Density functional theory (DFT) is a method to investigate the ground state electronic structure using quantum mechanical theory With this theory, the properties of a many electron system are det ermined by the spatially dependent electron density DFT is among the most popular and versatile methods available in many different area including condensedmatter phy sics, computational physics and computational chemistry 2 3 .1 HohenbergKohn T heorem s Although early DFT has conceptual roots in the Thomas Fermi model, its modern form originat es from t wo different Hohenberg Kohn theorems  HohenbergKohn Existence theorem: the groundstate energy from Schrdingers Equation is a unique functional of the electron density, n(r) This means that ground state properties of a many electron system are uniquely determined by an electron density that depends on only 3 spatial coordinates. Also, they noted that the
27 integrals densi ty is th e num ber of electrons in the system. T he Hamiltonian operator and the electron density are sufficient to determine the energy Thus, DFT theorem can be extended to excited states as well as ground state. HohenbergKohn v ariational theorem : The second Hohenberg Kohn theorem is that the density obeys a variational principle. The energy, determined from Hamiltonian with any given density must be equal to or great er than the groundstate energy. So, we can minimize the energy functional by choosing different densities The energy functional is determined by classical interaction by three different terms: the attraction between the density and the nuclei, the self repulsion of a classical charge distribution, and the kinetic energy of a continuous charge distri bution. 2 .3 2 KohnSham Self C onsistent Field Methodology Based on the HohenbergKohn theorem the electron density determines the ground state originated from many electron wave function. However, we still need to include the electronelectron interactio n term of Schrdinger equation. In 1965, Kohn and Sham  suggested a simple model for solving the electronelectron interaction using a Ha miltonian operator which expressed as a hartree product of one electron operators. T he energy functional can be divided into two parts as: i XC i known iE E E ( 2 7 ) Where, ion i i i i knownE r d r d r r r n r n e r d r n r V r d m E 3 3 2 3 3 2 22 ( 2 8 ) and i XCE is all of t he non classical corrections to the electronelectron repulsion energy. The terms in e quation ( 2 8) are, from left to right, the kinetic energy of non-
28 interacting electrons, the Coulomb interaction between the electrons and nuclei, the classical electron e lectron repulsion, and the correction to the kinetic energy deriving from the interacting nature of the electrons i XCE defines all other interactions not included by e quation ( 2 8), and is referred to as the exchangecorrelation energy In order to solve for i XCE Kohn Sham used a set of e quations involving only single electrons Therefore, e quation ( 2 3) can be modified by the KohnSham one electron operator. T he KohnSham e quation can be written as: r r r V r V r V mi i i XC H ion 2 22 ( 2 9 ) where Vion (r) is the electronion potential, VH (r) is the Hartree potential, VXC (r) is the exchangecorrelation potential, i is the KohnSham eigenvalue, and i is the wave function of state i Vion (r) represents the attraction between a single electron and nuclei and can be expressed as: nuclei n ionr r Z r V') ( ( 2 10) ) ( r VH represents the repulsion of a single electron from electron density and is given by: 3 2) ( r d r r r n e r VH ( 2 1 1 ) The exchangecorrelation potential is the functional derivative of exchangecorrelation energy as the oneelectron operator and is given as: r n r E r VXC XC ( 2 1 2 )
29 where, EXC(r) is exchangecorrelation energy. Even though e quation ( 2 9 ) looks very similar to e quation ( 2 3) there is one very important difference: e quation (23) considers the interaction of two wave functions of electrons. However, the Kohn Sham e quation considers a singl e electron interacting with the electron density Thus, a set of the interactions between two electron wave function s can be replaced by electrons interacting with the electron density So increasing the number of electrons in the system only increases the number of singleparticle s for equation (29 ). This will improve the computational scales from N4 to N3 by The K ohnS ham S elf Consistent F ield procedure can be summarized as figure 2 1 Figure 2 1. Flo w chart of the KohnSham SCF procedure. 
30 2 3 3 Exchange Correlation Functional The Exchangecorrelation term, VXC(r), re p resents both the difference between the classical and quantum mechanical electronelectron repulsion and the difference in kinetic energy between the fictitious noninteracting system and the real system How ever, in most modern functional s, the energies for this portion is replaced by values obtained under specific conditions instead of solving them explicitly. The exchangecorrelation energy functional, EXC, can be expressed as an interaction between th e electron density and energy density, XC dr r n r n r n EXC XC) ( ) ( ) (* ( 2 13) w here, XC is a sum of individual exchange and correlation contributions. Although the exact exchangecorrelation functional varies DF T determines EXC[n(r)] by making r easonable approximations. Needless to say, i t is important to choose a good exchangecorrelation approximation for the system. One of t he most common approximations used in DFT calculation is the local density approximat ion (LDA). In LDA, the energy density at some position r is computed exclusively from the value of the electron density at that position. LDA is derived from the homogeneous electron gas in which the density has the same value at every position. In a molec ular system, the electron density is normally far from homogeneous so the LDA approach has the limitations. A more sophisticated approximation to improve the exc hangecorrelation functional is to consider the change in the local density using the gradient of the density. The gradient correction in the density defines the generalized gradient approximation (GGA)  The gradient c orrected functional is added to the LDA functional, i.e.,
31 ) ( ) ( ) ( ) (3 / 4r r n r n r nXC LDA XC GGA XC ( 2 1 4 ) Although GGA treat the exchangecorrelation functional in a more sophisticated manner GGA does not always provide a more accurate description of all systems  Therefore, c hoosing the most appropriate exch angecorrelation functional for the system is important. 2 .3. 4 Pseudopotentials Although we simplified the Schrodinger equation using the Kohn Sham Theorem, all electrons must be included to solve the equation. Intuitiv ely, we expect the core electrons of an atom to interact with other atoms less strongly. T herefore, we can use a pseudopotential to simpli fy the calculation. The basic concept of a pseudopotential is to replace the exact potential created by the attraction between nucleus and core electrons and repulsion between core electrons with an effective potential. Th us the pseudopotential reduces the calculation time during self consisent field method. There are four different requirements for generating pseudopotentials 1 ) Boundary matching: the angular momentum of all electron and pseudowavefunction must match at the boundary between core and valence electrions 2 ) Smoothness: Pseudowavefunction should not have a maximum w ithin the core radius 3 ) Eigenvalue matching: the eigenvalues of pseudowavefunction must match the all electron values at the chosen atomic reference state. 4 ) Norm conversation: The valence electron density at core region should be equal in the pseudopotential and in the all electron cases.
32 2 3 5 Plane w ave Implementation Although we have reduced the number of electrons for solving the Schrdinger equation of system by using Kohn Sham e quation and a pseudopotential, the calculation for many electron systems is still challenging One of the methods to enable the calculation is the planewave approach for the wave function (). A plane wave is a constant frequency wave having constant amplitude normal to the phase velocity vector. For periodic system s, such as a solid, the wave function can be described as a plane wave  which has the following form: r iK j k je r U r. .. ( 2 1 5 ) where, Uj is a lattice periodic component, 1 i and K is the plane wave vector. T he plane wave vector K is the reciprocal of wavelength and the direction of vector indicate the direction of wave propagation. Expressin g the lattice periodic component as a sum of discrete planewave basis sets with the reciprocal lattice vectors of the crystal G e quation ( 2 15) can be rewritten as: r G K i G G K j k je C r). ( .. ( 2 1 6 ) where Cj, K+ G is the planewave coefficient. The dimension of the planewave basis set s should be infinite for an exact solution. T runcating the basis set is required for computational schemes. The basis sets can be reduced by defining an energy cut off given by cutE G K 22 In practice, plane wave basis sets are often used in combination with a pseudopotential so that the plane waves are only used to describe the valence charge density. The advantage of planew ave basis sets is that all functions in the basis are mutually orthogonal, so it
33 does not show a basis set superposition error Also, derivatives are computationally less expensive, because plane wave basis sets can use Fast Fourier Transforms (FFT). Another important advantage of a planewave basis is convergence to the target wavefunction. However plane wave basis sets are less suitable for gas interface calculation because it requires increased number of calculations over the vacuum space For a periodic system, the integrals over infinitely extended system are replaced by integrals over the finite first Brillouin zone in reciprocal space. Suc h integrals are performed by summing the electron density at a finite number of points in the Brillouin zone, called the kpoint mesh T he integrations are replaced by the summations of Fourier transform at a finite number of K points: BZ j j jK F w K d K F 1 ( 2 1 7 ) wh ere, F(K) is the Fourier transform of an integral function such as electron density or total energy, is the supercell volume, and wj is the weighting factor. Choosing the k point mesh appropriately is crucial for the convergence of the results, requiring convergence tests. Since the wave functions at closely spaced K points are nearly identical, it turns out that only a few K points are typically necessary to describe the whole reciprocal space. 2 .3. 6 Linear Scaling Because plane wave DFT requires solv in g the wavefunction of Kohn sham equation via density fitting, it scales as the cube of the number of atoms (or electrons) of the system. As the system size increases, determining the eigenvalues and eigenvectors will be the major cost. This makes it ver y difficult to reach system sizes larger than a few hundreds atoms, and is therefore a huge barrier for the study of many
34 problems in nanoscale materials. In order to overcome this limitation, locality is considered in real space. If we assume that an atom only interact with other particles (electrons and nuclei) within a given radius, then the number of interactions per atom remains constant regardless of the system size. In other words, the total cost for the calculation will scale linearly with the dimension of the system. One might also think of how the central role of t he long ranged Coulombic potential in the Hamiltonian is disregarded when atoms are separated more than the cutoff distance. This is the result of screening. The Coulombic interactions are effective within quite a short range in real space, the so called near sightedness principle. Based on this idea, the linear scaling of the computational expense for large system  can be achieved. The main challenge is to determine the radius for interaction in order to obtain a physically well motivated basis function. Typically, the atomic orbitals correspond to a set of functions which decay exponentially with distance from the nuclei. Later, Boys showed that the atomic orbital s could be approximated as linear combinations of Gaussian orbitals  For linear scaling calcu lations, Gaussian basis functions lead to huge computational savings because of the ease in the calculation of overlap, and other integrals. The physically motivated shape of atomic orbital may not be the best answer for the basis set since the pseudopotential used modifies the shape of orbital in the nuclear region. Therefore, Pseudo Atomic Orbitals (PAO s) are used for the standard choice of basis set.[18,19] These PAO s can be determined during the generation of pseudopotential as a form for describing the isolated atom. The PAO s decrease rapidly with distance, reaching zero at infinite radius. L ocality in real space is imposed to achieve linear scaling. The t ail of
35 the PAO is modified to reach zero at a given radius by c onfining the eigenfunctions of the pseudo atomic problems within a spherical boundary at which the potential becomes infinite  The radial confinement is an approximation and can be chosen to achieve an appropriate compromise between high precision and great computational efficiency. Whil e PAOs provide a single basis set, multiplezeta basis sets needs to be consider ed for increased variational freedom such as the change in chemical bonding, external fields or other perturbations to electronic structure. The use of multiple zeta basis sets multiple Gaussian functions, allow the effective atom size to respond to its environment. 2 .3. 7 Application to LiNbO3 S tudy For ferroelectric LiNbO3, the Vienna Ab Initio Simulation Package (VASP) , a plane wave DFT package, is used to the study the properties of domain walls. VASP offers the two most widely used pseudopotentials ; ultrasoftpseudopotentials (USPP)  and projector augmented wave (PAW)  For USPP, the norm conservation is relaxed, and then an augmentation charge density is added. The PAW approach focuses on the augmentation of the wavefunction instead of the density, so all electron properties are captured in the frozen core region. Therefore, PAW methods leads to a better description of the core electrons than USPP i n general PAW requires smaller core radii and a larger energy cutoff than USPP A compar ison between USPP and PAW method shows that well constructed USPP and PAW method give almost identical results and are in good agreement with hyper precise allelectron calculations. However, the PAW method yields more reliable results than the USPP if the system includes the materials having strong electric and magnetic moments For ferro electric LiNbO3 system, the PAW method yield better results. Therefore, PAW basis set is used for investigating the ferroelectricity of LiNbO3.  The projected augmented wave
36 (PAW)  pseudopotential within the generalized gradient approximation (GGA)  is used to evaluate the exchange and correlation interactions. The outmost shells of each ion, 2s1 on Li, 4p64d45s1 on Nb and 2s22p4 on O, are considered as active valence states for the interactions and are treated explicitly; the ener gy cut off for the plane waves is 400 eV.  A conjugate gradient and quasi Newton algorithm are used for the ionic relaxation.  The force criterion for complete ionic relaxation is 0.0001eV/. The residual minimization method direct inversion in the iterative subspace (RMM DIIS) algorithm  which optimizes several individual energy bands at the same time, is used for electronic energy minimization. The pseudopotential and methodologies used here are the same as used previously in studies of intrinsic defects in LiNbO3.  2 .3. 8 Application to Molecular B eam D eposition For studies in volving a gas phase such as molecular beam deposition, plane wave DFT is not computationally efficient enough to carry out the simulation. Instead, linear scaling DFT is a much more powerful method to overcome this limitation. VASP is not suitable because it uses plane wave based pseudopotentials and basis sets. Therefore, we use Spanish Initiative for Electronic Simulations with Thousands of Atoms ( SIESTA) software package for this study Also, B3LYP hybrid method that is based on the mixing of quantum chemical and DFT methods is used for the validation of SIESTA. When plane wave basis set s are implemented into pseudopotent ial s for the core region, such as PAW or USPP the number of calculation s can be decreased with the cut off distance of the pseudopotential s in reciprocal space. However, implementing plane wave basis sets of core electrons into the pseudopotential s of SIE STA do not give much benefit because SIESTA works with real space localized basis function. Therefore, norm conserving pseudopotential s are implemented in SIESTA. While t he
37 pseudopotential s in SIESTA reproduc e exchangecorrelation interactions between nuc lei and core electrons by using either LDA or GGA Gaussian type basis set s are implemented explicitly during the calculation of atomic orbital. For this study, t he nucleus and core electrons of each atom are represented with a Troullier and Martins  type norm conserving pseudopotential through the GGA method The valence electron KohnSham states are expanded in a linear combination of atomic orbital type basis sets.  A doublezeta polarized basis is employed and the valence basis set is localized within a soft confining potential in order to prevent the discontinuities in second derivative of atomic orbital. The basis set parameters for each atom species are optimized to their relative energies and compared to high level quantum mechanical calculations, unrestricted coupled cluster spin contamination corrected with correlationconsistent polarized valence triple zeta basis set ( u CCSD(t)/cc pVTZ ) [27 29] The optimized basis set is polarized with an additional d orbital on carbon and sulfur and an additional p orbital on hydrogen. The soft confinement potential is set to 2 0 Ry ( 272 eV) with an inner radius of 0.95 times the outer hard cut off. A split norm fraction of 0. 2 is used in the double zeta construction of PAO A mesh cutoff of 1 50 Ry (224 0 eV) is used to determine the fineness of the auxiliary grid basis for the expansion of the density. 2 .3. 9 C harge analysis M ethods The output of quantum mechanical calculations is the electronic charge density. A tomic charges in molecules or solids are not direct observables in quantum mechanical theory. Because t he electronic properties of material s can be fundamentally understood by the charge distribution between components, u nderstanding the variation of charge helps to trace the variation of properties during reaction. There are various
38 schemes to separate electrons from small portions of the system such as atoms or molecules including Mulliken analysis  n atural bond orbital analysis  Lwdin analysis  and Bader analysis  Two different charge analysis schemes are implemented for this study : Mulliken and Bader [34 36] analysis the two most common schemes in order to partition electrons between atoms or molecules in the system Mulliken charge analysis: Mulliken charge analysis is based on the direct population of the orbital s. This scheme can be applied when basis functions of atomic orbital are used in the form of electronic wavefunction of the system. The Mulliken charge on any atom is calculated by the charge associated with the basis functions centered on the particular atom. Even though this method is fast and useful for determining partial charges on atoms it is sensitive to the choice of basis set. So, Mulliken charge is very dependent on the number of zetas and their cutoffs. Also, t he charge assigned to an atom becomes arbitrary for an infinite basis. For plane wave basis functions, which are not associated with any particular atom, Mulliken charge analysis is not applicable. Therefore, Mulliken charge analysis is mainly used in the Gaussian and SIESTA calculation which uses Gaussian type atomic orbital basis set. Bader charge analysis : Th e scheme proposed by Bader is based on the electronic charge de nsity instead of the electronic orbitals For Bader charge, 3D space is divided into subsystems each of which normally contains one nucleus. Then, the critical points where t he gradient of charge density is zero are determined by constructing the zerof lux surface. Sets of gradients of electron density are traced in an infinitesimal step until each gradient reaches a maximum at a nucleus. The electron density is integrated within the region between the nucleus and zero flux critical point.
39 Since this sch eme is based on the charge density, it can obtain partial charges and dipole moments of individual atoms in molecules or crystals. Also, Bader analysis is suitable for plane wave basis calculations because it does not require information on the electronic orbitals However, Bader analysis has convergence problems in some systems. This scheme is used for determining the charge of plane wave basis set implemented in VASP. 2 .3. 10 Limitations on DFT Despite the improvements in DFT, it still has limitations. T hes e limitation s aris e from two major sources First, dominati on of the Coulomb term pushes electrons apart Second, describ ing the interaction of degenerate states u sing the electron density i s difficult The se leads to two major limitations on DFT including delocalization error and static correlation error  Th e delocalization error refers to the tendency of the approximate functional to spread out the electron density. When an electron is delocalized over two center s, DFT predict s a lower energy for overly dispersed fractional charges This behavior can cause o verestimation of the charge distribution during chemical reaction, since the electrons delocalized over more than one center.  Thus DFT predict s highly polariz ed system s because fractional charges appear at the edges of the molecule.  DFT can also predict unphysical charge transfer between the molecule and the metal in molecular electronic device and cause an overestimation of electronic conductance. [40, 41] Thus the delocalization error explains the underestimation of the band gap for DFT calculation.  The static correlation error the overestimation of the d issociation energy arise s in situation s of degeneracy or near degeneracy of electron spin states in strongly
40 correlated systems or in transition metal chemistry DFT do es not accurately describe the interaction between the degenerate spin states and thus gives a large error in des cribing fractional spin states. The fractional spins cannot be explained by a static correlation function based on constancy of electron correlation. So, this leads to a difficulty in using the electron density to describe degenerate states in transition m etal s.  As a result, density functional theory is known t o be weak at describing intermolecular interactions especially van der Waals forces (dispersion) charge transfer excitations transition states, global potential energy surfaces and some other strongly correlated systems such as the band gap in semiconductors Th e se known problems make DFT unsuitable for the systems which are dominated by dispersion including interaction with noble gas atoms or biomolecules Although there is still much work to do, diverse corrections [42 44] have been suggested to the functional or the inclusion of additive terms to improve many related predictions 2 4 Molecular Dynamics (MD) Simulation 2 4 .1 Introduction Electr onic structure calculations generally give the highest available materials fidelity. Unfortunately they are not suitable for the systematic study of system such as domain walls because the computational load for the system sizes required of systems is prohibitively high. Mo lecular dynamics (MD) is an alternative method for studying molecular behavior. MD simulation is a specialized discipline for predicting the motion of atoms based on statistical mechanics. MD simulation traces the atomic behavior by solvi ng Newtons e quation of motion. Newtons second law can be expressed by:
41 i i im a F ( 2 1 8 ) where, Fi is the force on an atom i mi is mass and ai is the acceleration of atom i respectively. The force on atom i can be expressed by the gradient of the interatomic potential energy with respect to position. i i i iz V k y V j x V i V F (2 19) where, V is the potential energy. The potential energy is calculated by different interatomic interactions in the simulation system. 2 4 2 Interatomic Interactions For ionic systems, the interatomic interactions can be described by two different terms; a long range and short range contribution. The detail s of the interactomic interac tions are discussed below. Longrange interactions: Ch arged ion s interact with another charged ion t hrough Coulombic interactions in a long range. It is described as : ij j i Coulomb ijr q q r V04 ( 2 20) where ijr is the distance betw een the position of two species whose charges are qi and qj. Because this equation is proportional to 1 ijr the potential decrease rapidly when the dis tance between ions increases. T wo different ways have been developed to accomplish the convergence of the electrostatic energy and forces within a certain cutoff distance : Ewald summation and direct summation The Ewald summation method is the traditional approach to solve the electrostatic interaction between ions in an ionic cr ystal system. However, its complexity and poor compatibility with non charge
42 neutral system  cause d limitation s f or applications. In this study, therefore, t he direct summation method of Wolf et al.  is used for the calculati on of the Coulombic interactions Ewald sum : The Ewald sum is a method by which the Coulombic interactions at a certain distance, Rc can be truncate d without causing a convergence problem. Ewald sum assumes that each point charge is surrounded by a charge dis tribution of equal magnitude and opposite sign. This imaginary charge is normally considered as Gaussian type distribution. T he total sum of these charge can be divided into two different terms : screening interaction between the neighboring charges and charge distribution of the lattice The screening interaction between charges can be summed in real space, while the long range charge distribution term is calculated in Fourier space, so call reciprocal space. The r eal space summation can be easily achiev ed. However, the summation of longrange charge distribution in reciprocal space is less intuitive and requires an assumption that the system is infinitely periodic. A cut off distance Rc is used to improve the simulation efficiency without significant lo ss of accuracy. Normal ly a cut off of 2 or 3 times the lattice parameter is enough to obtain a good description of long range interactions However, the computational load in the Ewald summation method increases as N2 for N ions. Moreover, t he method is n ot physically transparent. Direct sum: Although the E wald summation method is more precise, the computational load is considerable and it causes convergence issue with cut off distance, Rc. Thus, t he direct summation was developed by Wolf et al. and is w idely used to calculate the long range Coulombic interatomic interactions The convergence issue with charge neutrality is resolved by placing neutralizing charges on the surface of
43 truncation sphere, see figure 2 2. Thus, any net charge contained in a sph er ically truncated system is compensated at cutoff distance, Rc. Figure 2 2. 2D Schematic of the charge neutralization scheme used in the direct sum method. For any atom i, interacting with atom j, an equal and opposite charge is placed at the perimeter of the circular cut off to make the system charge neutral. (Reproduced from reference  ) T he energy for the dir ect sum can be correctly calculated by subtracting chargeneutraliz ation term from Coulombic sum between interatomic charges under cutoff distance, Rc. The direct sum can be represented as: N R r i j c c i i ij j i c Coulc ijR R q q r q q R E ( 2 21) w here, qi( Rc) is the net charge within the cutoff sphere of ion i The calculated energy oscillates significantly, in a slightly damped manner. Thus da mping the Coulombic potential can be used to quickly damp out the oscillations with increasing Rc. The direct summation method has su ccessfully predicted the long range Coulombic energy for disordered or noncharge neutral systems such as interfaces or
44 defects as well as bulk perfect crystal. Therefore, the direct summation method is mainly used to calculate the long range interatomic interactions in this study. Short range interactions: The short range interatomic interaction describes the repulsion between ions at a short distance. There are many different models for describing short range interatomic interactions including Buckingham  LennardJones  Morse  Inverse Gaussian TangToennes  and Rydberg  In this study, the Buckingham interatomic potential whos e ability to describe ionic crystals is well established is used. The functional form for Buckingham potential is: 6expij ij ij ij ij ij Buckr C r A r V ( 2 2 2 ) where rij is the distance between two different ion species i and j and Aij, ij and Cij are potential parameters. Th e exponential term represents the increase in repulsive energy with the decrease in the distance between atoms i and j The second term expresses the Van der Walls attractive contribution due to the dipole interactions. 2 4 3 Shell Model for Describing P olarizability of E lectrons In addition to the interatomic interaction potential, a shell model, which describes each ion as a core consisting of the nucleus and inner electrons and a shell of valence electron, is used for capturing the electro nic polarizability of the ions. The core and shell carry partial charges, the sum of which is close to the full ionic charge of each atom species. The positions of the cores and shells are determined in a fully dynamic manner with fictional masses of 10% of the mass of the individual ions assigned to the shells. The core and shell of each ion interact with the cores and shells of other ions via long range Coulombic interactions. The short range Buckingham interactions describe t he
45 shell shell interactions between the ions T he core and shell of each ion are connected by a harmonic spring: 2 22 1 k V ( 2 2 3 ) where is core shell displacement and k2 is the harmonic spring constants. 2. 4 .4 Three body I nteractions Another term which is considered for this study is a threebody interaction potential. The threebody potential represents the repulsion between bond pairs. Thus, the harmonic form is normally chosen in order to penalize the deviation from the expected angle for the coordination envi ronment. The functional form for harmonic threebody interaction is described: 2 0 2) ( 2 1 k Vijk ( 2 2 4 ) where 2k is the threebody harmonic constant and is the angle between three ionic species and 0 is the expected angle for the coordinate system. 2 4 5 Potential for LiNbO3 study The interatomic potential used for LiNbO3 system was developed by Jackson and coworkers.  T hree different interatomic interactions, Li O, Nb O and O O are included in the Buckingham potentials. The interatomic interactions are cut off using a shifted force method at a distance of 10.2, which corresponds to t wice the lattice parameter along . The coreshell model is used for oxygen ions with the shell charge. Full details are given in Table 21.
46 Table 2 1. The potential parameters for LiNbO3. LiNbO 3 Interaction A (eV) () C (eV/ 6 ) Li O 950.0 0.2610 0.0 0 Nb(Ta) O 1425.0 0.3650 0.0 0 O O 22764.0 0.1490 27.88 LiNbO 3 Shell parameter Shell charge, Y Spring constant, K 2 (eV/ 2 ) O 2 2.9 70 LiNbO3 Three body parameter Force constant, k ( eV/rad) Equilibrium angl e, 0 (deg) O Nb(Ta) O 0.5776 90 2. 4 .6 Integration Algorithms The force of an each atom can be obtained from t he interatomic potential s using e quation ( 2 19) The velocity and position of atoms at a time t, can be traced by the atomic force. The evolution of the atoms i s achieved by integrating the position, velocity and accelerations of the atoms for small increment s of time ( ) There are various integration algorithms available including Verlet  leapfrog  Beemans  and Predictor Corrector  and symplectic integrators  Each method has its advantages and disadvantages with regards to accuracy and computational load depending on the condition of system. In this study, t he Verlet algorithm  and a 5thorder Gear predictor corrector method  are used. Verlet algorithm : The v erlet algorithm is an integration scheme which offers very good stability and time reversibility for MD stud ies It reduces the level of errors using the positio n of previous steps instead of velocity. The position at the next time step is predicted by the difference in position at the previous and current time steps. This can be derived from two Taylor expansions of the position vector at two different times.
47 t d t t c t t b t t a t t t v t x t t x t d t t c t t b t t a t t t v t x t t x 5 4 3 2 5 4 3 2120 1 24 1 6 1 2 1 120 1 24 1 6 1 2 1 (2 25) Where x is the position, v is the velocity, a is the acceleration and b c and d are third, fourt h and fifth derivatives of the position. T he higher term of Taylor expansion can be neglected for small time step, adding two Taylor expansions of e quation ( 2 25) gives: t a t t a t t t x t x t t x 4 22 (2 26) Equation ( 2 26) can be rearranged as: 2 2 22 t t t x t x t t x t t x t a (2 27) Therefore, the position at next time step is calculated by positions at previous and current step and the acceleration from the interatomic potential This algorithm is implemented by SIESTA during deposition study. Predictor Corr ector algorithm: Another algorithm used for this study is t he predictor corrector algorithm Th is method uses a Taylor expansion series to predict the positions, velocities, accelerations, and higher order derivatives of position at time t + Different order of predictor corrector algorithm can be used. The fifth order Gear predictor corrector integration method is implemented in our inhouse MD code. The method estimates the position of next step in two steps. In first step, the Taylor exp ansion of positions at time t + is determined to the fifth order term derivative term. Then, the higher order derivatives of it are also list ed in a similar manner, see e quation ( 2 28)
48 t d t t d t d t tc t t c t d t t c t t b t t b t d t t c tt b t t a t t a t d t t c t t b t t a t t vt t v t d t t c t t b t t a t t v t t r tt rp p p p p p 23 2 4 3 2 5 4 3 22 1 6 1 2 1 24 1 6 1 2 1 120 124 1 6 1 2 1 (2 28) where r v and a are the position, velocity and acceleration of each atom at time t + respectively. b c and d are third, fourth and fifth derivatives of the position at the incremented time. The superscript p represents that they are the predicted values. The position, velocity, acceleration and higher order derivatives at time t + can be predicted using interatomic potentials. In the second step, th e predicted acceleration values are corrected. The correct accelerations are calculated by t he forces on the atoms. Using this force, the correct accelerations ( t t ac ) can be obtained. Therefore, the error between the predict ed accelerati on t t ap and the correct acceleration t t ac can be given as : t t a t t a t t ap c ( 2 2 9 ) T he predicted position, velocity and acceleration can be corrected using the error. t t a c t t d t t d t t a c t t c t t c t t a c t t b t t b t t a c t t a t t a t t a c t t v t t v t t a c t t r t t rp c p c p c p c p c p c 5 4 3 2 1 0 ( 2 30)
49 where, c0, c1, c2, c3, c4 and c5 are empirically determined numerical constants, respectively. The values for a fifth order predictor corrector are c0 = 3/16, c1 = 251/360, c2 = 1, c3 = 11/18, c4 = 1/6, and c5 = 1/60 [57, 59] These corrected values are used to predict the positions and other higher order derivatives at the next step Th e t crucial factor for sampling the simulation trajectory. The average phonon period for a typical atomic vibration is about 200 fs. A typical time step for capturing correct atomic vibration is 1 fs For this study, however, t he time step 0.05 fs is used to capture internal vibrations of ions between core and shell 2. 4 .7 Geometry O ptimization Figure 2 3 shows the potential energy surface of an arbitrary system. For the potential energy surface, point A represents the global minima and point B represent s a local maxima or transition state and point C represent the local minima. T he task f or geometry optimization is to find the global minima for the given structure. Although the optimized geometry is desired, the structure can be easily trapped by local minima of potential energy T here are various methods  available to find the local minimum including the bisection, golden search, steepest descent, conjugate gradient and NewtonRaphson method. In this study, steepest descent, conjug ate gradi ent and NewtonRaphson method are mainly used to find out optimized geometry structure.
50 Figure 2 3. Illustation of potential energy surface of a arbitrary function. Pict ure taken with permission from http://www.chm.bris.ac.uk/pt/harvey/msci_pract/pics/pot_surface.jpeg The s teepest descent is a first order method which uses the first derivative of a function. The first derivative goes to zero at either a minima or maxima. Figure 24 shows a schematic view of ste epest descent optimization method. the optimization scheme moves to a minimum along the negative gradient in the steepest direction. If the energy of system increases, line search among the last three steps will be performed to determine where to go for the next step. Steepest descen t has the advantage of relaxing a poor initial geometry quickly and is very useful for preliminary optimization. However, it is very slow to achieve optimization because it can oscillate around the minimum energy path.
51 Figure 24. Schematic of steepest descent method. Further improvement to steepest descent s resulted in the conjugate gradient method The conjugate gradient method follow s a modified direction determined by gradients at the current and previous searching directions instead of following the current gradient direction. The conjugate gradient method shows much better convergence than the steepest descent method. However, it has a higher computational load for each step because it needs to store the trajectory of previous st ep and also use the second derivatives. The conjugate gradient method can be expressed as: 1 i i i idg d ( 2 31) where i is determined by gradients at the current and last few points and id is the search direction at step i The NewtonR aphson method is a secondorder method in which a function at a given point is expressed as Taylor expanision The displacement directions for the next step follow the eigenvectors of the current Hessian. The step size is determined by the eigenvector. The functional form can be express as: g H H g x x1 0/ ( 2 32) and i i if x/' ( 2 33)
52 where H is the Hessian, if is the projection of gradient along the i th eigenvector whose eigenvalue is i. This method converges very efficiently when the initial point is close to a minimum point. step size. These will prevent the system from moving out of the potential energy surface by taking too large step size. The Hessian is computationally very expensive. So, approximated Hessian s, updated by energy and gradient can be also used. However, this method is not practical for optimizing large simulation system. Also, it can be easily trapped by local minima if the initial point is far away from the global minim um Due to the different properties of each method, a combination of different optimization methods is implemented in this study. 2. 4 .8 Thermodynamic Conditions ( Ensembles) The ensemble describes the specific thermodynamic conditions which can replicate the environmental condition. T he environmental condition of the system can be described by three different thermodynamic state variables : pressure (P) volume (V) and temperature (T) The most common experimental condition is constant pressure and temperature. For simulation sys tems wi th fixed number of atoms the experimental condition can be normally reproduced by the isothermal isobaric ( NPT ) ensemble, in which P and T are fixed. Other ensembles are also available for simulation study depending on the thermodynamic conditions. When t he syst em is under constant volume instead, the system is in the canonical (NVT) ensemble. Similarly, for constant pressure and enthalpy the condition for the system is isobaric isoenthalpic (NPH) ensemble In this study, the most of MD simulations are performed under NPT ensemble because it corresponds to the most experimental condition. There are various
53 algorithms to maintain the constant pressure and temperature. In this study, the constant pressure of the system is obtained by Parrinello and Rahman method [61, 62] This method originates fro m an extended Lagrangian method  which considers a set of pistons attached to the system. The variation of the force on the pistons is compensated by the variation of simulation volume to maintain a constant pressure. This method was extended by P arrinello and R ahman in order to allow it to change the shape and size of simulation cell. A t hermostat is used to obtain the constant temperature of the system. There are also different thermostats including velocity rescaling  Nos Hoover [64, 65] Berendsen  and generalized Langevin [67, 68] Velocity rescaling is the most intuitive and simplest form of thermostat. T he velocities of the atoms are rescaled to match the desired temperature. This can be given by: i d i sT T v v ( 2 34) where vs is the rescaled velocity, vi is the initial velocity, Td is target temperature and Ti is the initial temperature. In this study, the velocity rescaling method is used in the MD studies
54 CHAPTER 3 STRUCTURAL ANALYSIS In this chapter, the crystal structures of the materials under investigation are introduced. In this thesis, the primary focus is on elucidating interfacial phenomenon of the trigonal ferroelectrics, LiNbO3 and LiTaO3 and the conductive polymer, polythiophene. For better understanding of interfaces, the crystallography of trigonal ferroelectrics and polythiophene are introduced in details. 3 .1 Trigonal Ferroelectrics 3.1.1 Introduction on LiNbO3 LiNbO3 has a high spontaneous polarization, 70 2, and high Curie temperature, ~1480 K [69, 70] Recent successes using the vapor transport equilibration (VTE) method  have enabled the growth of stoichiometric LiNbO3; however, congruent LiNbO3 is Li deficient [72, 73] Although various defect models including pseudoSchottky defects, Schottky defects, and Frenkel defects have been suggested only a Nb antisite, compensated by four Li ion vacancies ( NbLi +4 VLi or AntisiteI) [74 76] and five Nb antisites compensated by four Nb vacancies (5NbLi +4 VNb or AntisiteII)  were able to successfully explain the stoichiometry of congruent LiNbO3. Atomistic modeling by Donnerberg et al.  showed that the compensation by Nb vacancies is energetically less favorable than the compensation by Li vacancies. Recent DFT calculations  showed that AntisiteI cluster and Li Frenkel pair are the most energetically favorable defect pairs in LiNbO3 under Li deficient (congruent) and Li rich (stoichiometric) conditions respectively. An analysis of AntisiteI clusters determined that several arrangements of Li vacancies around a Nb antisite are nearly energetically equivalent.  The ability of calculations with empirical potentials to describe point
55 defect properties was demonstrated by Araujo and coworkers [7 9] Moreover, two different types of domain walls, parallel and perpendicular to the c glide plane in LiNbO3 were suggested by Gopalan et al.  T he structure and energetics of them are characterized here  Before it is possible to understand the domain wall structure, it is necessary to understand the crystallography of trigonal ferroelectrics in some detail. In this section, LiNbO3 is only considered. The structure of LiTaO3 is identical. 3.1. 2 Crystallography of LiNbO3 Figure 3 1. Crystallographic direction of rhombohedral axes (middle) and hexagonal axes (bottom) in a unit cell of LiNbO3. The Cartesian axes are given at the bottom left of the figure.
56 To understand the crystallography and prop erties of LiNbO3, it is necessary to use three distinct coordinate systems: rhombohedral, hexagonal, and Cartesian coordinates, see figure 3 1. The Cartesian coordinates are denoted as x, y, z; the rhombohedral coordinates are denoted as a, b, c; the hexagonal coordinates are denoted as Y1, Y2, Y3, Z. The rhombohedral coordinates can be represented as 3 / 1 2 / 1 2 / 3 z y x a z y x b3 / 1 2 / 1 2 / 3 and z y c 3 / 1 The hexagonal coordinates can be represented as 2 / 1 2 / 3 2 / 1 2 / 3 ,3 2 1 x Y y x Y y Y Z=z. In the paraelectric R 3 c (#167) state, the Li and Nb ions occupy 6a and 6b Wyckoff positions respectively, while the oxygen ions occupy 18e Wyckoff positions, see Table 3 1. While the 6a and 6b site positions are fully defined by the crystallographic tables  the positions of the 18e oxygen sites are in addition characterized by a positional parameter XO, which defines the displacements of each oxygen ion along its associated hexagonal axis, the precise value of XO being different in different R 3 c systems. Neutron powder diffraction studies have determined the positional parameter for LiNbO3 to be XO=0.0591 0.0607  Previous DFT studies yielded XO=0.036  and 0.049 in the local density approximation (LDA) calculations with Perdew Wang 91 (PW91) and Perdew Burke Ernzerhof 96 (PBE) basis sets and XO=0.048 in a gen eralized gradient approximation ( GGA ) calculation with PBE basis set Our DFT study using the GGA pseudopotential with projector augmentedwave (PAW) basis sets gives XO=0.039 which lies within the same range. In agreement with a previous atomic level simu lation with different interatomic potential study  our simulations with empirical potentials yield XO=0.034, which is somewhat smaller than the experimental values, but very similar to the DFT values.
57 Table 3 1 Wyckoff positions for R 3 c (space group 167) state  Multipli city Hexagonal axes Coordinates 18e (Oxygen) x O 0 0 x O x O x O 0 x O 1/2 + 1/4 x O 0 1/2 + 1/4 x O x O 1/2 + 1/4 2/3 + x O 1/3 1/3 + 1/4 2/3 0 1/3 + x O 1/3 + 1/4 2/3 x O 1/3 x O 1/3 + 1/4 2/3 0 1/3 x O 5/6 + 1/4 2/3 x O 1/3 5/6 + 1/4 2/3 + x O 1/3 + x O 5/6 + 1/4 1/3 + x O 2 /3 2/3 + 1/4 1/3 0 2/3 + xO 2/3 + 1/4 1/3 x O 2/3 x O 2/3 + 1/4 1/3 0 2/3 x O 1/6 + 1/4 1/3 x O 2/3 1/6 + 1/4 1/3 + x O 2/3 + x O 1/6 + 1/4 6b (Nb,Ta) 0 0 0 0 0 1/2 2/3 1/3 1/3 2/3 1/3 5/6 1/3 2/3 2/3 1/3 2/3 1/6 6a (Li) 0 0 0 + 1/4 0 0 1/2 + 1/4 2/3 1/3 1/3 + 1/4 2/3 1/3 5/6 + 1/4 1/3 2/3 2/3 + 1/4 1/3 2/3 1/6 + 1/4 The crystallography of the ferroelectric phase is more complicated. On losing inversion symmetry along the Z direction during the phase transition t o the ferroelectric state, the space group drops to R3c (#161). The up polarization is obtained when the Nb atoms shift up from the octahedra centers and the Li atoms sits above their oxygen
58 planes, see figure 3 2. Correspondingly, down polarization is obtained when the Nb atoms are shifted down and the Li atoms sit below the oxygen planes. Figure 3 2. Schematic of the paraelectric (center) and ferroelectric phase with up polarization (right) and down polarization (left). The ferroelectric distortions as sociated with the displacement of the Nb ions (green) from the center of the In this structure, the Li and Nb ions both sit at 6a Wyckoff positions; the oxygen ions occupy the 18b positions, see Table 3 2. Positional difference between cations comes from the ZLi and ZNb(Ta) which determine the position in Z direction. For LiNbO3 and LiTaO3 systems, ZNb(Ta) is defined as 0 for simplicity of analysis X ray diffraction studies  yield a positional parameter for Li of ZLi=0.283 for up polarization state and ZLi=0.217 for down polarization states, which correspond to displacements of Z =0.033 from the positions in the paraelectr ic phase; the DFT and empirical simulations yield
59 almost identical values. In characterizing the oxygen displacements, it is important to distinguish between the Cartesian coordinates and the hexagonal coordinates, the latter being more useful for understa nding the inplane displacements. In hexagonal coordinates, the crystallographic positions of the oxygen ions in the ferroelectric state are characterized as ( u 1/3+ v, 1/12 w ) for up polarization state and ( u v, 1/3 v, 1/12+ w ) for down polarization An Xray study [87, 89] determined (u, v, w) = (0.0492, 0.0113, 0.0186), see Table 3 3 Values obtained from a more detailed crystallographic analysis  from DFT  and from simulations with the same empirical potential as we are using  yielded similar magnitudes. Our DFT and empirical calculations yield results consistent with the published val ues. The difference in oxygen position between v, ). Because the X and Y directions are at 120 to each other, the two polarization states actually have the same Cartesian x coordinate, but different y coordinate. Bec ause of the threefold rotation symmetry, the oxygen ions move along three equivalent Y directions, one of which is aligned to be parallel to the Cartesian y axis. Because it is easier to understand the structure in the rhombohedral unit cell, figure 3 3 shows the structure along the rhombohedral  direction (corresponding to the [42 6 1 ] direction in the hexagonal coordinate system). In this direction, and in the  and  directions, the alternating oxygen octahedra tilt in opposite directions plane displacement tilting can be expressed as a-a-awith =24.1 
60 Table 3 2 Wyckoff positions for R3c (space group 161) state  Multipli city Hexagonal axes Coordinates 18b (Oxygen) x O y O z O y O x O y O z O y O x O x O z O y O x O 1/2 + z O y O x O y O 1/2 + z O x O x O y O 1/2 + z O 2/3 + x O 1/3 + y O 1/3 + z O 2/3 y O 1/3 + x O y O 1/3 + z O y O +2/3 x O 1/3 x O 1/3 + z O 2/3 y O 1/3 x O 5/6 + z O y O +2/3 x O 1/3 + y O 5/6 + z O 2/3 + x O 1/3 + x O y O 5/6 + z O 1/3 + x O 2/3 + y O 2/3 + z O 1/3 yO 2/3 + xO yO 2/3 + zO y O +1/3 x O 2/3 x O 2/3 + z O 1/3 y O 2/3 x O 1/6 + z O y O +1/ 3 x O 2/3 + y O 1/6 + z O 1/3 + x O 2/3 + x O y O 1/6 + z O 6a (Nb,Ta) 0 0 0 + z Nb(Ta) 0 0 1/2 + z Nb(Ta) 2/3 1/3 1/3 + z Nb(Ta) 2/3 1/3 5/6 + z Nb(Ta) 1/3 2/3 2/3 + z Nb(Ta) 1/3 2/3 1/6 + z Nb(Ta) 6a (Li) 0 0 0 + z Li 0 0 1/2 + z Li 2/3 1/3 1/3 + z Li 2/3 1/3 5/6 + z Li 1/3 2/3 2/3 + z Li 1/3 2/3 1/6 + z Li As a result of this tilting in three spatial directions, there are six distinct oxygen layers along the z direction For the up polarization state, the three bottom oxygen layers rotate counter clockwise about the z axis while the three top oxygen layers rotate clockwise, see figure 3 4(a); this six oxygen layer pattern repeats. The rotation angle
61 between the first and second oxygen layers is slightly different from the angle between the second and third oxygen layers. The first oxygen layer is rotated by 7.7 clockwise direction along z direction from the nondisplaced center position. Then, the second oxygen layer is rotated 63.6 shows slightly lager rotation angle, 67.8 clockwise direction. The three upper oxygen layers rotate in the opposite direction. The first oxygen layer shows rotation angle, 7.7 clockwise direction from nondistorted lattice. Then, the second and third oxygen planes are each rotated by 63.6 and 67.8 phase rotations cancel in both th e x and y directions, leading to a net zero inplane polarization. The down polarization shows corresponding displacements. Although the rotation angles between each layer are different from those in the up polarization state, the cancelation of the clockwise rotation of the upper three oxygen layers and the counter clockwise rotation of the bottom three oxygen layers again leads to net zero inplane polarization. Figure 33.  projection of the rhombohedral structure. Alternating oxygen octahedra tilt in op 24.1
62 Table 3 3. Comparison of positional parameters of LiNbO3 between empirical calculations, DFT calculations and experimental data. For the empirical study, the two different values correspond to 0 K (upper) and at 293 K (lower; in the parenthesis). The 293 K values represent the average values from molecular dynamics simulation. It is not able to stabilize the paraelectric phase at 293K using molecular dynamics simulation. Lattice parameter X O,para Z Li,ferr o u V w a () c () Experiment [82, 87, 89] 5.148 13.863 0.0591~0.0607 0.283 0.0492 0.0113 0.0186 LDA(PW)  5.067 13.721 0.036~0.049 0.285 0.0427 0.0125 0.0183 GGA(PBE)  5.200 13.873 0.048 0.282 0.0479 0.0097 0.0199 GGA(PAW) (present) 5.132 13.884 0.039 0.282 0.0388 0.0128 0.0207 Empirical  5.156 (5.187) 13.683 (13.710) 0.034 (0.053) 0.352 (0.288) 0.0525 (0.0563) 0.0310 (0.0294) 0.0495 (0.0133) Empirical (present) 5.169 (5.185) 13.685 (13.738) 0.034 0.283 (0.281) 0.0527 (0.0504) 0.0308 (0.0313) 0.0208 (0.0 207) Because the inplane positions of the oxygen ions differ in the up and down directions, when the polarization is flipped from up to down, the positions of the oxygen ions change along one of the three equivalent hexagonal directions. Thus, the three different sets of 6 oxygen ions can be denoted as O1, O2 and O3 and their corresponding hexagonal axes as Y1, Y2 and Y3. The arrows in Figs. 3 4 (b) and (c) represent the directions of the displacements of the three types of oxygen ions when the system swi tches from up to down polarization. The domain walls, parallel to the crystallogr aphic c glide planes (yz plane) are termed as Y walls. Similarly, the domain walls, parallel to the crystallographic xz planes and perpendicular to c glide plane are termed X walls. As we shall see, this complex crystallography of the single crystal has a profound effect on the structure and properties of the domain walls.
63 Figure 3 4 Schematic of (a) Y1, Y2 and Y3 axis along with three oxygen ions, O1, O2 and O3. (b) Crys tallographic direction of Y1, Y2 and Y3 axes parallel to the Y wall in the y z plane. (c) Crystallographic direction of Y1, Y2 and Y3 axes parallel to X wall i.e., parallel to an x z plane and at a 30 z plane. The top oxygen layer of (b) and (c) is a duplicate of the bottom oxygen layer. 3. 1. 3 Temperature E ffects on LiNbO3 For empirical calculations, t he ferroelectric polarization of LiNbO3 can be understood in terms of the electronic polarizations associated with the coreshell interactions, the displacements of Li and Nb ions from the high symmetry positions of the paraelectric phase, and the associated small oxygen displacements. Figure 35(a) shows the temperature dependence of the polarization as determined with the Jackson potential. The polarization decreases quite gradually up to 1400K, above which it rapidly drops to zero at 15001510K, which we take as an estimate of the Curie temperature for this potential. Th is value is in good agreement with the experimental Curie temperature (a) (b) (c)
64 of 1483K.  This result is also supported by an analysis of the distribution profile of the displacements of the Nb and Li atoms from the centers of the oxygen octahedra and from the oxygen planes. Well below the Curie temperature, the main effect of the incr eased temperature is a broadening of the distribution of cation positions about their displaced values with only a small decrease in polarization. At temperatures above 1400K, there is a significant shift of the cations to positions that are symmetric relative to the oxygen sublattice, with an associated significantly decreased polarization. Figure 35. (Color online) (a) Temperature dependence of polarization of LiNbO3 from experiment and MD calculation; (b) Temperature dependence of averaged displac ement of Li and Nb ions from an oxygen plane and the center of the O for various temperatures. (a) (c) (d) (b)
65 Figure 35(b) shows the average displacements of the Li and Nb ions from their cryst allographically expected sites, both of which decrease to zero at the Curie temperature. We see that the Nb displays displacive dynamics, figure 3 5(c), while the Li ions show order disorder behavior, figure 3 5(d), consistent with the results previously o btained with the Tomlinson potential  These results differ from the previous results using the Tomlinson potential  which showed that the Nb displacements went to zero well below the Curie temperature; this effect is not seen for this potential. 3 1. 4 Domain W alls in LiNbO3 Ferroelectric materials belong to a class of crystals whose low symmetry supports a spontaneous polarization along one or more crystal axes in the absence of an electric field  This spontaneous polarization consequently appears as an apparent surface charge density with depolarizing field (ED). The ferroelectric phase will try to minimize its energy associated with t he polarization by twinning. This causes the crystal to divide into many oppositely polarized regions called domains. The type of domains formed in the crystal depends on the crystal structure. For LiNbO3, the longitudinal polarizations, come from two deg ree of freedom (up and down), make only 180. Figure 3 6 shows a schematic representation of the formation of 180 domains. Alternate areas with opposite charges are generated on the surface to minimize ED, resulting in an overall zero polarization in the crystal. The domain walls will affect the structure and properties of ferroelectric materials.
66 Fig ure 3 6. Schematic of (a) surface charge associated with spontaneous polarization, (b) formation of 180 domains to minimize electrostatic energy  Figure 3 7 (a) Schematic view parallel to the (0001) plane showing two structurally distinct domain wall structures; (b) planar structure of Y wall on )0 2 11 ( p lane; (c) planar structure of X wall on ) 0 1 10 ( plane. The domain walls in LiNbO3 separate regions with polarizations in the +z  and z [000 1 ] directions. Scrymgeour et al.  determined that two different crystallographic planes can define the domain wall. In figure 3 7(a) the X walls and Y wall are projected onto the (0001) plane. The Y wall lies parallel to c glide planes, thus containing only cations or anions, but not both. The X wall lies perpendicular to c glide planes, and thus contains both cation and anions. Figures 37(b) and 37(c) show slices (a) (b) (a) (b) (c)
67 through these domain walls on one side of the wall, parallel to (01 1 0) for the Y wall and parallel to (10 1 0) for the X wall. The t wo structurally distinct domain wall planes lead to two variants of domain walls. Every Y wall lies parallel to one of the Y axes; every X wall lies parallel to one of the X axes, symmetrically between two of the Y axes. For the Y wall parallel to the Y1 axis, see figure 3 4(b), we can expect the O1 oxygen ions to be displaced parallel to the wall (i.e., along the Y1 axis), a nd the O2 and O3 ions to be displaced along the Y2 and Y3 axes at 120 to the domain wall. Because, the X walls lie along planes at 30 walls, they lie perpendicular to one of the Y axes ( figure 3 4(c)). For the X wall plane per pendicular to Y2 axis, as shown in figure 3 4(c), the O2 ions mainly displace parallel to the Y2 axis, i.e., perpendicular to the domain wall. By contrast the O1 and O3 displace mainly at 30 from the X wall. As we will see in C hapter 4, because of their different crystallographies, the two domain walls show different atomic structures and different patterns of polarization in the domain wall region, and thus have different domain wall energies. 3 2 Conductive P olymer s 3. 2 .1 Introduction A ll carbon based polymers were generally accepted as insulators for a long time. Plastics, which are a particular type of polymers, w ere mainly utilized as insulating material s, such as inactive packaging After electronic conductivity i n polymers were discovered about 30 years ago  polymers have been incorporated into displays and batteries. However, low stability poor processibility and high manu facturing costs limits their applications to small scale applications. They are extensively studied for their
68 potential application to additional electronics devices such as o rganic solar cells, printing electronic circuits, organic light emitting diode, actuator, electrochromism, supercapacitors, biosensors, flexible transparent displays, electromagnetic shielding and transparent conductor s[98 101] The new class of polymer known as intrinsically conductive polymer or electroactive polymers changes the narrow perspective on polymer s. As a new technology, the potential uses of conductive polymers for industry are quite significant. There are many known conductive polymers including polyacetylene, p olyphenylene, pol yphenylene sulfide ( PPS ) p phenylenevinylene ( PPV) p olypyrrole, p olythiophene, and p olyaniline. Table 34  compar es the properties of different conductive polymers. As we can see from the T able 34, polythiophene received wide attention for indus trial application s due to their electrical properties, low cost, light weight, and high processability. These properties have led polythiophene to be used for electronic and optical devices including light emitting diodes, electrochromic devices, field eff ect transistors, organic photovoltaics, sensor films, recording materials and rechargeable batteries [103, 104] However, conventional growth methods lead to thermal and mechanical instabilities in thiophene thin film s. [97, 105] S urface polymerization by the ion assisted deposition (SPIAD) method has been proposed to increase the stability of polymer films by controlling their chemistry and morphology on the nanometer scale In this thesis, polythiophene is examined to gain a n understanding of polymerization and determine the crosslinking mechanism s that dominate during the SPIAD process, which is not commonly used in industry at this time.
69 Table 34. The co mparison on conductivity, stability and Processing attributes of representative conducting polymers  Polymer Conductivity 1 cm 1 ) Stability (doped state) Processing possibilities Polyacetylene 10 3 ~ 10 5 P oor L imited Polyphenylene 1000 P oor L i mited PPS 100 P oor E xcellent PPV 1000 P oor L imited Polypyrrole 100 G ood G ood Polythiophene 100 G ood E xcellent Polyaniline 10 G ood good 3. 2 .2 Structure of P olythiophene Figure 38. The comparison of calculated thiophene structure using (a) Gaussia n with B3LYP/6 31G basis set and (b) SIESTA with GGA/PBE and double zeta polarized basis set Polythiophene is the polymerized form of thiophene, which has the formula unit, C4H4S. Figure 38 compares the thiophene structure calculated using two different calculations. Figure 38(a) shows the calculated structure with a hybrid method based on quantum chemical and DFT methods B3LYP/6 31G using Gaussian software package. The thiophene structure in figure 3 8(b) is obtained using DFT within the SIESTA package with GGA/PBE pseudopotential s and a double zeta polarized basis set. Normally, quantum chemical method is believed more accurate than normal DFT (a) (b)
70 method. The two different calculation methods predict very similar structures, with acceptable differences in bond lengths and angle s. Figure 39 The structure of terthiophene; three thiophene rings are arranged in an antiparallel to neighbor The conjugated chain of polythiophene plays an important role in the properties of conducting polymer and nonlinear optics.   Because the polythiophene displays fluorescence and Raman spectral features similar to terthiophene films  terthiophene molecules are examined here The structure of terthiophene was previously investigated by X ra y and a quantum chemistry study  The structure of polythiophene is determined to show trans position with a continuous antiparallel ordering of thiophene rings, see figure 3 9(a)  Figure 39 also shows a comparison of the middle terthiophene ring from (b) an X ray study and (c) hybrid quantum mechanical calculation with B3LYP/631G basis set (c). The two different methods predict similar bond length and angles as the single thiophene structure. The X ray examination determined the density of terthiophene crystalline as 1.503 g/cm3. (a) (b) (c)
71 The film densities for our simulation study are chosen to be similar to the experimental result. 3. 2 3 Conduction M echanism in P olythiophene One early explanation of the conducting mechanism in polymers was a half filled valence band created by continuo bonds along the polym er chain  However, polymers favor bond alteration, creating alternative short and long bonds. The energy gap created by bond alternating is a bout 1.5 eV. This band gap can be controlled by doping with either donors or acceptors. Thus it is necessary to examine the way that charge is stored along the polymer chain and can be transferred from pol ymer chain to polymer chain. Charges can be stored in polymer chain s during oxidation or reduction   Different types of oxidizing and reducti ng dopants are available that do not decrease without l owering its stability. In an oxidation process positive charge can be localized in a small section of the chain. The localized charge distorts the geometry of the polymer chain, but decreases the ionization energy and increases the electron affinity of the polymer chain. Eventually, the polymer chain becomes electronic ally conducting by charge delocalization  Figure 310 shows the changes in polythiophene chain during oxidation. Initially, polythiophene has a consecutive backbone chain of heterocyclic thiophene, see figure 3 8(a). With oxidizing dopants such as i odine, arsenic pentachlor ide, iron(III) chloride and nitrosyl ( NOPF6 ) an electron is removed from a porbital of the backbone and free radical and spinless positive charge in the polymer chain are created see figure 3 8(b ). T he radical and positive charges of the polymer chain a re coupled to each other since the separation of these defects requires a considerable amount of energy. T he
72 combination of radical and site charge is called as a polaron. Upon further oxidation, the free radical of the polaron is removed and creates two spinless posi tively charged defect state, the so called bipolaron see figure 3 10(c). B ecause the system can lower the energy by creating a bipolaron, two polarons are replaced with a bipolaron at higher doping level s. Figure 310. The variation of polythiophene polymer chain with ptype doping.  Figure 311 shows the variation of the electron bands during doping. The polythiophene will be in neutral state without doping of oxidative dopants. With a single dopant, the polaron state will be created. The polaron will create new electron states between conduction band and valence band, see figure 3 11(b). A single unpaired electron will occupy the lower energy state of this defect level. A dditional oxidation will combine two polarons into a bipola ron. As we can see from figure 3 11(c), in the bipolaron state two spinless positively charged energy level s, are generate d by two acceptors. With continuous doping, the multiple bipolaron energy state will overlapped and a bipolaron band will be created, see figure 3 11 (d). (a) (b) (c)
73 Fig ure 3 11. The variation of electron band during oxidation.  The nature of charge transfer between polymer chains is not completely understood. However, the main mechanism for charge transport is accepted as a combinatio n thermally activated hopping and tunneling between charged particles such as the polaron and bipolaron. 3. 2 4 Growth of P olythiophene The p roperties of conducting polymer thin films can vary depending on the growth process. The optimized process can impr ove the device performance by controlling the charge injection, mobility and recombination efficiency Thus, development of a processing method which can control film chemistry and morphology at the nanoscale, is a challenge for the application of conducti ng polymer s [112, 113] Direct thermal deposition and solution based growth methods are typical growth techniques for conductive polymer thin films For d irect thermal deposition, lower molecular weight oligomers are deposited on the substrate surface. During deposition, the small molec ules of gaseous species sublime back to the chamber. For s olutionbased growth methods printing, casting or spin coating techniques are used for the growth of polymer. This method cause s inadequate or uncontrolled ordering of polymers since the self (a) (b) (c) (d)
74 assem bly of molecules is the only way to decide the morphology of polymer. Therefore, both growth methods lead to difficult i es in control ling the film thickness and morphology. Surface polymerization by ionassisted deposition (SPIAD) [105, 114 117] which deposit s hyper thermal polyatomic ions and thermal neutrals in vacuum simultaneously, has been developed to avoid thermal and mechanical instability from conventional growth technology by controlling surface morpholog y B oth mass selected and nonmassselected ions can be used for SPIAD. Selective chemical bonds with self assembled monolayers are achieved by m ass selected organic cations [118, 119] N onmassselected ion deposit ion is suitable for prototype manufacturing processes for new conducting polymers [120, 121] The SPIAD method al lows more control over film properties includ ing ion energy, ion structure, ion kinetic energy, charge state, the ion/neutral ratio, and substrate temperature  SPIAD method can not be only used for optimizing morphology, film thic kness, but also fine tuning of film optical band gaps and other optoelectronic properties by modifying the chemical and morphological structure of the film [105, 115, 123] For example, SPIAD films of polythiophene and polyphenyl display narrow er band gaps with reduced barriers for hole injection compared with their evaporated film counterparts [114, 117] The poly thiophene film s produced by SPIAD also display fluorescence in the UV and visual region. This property makes polythiophene a candidate for organic photovoltaic device applications. P revious SPIAD studies reporte d mechanisms for the ioninduced surface polymerization reactions are pronounced at the specific ion/neutral ratios and ion energies [105, 115] However, SPIAD failed to form a single oligomer state by generati ng a distribution of species [114, 115] Also, different experiment al result s are achieved by using atomic vs.
75 polyatomic ions Polyatomic ions behaved both as catalyst s and reagent s by energetically inducing polymerization and forming adducts with the neut ral reagent, respectively  Thus, t he util ization of SPIAD for more sophisticated applications require s comprehensive understanding of various mechanisms cont ribut ing to the overall film formation M anag ing the polymerization process during growth of polythiophene make s it possible to control the optical properties of p olythiophene for applications in nonlinear optics and photovoltaic devices.
76 CHAPTER 4 DOMAIN WALLS IN FERROELECTRIC LITHIUM NIOBATE 4.1 Introduction Ferroelectric materials have potential use in diverse applications such as nonvol atile Ferroelectric Random Access Memory (FRAM)  electro optic modulators  and frequency converters  Ferroelectric LiNbO3 has an R3c structure, has excellent piezoand pyroelectricity propeties, is photo refractive, and displays nonlinear optical propert ies [69, 127] Thus lithium niobate (LiNbO3) has many applications in optoelectronics  nonlinear optics  and micro electromechanical devices  As we discussed in Sec. 3.1.1, the hig h temperature paraelectric trigonal ( R 3 c) structure  becomes ferroel ectric at the Curie temperature by los ing its inversion symmetry along the c axis, resulting in space group R3c Since all of the ion displacements associated with the ferroelectricity are parallel to the c axis (i.e., along ), the paraelectric and ferroelectric phases look i dentic al in the top view, figure 4 1 Gopalan et al.  identified two crystallographically different 180 domain walls between domains of opposite polarization. As discussed in Sec. 3.1.2 an X wall lies parallel to a m ixed anioncation plane, while a Y wall lies parallel to alternating planes of cations only and anions only. Gopalan et al derived a relationship between the domain wall width and energy based on Ginzburg LandauDevonshire theory  However, the detailed atomic structure and energetics of these domain walls have not been determined. The objective of this chapter is therefore, to characterize the structure and energetics of these domain walls. The approach we will take is atomic level simulation
77 with empirical potentials, validated against electronic structure calculations at the level of density functional theory. Figure 41. Structure of LiNbO3 (a) Edge view of polarization down (left), paraelectric (center) and polarization up (right) states (b) Planar view of (0001) plane. Three equivalent hexagonal axes are drawn X1, X2, X3 and Y1, Y2, Y3, with X1 and Y1 paral lel to the Cartesian x and y axes. 4. 2 Domain W all P reparation To study these domain walls, we initially prepare two single crystals of the same size: one with polarization up, the other with polarization down. We join these two structures and apply periodic boundary conditions, thereby forming two domains, separated from each other by two crystallographically identical (X or Y) 180 walls. We determine the zerotemperature equilibrium structure of each system by allowing all atoms to move to zerofo rce positions, thereby minimizing the total energy of the system. The system size for DFT calculation is eight unit cells with 240 atoms and 1440 electrons. For empirical study, we increase the system size into 36x4x3 unit cells
78 including 12960 atoms in or der to prevent interactions between domain walls. The simulations are performed under constant pressure conditions for both empirical and DFT calculations in order to capture the change in a volume as well as atomistic displacements. Both methods show less than 1% change in lattice parameters and volume from bulk ferroelectric values. 4.3 Domain W all E nergy Although the X walls and Y walls are defined to be parallel to the (10 1_0) and (11 2_0) planes respectively, the actual position of the walls is not determined crystallographically but from minimization of the total energy of the system. The solid and dotted lines in figure 4 2 ( a) denote planes of high symmetry upon which the domain walls might be expected to lie. For the X walls, the two high symmetry positions are along the ion planes and between the ion planes. For the Y walls, the three highsymmetry positions are at the cati on plane, at the anion plane, and between two anion planes. Both the empirical (E) and DFT (DFT) calculations show that the equilibrium positions of the Y wall lies halfway between the anion planes. The calculations yield a Y wall energy of EY DFT=160 mJ/m2; and EY E=230 mJ/m2. That is, the domainwall energy is defined by the difference in energies of a system with domain walls and a single crystal with the same number of ions, divided by the total area of domain wall in the system. The two methods also agr ee that the equilibrium position of the X wall lies equidistant between two ion planes, rather than at an ion plane. The X wall displays two different variants (XI and XII) with very similar energies (EXI E=260 and EXII E=255 mJ/m2), but different polarizati on profiles; the structural differences between these two variants,
79 and their origin will be discussed below. The DFT calculations show the same two different structural variants. However, it is not possible to extract individual domain walls energies from DFT, but only the averaged energy value: (EXI DFT +EXII DFT) = 200 mJ/m2. Although, as is to be expected, the DFT and empirical values do not precisely agree, both methods do predict the Y walls to have a lower energy than the X walls. This conclusion is a lso consistent with the conclusions of Scrymgeour et al.  based on calculations using Ginzburg LandauDevonshire theory. It is also possible to estimate the energy barriers for the motion of the domain walls by determining the energy of likely transition states. In particular, we can force the domain wall to be at a specific plane by setting the positions of the ions in the plane to be the average of the two opposite polarization states. For the Y wall, it seems logical that the cation plane should be the transition state because it is symmetric, lying half way between two equivalent equilibrium positions. For this highly symmetric state, the empirical simulations yield an energy for the transition state of EE* Y=485 mJ/m2, see figure 4 2 Assuming that this state is actually the transition state, the empirical calculations thus yield an energy barrier of 255 mJ/m2, w hich is almost identical to the Y wall energy itself. However, The DFT calculations yield an energy for the transition state of EDFT* Y=310 mJ/m2, corresponding to an energy barrier of 150 mJ/m2, which is only a half of the value obtained from the empirical calculations. The transition states are slightly different for the two X wall variants. The energies of the two X walls when centered at ion planes are calculated to be EE* XI=375 mJ/m2 and EE* XII=394 mJ/m2. Thus the estimated energy barrier for X wall mo tion is 115 mJ/m2 for the XIwall and 140 mJ/m2 for the XII wall. The DFT calculations yield an average
80 transitionstate energy of (EDFT* XI +EDFT* XII)=230 mJ/m2, corresponding to an average energy barrier of 30 mJ/m2. Again, there is a difference in the precise values from the empirical and DFT calculations; however, they both agree that the Y wall has a higher energy barrier for domain wall motion than the X walls. Figure 4 2 Energies of domain walls in equilibrium and in transition states. The equilibrium position is determined as equidistance between two anion planes for Y wall and center between the cationanion mixed planes for X walls. The Y wall has a lower energy, but higher energy barrier for the wall motion than X wall. The lines are merely guides to the eye. 4. 4 Atomic D isplacement around D omain W alls The domain walls break the crystallographic symmetry of the single crystal. Because of their proximity to two regions of opposite polarization, the atomic structure in the planes close to the dom ain walls should be a close approximation to the paraelectric structure, converging to the ferroelectric structure at some distance from the domain wall. Moreover, the domain walls produce displacements, not only parallel to the polarization direction (z), but also in the x y plane, with both components normal to and
81 parallel to the domain wall. Following the notation of Scrymgeour et al.  we denote displacements and polarizations normal to the domain wall as normal or n, and those parallel to the domain wall (but orthogonal to the bulk unixial polarization) as transverse or t. Close to the domain walls, all three components of polarization can be manifested. 4. 4 .1 Y walls, P arallel to cglide P lanes The Y wall lies between two anion planes separated by only 0.64 The oxygen ion displacements in these planes are therefore large, as shown in figure 4 3 F or the Y1 domain wall, the O2 ions are closest to the domainwall plane (~0.30) and show the largest Z displacement from their ferroelectric positions, 0.23. Because this displacement is almost as large as the displacement from the paraelectric to ferroe lectric positions in bulk, we can understand this as the O2 ions moving back to positions close to their positions in the paraelectric phase. The next closest, O3 ion, is ~0.35 from the Y wall and shows a slightly smaller Z displacement of 0.21. The O1 ion which sits further away yet, 0.65, shows an even smaller Z displacement. The oxygen ion positions finally recover their bulk ferroelectric positions at a distance of ~5.8 away, as defined by the displacement from the ferroelectric positions dropping below 10% of its maximum value. Since at equilibrium the Y wall lies between two anion planes, the closest cation planes are much further away, 1.3 and thus undergo significantly smaller Z displacements than the anions. Interestingly, the symmetries of the cations are also broken by the Y wall. However, in this case the symmetry breaking does not arise from the wall itself, but as a secondary effect caused by the breaking of the symmetry of the O1, O2 and O3 ions. F or example, the top Li in figure 3 4 (c) is in an oxygen cage of
82 different orientation than the oxygen cage enclosing the lower Li ion. Therefore, although they are the same distance from the Y wall, these two distinct LiI and LiII ions undergo different Z displacements (0.02 and 0.00). Similarly, the distinct NbI and NbII ions also undergo Z different displacements (0.019 and 0.016). The maximum cation displacements of both Nb and Li ion are, unsurprisingly, in the planes closest to the Y wall. These cation displacements are in the direction opposite to the displacements associated with the ferroelectricity, indicating that the cations are displaced back towards their paraelectric positions. The cations recover their ferroelectric positions ~11.7 from the wall, which is twice the distance of that required by the anions. Figure 43 Displacement pattern in z direction of (a) oxygen ions and (b) cations around Y wall. Note the difference in scales of displacements. Because the oxygen ions near the Y wall shift oppositely poled ferroelectric states to approximately paraelectric symmetric position at the domain walls, the oxygen ions on opposite sides of the domain wall displace in opposite transverse directions. If an oxygen ion moves along the positive y axis along on one side of the wall, t he crystallographically equivalent oxygen ion on the opposite side moves along the negative y axis. For displacements normal to the wall, symmetric oxygen ions on opposite sides displace equally either toward or away from the Y wall regardless o f (a) (b)
83 ferroelec tric state. Figures 4 4 (a) and (b) shows the transverse and normal displacements of the oxygen ions in the region around a Y wall; Figs. 4 4 (c) and (d) show the corresponding displacements for the cations. Here, the positive and negative signs of the trans verse displacement correspond to movement in the +y and y directions. F igure 4 4 (c) shows the transverse displacements of cations at the nearest cation plane of the uppolarized region; the transverse displacements for down polarization side are equal and opposite. For normal displacements, the positive and negative signs denote displacement movement away from and towards the Y wall respectively. The transverse displacements of the three different nonsymmetric oxygen ions can be analyzed in a similar manner. As figure 4 4 (a) shows the O1, O2 and O3 ions on the one side of the wall displace in opposite directions from their counterparts on the other. Because the closest anion to the wall, O2, shows the largest displacement in z direction, it also has the largest transverse displacement, O2T~0.066. As we can see from figure 4 4 (b), the normal displacements of oxygen ions show the same pattern on the two sides of the domain wall. For the Y1 axis, the O1 ions have the same x positions, but different y posit ions in the two different ferroelectric states; they thus displace along only in the transverse direction. Specifically, the simulations show that the Y wall pushes the O1 ion away from the wall; the normal displacement is O1N~0.02 at the nearest oxygen plane. The normal displacements of the O2 and O3 ions are almost equal and opposite to each other (O2N ~ 0.15 and O3N ~ 0.13) because the Y2 and Y3 axes are oriented at 60 wall. These displacements ar e larger than the ferroelectric
84 displacements, 0.11, resulting in additional normal displacements of oxygen ions which break the paraelectric symmetries of oxygen ions near the Y wall. Figure 44 Atomic displacements of O1, O2 and O3 ions in (a) transverse, and (b) normal direction directions near a Y wall. Schematic view of corresponding cation displacements along (c) transverse and (d) normal directions. Relative larger + signs (red) and smaller + signs (black) represent larger and smaller cation displacements. For transverse displacements, the plus and minus sign indicate the positive and negative displacement along Y1 axis. For normal displacement, the plus sign indic ates that the cation moves away from the Y wall, while the minus sign designates that cation displaces towards Y wall. The broken symmetry of the oxygen ions also causes the cations to break the uniaxial symmetry around the Y wall. Figures 4 4 (c) and (d) show schematic views of (b) (b) (a) (c) (b) (d)
85 the transverse and normal cation displacements for the cation plane closest to the Y wall, at which the displacements are largest. The LiI and NbI show relatively larger transverse displacements than the LiII and NbII ions. As we ca n see from figure 4 4 (c), the Li and Nb ions displace in opposite transverse direction, although LiI and LiII ions displace in same transverse directions, as do the two types of Nb ions. Figure 4 4 (d) shows that all of the cations near Y wall are displaced in the normal direction away from Y wall. Because of their weaker electrostatic interactions with the oxygen ions, the Li ions displace more than the Nb ions. Because of the differences in the symmetries of the displacement pattern along the transverse and normal direction, the normal displacements of cations converge to zero at relatively longer distances from the Y wall (~11. 67 ) than the transverse displacement (~6.48 ). 4. 4 .2 X walls, P erpendicular to cglide P lane The structure of the X wall is ve ry different from that of the Y wall for two reasons. Most importantly, the X wall lies parallel to planes that contain both anions and cations. More subtly, because it lies at 30 wall, the X wall is normal to one of the Y directions, and at 30 on the oxygen sublattices, this creates t wo X wall variants, depending on whether the other oxygen displacements point into the X wall or point out of the X wall. In particular, as figure 4 5 (a) shows, when the Y2 direction is perpendicular to the X walls, the O2 ions mainly displace in the normal direction into the domain wall, leading to the XIwall. An O2 ion sitting near the XIwall moves in the +Y2 direction for the up polarization and in the Y2 direction for the down polarization state. Both of these displacements represent normal displacem ents toward the XIwall. Likewise, figure 4 5 (b) shows that when the Y3 direction is perpendicular to the X wall, the O3 ions displace along the
86 normal direction away from the domain wall. This direction is opposite to that for XIwall, thus creating a dif ferent variant, the XIIwall. The positive displacement of the O3 ion along Y3 axis near the XIIwall for up polarization and negative displacement for down polarization indicates the movement of O3 ions away from the XIIwall. Figure 45 Crystallogra phies of the XIwall and XIIwall. Figures 4 6 (a) and (b) show the displacement patterns of the anions and cations in the Z direction around the XIwall. The Z displacement is largest in the planes closest to the domain wall, decaying to zero at ~5 from the wall. The corresponding cation displacements are much smaller, with the cations moving closer to their paraelectric positions. Figures 4 6 (c) and (d) show the anion and cation displacements around the XIIwall. The Z displacements of the oxygen ions ar e quite similar for the two X wall variants, except that the maximum Z displacement of the XIIwall is slightly larger (0.19 vs. 0.17) However, comparing Figs. 4 6 (b) and (d), there are significant differences in the cation displacements. For the XIwal l, the maximum Z displacement of both cations (a) (b)
87 is observed at the nearest plane. For the XIIwall, the Li and Nb ions move by the same amount in the Z direction at the closest plane, thereby compensating for the displacement of oxygen positions in the Z dir ection. As a result, the Li ion closest to its paraelectric position is at the second nearest plane. In contrast to the Y walls, there is only one type of Z displacement pattern for Li ion and Nb ion in the X walls Figure 46 The Z displacements patterns around XIwall and XIIwall. (a) The anion displacement pattern around XIwall (b) The cation displacement pattern around XIwall (c) The anion displacement pattern around XIIwall (d) The cation displacement pattern around XIIwall The two X wall var iants also create two distinct patterns of normal and transvers e ion displacements, see figure 4 7 (a). The negative sign in the normal displacement for the O2 ions indicates that they move toward the XIwall; the O1 and O3 ions move away from XIwall by eq ual amounts along their internal Y ax es. In a similar manner, figure 4 7 (b) shows the displacements of the anions near XIIwall. The maximum normal (b) (a) (c) (d)
88 displacement (~0.17 ) at the plane closest to the wall disappears ~5 away from both of X walls. Figure 4 7 Normal displacements of three different oxygen ions near (a) XIwall and (b) XIIwall. (c) transverse displacements of three different oxygen ions around XIwall. (d) normal and transverse displacement of cations around XIwall. Normal displacement o f cations around XIIwall at (e) the nearest and (f) the second nearest plane. (a) (d) (b) (e) (c) (f)
89 Figure 4 7 (c) shows the transverse displacement of the anions around the XIwall; the displacements for the XIIwall are very similar. All of the transverse displacments of the anions are cancelled by equal and opposite displacement patterns. The displacement of cations around the two X walls are shown in figure 4 7 (d) through (f). Unlike the Y wall, only one normal displacement pattern for each cation is observed for both X wall s. Figure 4 7 (d) shows the cation displacement around the XIwall. The first +/ sign for each cation represents the normal displacement, while the second +/ sign designates the transverse displacement of the cations. The Li and Nb ions, which move in opposite normal direction, show the largest displacement at the second nearest plane from XIwall. Both Li and Nb ions in the first plane move away from the XIIwall, figure 4 7 (e). However, the Li ions in the second plane move into the XIIwall while the Nb ions shift away, figure 4 7 (f). The transverse displacement of the cations are very similar to the anion displacement. Unlike the normal and Z displacement pattern of cations, the nonsymmetric cations, LiI, LiII and NbI, NbII, show opposite displacements i n the transverse direction. Therefore, the equal and opposite transverse displacement of each cation along Z direction are cancelled out. 4. 5 Polarization at Domain Walls The switching of the polarization at a ferroelectric domain wall is generally considered to be relatively abrupt, taking place over only one or two unit cells  However, the relatively long distances (up to ~10) over which the atomic displacements return to bulk values suggest that the polarization may be governed by a similar length scale. In this section, we characterize the transition in the uniaxial polarization, and in the components of polarization normal to and transverse to the domain wall. As we shall see, these more complex polarization patterns manifest some
90 features reminiscent of those present in ferromagnetic domain walls in which the spin state switches from one to another over, typically, many tens of unit cells.  However, because of the much greater strength of electrostrictive interactions compared to magnetostrictive interactions, the nonuniaxial components are relatively small and characterized by shorter length scales than their magnetic counterparts. 4. 5 .1 Y wall s We have determined the polarization along the direction to the normal Y wall using both empirical potentials and DFT. As is customary,  for the DFT calculation only the ionic contributions to the polarization are considered. For the calculation, we use the Born effective charge developed by Veithen and Ghosez  for LiNbO3 ferroelectric system. As figure 4 8 (a) shows, the two methods give very similar profiles in the uniaxial contribution to the polarization. The asymmetric pattern of up and down polarization arises from the shift in the Y wall from t he anion plane to the equilibrium positions; as a result, one of the domains is wider than the other. We note that using the full ionic charge in analyzing the DFT results, rather than the Born effective charges, leads to almost complete overlap with the polarization profiles determined from the empirical potentials. We take the good agreement between the DFT and empirical results as validation for the larger scale calculations using the empirical potentials, which show a change in Pz over a very short dist ance, see figure 4 8 (b). The non uniaxial displacements of the ions generate nonuniaxial contributions to the polarization. Figure 4 8 shows (c) the normal and (d) the transverse polarization around the Y wall. Both are anti symmetric around the domain wall, with zero polarization contributions at the center of domain wall. The normal and transverse
91 2 2 at 2.58 decaying to zero ~9 away from the wall. Figure 48 (a) Comparison of z direct ional polarization (Pz) around Y walls, determined by DFT and empirical simulation methods. (b) Uniaxial polarization profile of larger simulation cell with Y wall. Variation of inplane polarization near Y wall; Normal(c) and transverse. (d) polarization profile at each side of Y wall show exactly symmetric with opposite sign. A 180 ferroelectric wall is typically understood as Ising type see figure 4 9(a) ; by contrast a 180 magnetic wall is typically understood as Bloch or Nel type.  The presence of signific ant polarization components Pn and Pt of up to ~2.66% of Pz gives rise to a Blochlike and a Nel like contribution to the predominantly Ising like Y wall.  The polarization vector, Pz+Pt leads to a Bloch like r otation by angle B in the t z plane of the domain wall figure 4 9 (b), whereas polarization vector, Pz+Pn leads to a Nel like rotation by angle N in the n z plane perpendicular to the domain wall figure 4 9 (c). The maximum rotation angles of B~1.52 a nd N~0.55 occur at a distance
92 xn~ 1.3 from the Y wall where the Pn and Pt components are maximum. Thus the Y wall has a relatively strong Blochlike rotation behavior and a weak Nel like rotation behavior. Figure 49 Schematic representation of ( a) three different polarization, Pz, Pt and Pn; (b) Bloch like wall is created by vector sum of Pz and Pt and (c) Nel like wall is created by vector sum of Pz and Pn. B and N is calculated by maximum rotation angle of Pz at the position where Pz returns to bulk polarization value. 4. 5 .2 X walls As in the case of the Y wall, the empirical and DFT calculations yield very similar polarization patterns around the X walls, see figure 4 1 0 (a). To understand the polarization pattern around the X walls without the limitations of small system sizes, we again performed larger scale simulations with the empirical potentials. For the geometry (a) (b) (c)
93 used in thus study, the XIwall lies at center of the system, while the XIIwall lies on the periodic boundary at 70. As w e can see from the figure 4 10 both X walls show very sharp transition regions between polarization states. The bulk polarization of 2 is reached at ~6 away from XIwall and at ~9 away from XII wall. Figure 41 0 (a) Comparison on Z directi onal polarization around X walls between DFT and empirical simulation methods. (b) Uniaxial polarization profile of larger simulation cell with X walls. Two different X walls are created at the center(XIwall) and at the boundary(XIIwall). (c) The variati on of inplane polarization around X walls; Normal polarization profile near XIwall(c) and XIIwall(d). Because the atomic displacements in the transverse direction are cancelled by equal and opposite movement of each ion, figure 4 7 (c), there are no transverse polarizations at the X walls. Figures 4 10(c) and (d) show the normal polarization profile around the XIwall and XIIwall respectively. The different normal polarization patterns near the two X walls are the result of different displacements of th e cations. With the (a) (b) (c) (d)
94 largest normal displacements at the second nearest plane, the maximum polarization of 2 out of XIwall is observed at 5.2 away, slowly decreases in a monotonic manner, reaching 10% of its maximum value about ~15 away. The c hange in normal displacement pattern of cations between the first and second nearest plane to XIIwall 2 at 1.49 and leads to oscillatory convergence over a distance of ~9. Considering the maximum B loch like and Nel like rotation angles near the X walls, the angles are B~0 and N~0.56 for the XIwall and B~0 and N~0.65 for the XIIwall. Thus, both of X walls show only Nel like rotation. 4.6 Domain W all W idth The Ginzburg LandauDevonshire ( GLD) [136 138] thermodynamic theory successfully explains the relationship between the variation of the polarization and domain wall width in the transition region for both perovskite materials  and for trigonal ferroelectrics such as LiNbO3 and LiTaO3 In the transition region, the total free energy of a ferroelectric material can be understood as arising from four different energy components [139, 140] : a free energy associated with the second order phase transition, a free energy associated with elastic lattice displacement, a free energy derived from structural distortion by nonuniform polarization, and a free energy associated with gradient in the polarization. Solvin g the free energy equation for energy minima, Scrymgeour et al  show ed the relationship between three components of polarizations, Pn, Pt, and Pz as: 0 0tanh ) (n n zx P x P ( 4 1 ) and
95 5 3 3 2 1 3 2 10 0 1 ) ( ) ( ) (z z z t t t n n n n z n t n nP P P x P x P x P ( 4 2 ) where, ik ( i=n,t; k=1,2,3 ) are the coefficients of the polynomials. Our atomic scale simulation study also shows that non uniaxial polarizations are strongly correlated to the uniaxial polarization. While Pz shows a kink type polarization variation, reaching zero at the center of the wall, Pn and Pt are zero both at the center of the wall and far away from the wall. The si mulations discussed above show that inplane polarizations reach their maximum value at approximately distance from the domain wall at which Pz returns to bulk value. Thus, Pn and Pt returns to zero bulk value at relatively longer distance than Pz. Because no complete model has been developed for the nonuniaxial polarizations of trigonal systems, we only show the domain wall width along with Pz component in this study. A fit of Pz to a single hyperbolic function, equation ( 4 2 ), yields o ~2.160.05. Figure 4 1 1 show s normalized polarization profiles of Pz along xn around the Y wall. The solid line is a fit to simulation results for Pz. Figure 41 1 A fit of Pz to hyper tangent function around Y wall. Eq uatio n ( 4 2 ) yields the domain wall width o ~2.160.05.
9 6 The two X walls yield essentially identical domain wall widths: o ~2.100.03 for the XIwall and o ~2.120.03 for XIIwall, as shown on figure 4 1 2 The fits of Pz to the hypertangent function lead very similar domain wall width for three different types of domain walls. The absence of analytic results for nonuniaxial contributions to the polarization indicates that there is still scope for fur ther developments in the GLD analysis of these domain walls. Figure 41 2 Fits of Pz of (a) XIwall and (b) XIIwall to Equation ( 4 2 ). The fits of Pz yield almost same pattern between Y wall and two X walls. Further analysis on inplane polarization i s required for the comparison. 4. 7 Threshold F ield for D omain W all M otion Ishibashi  an d Suzuki  showed that the width of the domain wall controls domain wall motion with the activation energy for wall motion being inversely (a) (b)
97 proportional to its width. In their analysis, the activation energy for domain wall motion, can be expressed as a free e nergy difference between equilibrium state and transition state, activationF : 2 / 7 4 2 3 0 2 12 exp e a w a w w P F F Fo o s m Equilibriu Transition activation ( 4 3 ) where 1 and 2 are the constant in the GLD polynomial, 0w is domain wall half w idth, SP is spontaneous polarization and a is lattice constant. From this the threshold field for domain wall motion can be determined as the activation energy dividing by SP w0  In the S ec 4.6 we obtained domain wall widths by fitting the Pz profiles. The lattice parameter for Y wall shows a ~5.1684 and spontaneous polarization 2/ 4 62 ~ cm C PS The coefficients for the GLD polynomial constants can be determined from the polarization and CurieWeiss law: 2 1 sP and 1 11 r ( 4 4 ) Table 4 1 Comparison of constant s in e q uation ( 4 3 ) Parameter s This works Literature Units A 5.1684 5.1474 x1010 m 1 5.04 2.012 x109 Nm2/C2 2 12.95 3.608 x109 Nm6/C4 Ps 62.4 75 C/cm2 The comparison between the calculated constant values and literature value for the parameters entering e quation ( 4 3 ) is listed in Table 4 1 Using Table 41 and e q uation ( 4 3 ), we obtain the activation energy and threshold field for the different
98 domain walls. Because the three domain walls have esse ntially identical widths, equation ( 4 3 ) yields almost the same value in activation energy and threshold field. The calculated ac tivation energies are acti~133 mJ/m2 and Eth~1.0x109 V/m with the GLD parameters calculated, and 77 mJ/m2 and Eth~5.0x108 V/m with the constants from the literature. These are comparable to the estimates from the atomistic simulations of of acti MD ~ 255 mJ/m2 and acti DFT ~150 mJ/m2 for Y wall, and acti MD ~ 115 mJ/m2 for XIwall and 140 mJ/m2 for XIIwall and acti DFT ~30 mJ/m2 for an average value of two X walls Because Pz yields essentially the same domain wall width for all three domain walls, the threshold fields are similar. However, the normal and parallel polarization contributions are different in each case, which could have a significant effect on the threshold field. In the absence of a theory of these nonuniaxial polarizations for the trigonal strucutere, we cannot estimate their effect, however. Thus further development of the model for the nonuniaxial component of polarization is necessary to improve the GLD phenomenological analysis. 4. 8 Conclusions and Discuss ion The structure and energetics of domain walls in LiNbO3 has been studied by empirical potentials and verified by DFT. The Y wall is energetically favored but should be less mobile because the activation energy for migration is higher. The broken symmetr ies of both cations and anions near the wall lead to nonuniaxial polarizations. This expands the characteristics of ferroelectric domain wall from simple Ising type to include Blochlike and Nel like type components. The fits to the GLD formalism of uniaxial polarization profiles of different domain walls lead to very similar domain wall widths. Although the absolute value of the inplane polarization is small, less than 3% of
99 uniaxial polarization, the inplane polarization has an important role for understanding ferroelectric domain wall behavior in LiNbO3. Moreover, since other material systems also appear to show such inplane wall components,  it is possible that mixed Bloch like and Nel like character to ferroelectric walls may be ubiquitous. These nonuniaxial polarization components can change the relative electrostatic energies of different wall orientations, and can explain why some wall orientations are seen and others not 
100 CHAPTER 5 DOMAIN WALL DEFECT INTERACTIONS IN FERROELECTRIC LITHIUM NI OBATE 5 .1 Introduction T his chapter focuses on the interaction between de fects and domain walls in LiNbO3 and their effect on domain dynamics. The effect of the type and charge state of a defect or defect pair on domain wall pinning is determined. 5 .2 Methodology and Validation 5 .2.1 Empirical S tudy with T wo D ifferent M ethods There are two widely used methods available to describe the defect formation energies (DFEs) of charged defects: the supercell  method and the Mott Littleton (ML) method  Supercell method: The imposition of periodic boundaries in the supercell method intrinsically involves defect defect interactions, and thus corresponds to a high concentration of defects. Because the supercell method is simply the extension of a bulk calculation, it can be calculated by: ) ( ) ( perfect E defective E Etotal total f ( 5 1) ML method: The ML method embeds a single defect into an atomistic region which, in turn, is embedded in a dielectric continuum. It thus corresponds to the infinite dilute limit of defects DFEs are obtained from the ML method by dividing the surroundings into two different regions; region 1 and region 2a, see figure 5 1 Atoms outside of the regions belong to region 2b. The ions in region 1 are assumed to be strongly perturbed by defects and therefore are relaxed in an explicitly atomistic manner. In contrast, the ions in region 2 are assumed to be weakly perturbed and therefore its relaxation is performed in an approximate way.
101 Figure 51 Three different regions divided by mott l ittleton (ML) method Normally, both the radius of region 1 and the difference of the radii of regions 1 and 2 should be greater than the short range potential cutoff distance to achieve convergence of ML method The ML method considers the total energy of the two region system as the sum of contribut ions of energies of each region and between them: ) ( ) ( ) ( ) (22 12 11 E x E x E x Etotal ( 5 2) where E11( x) represents the energy of region 1 as a function of the Cartesian coordinates, x, E22( ) represents the energy of region 2 as a function of the E12( x, ) is the energy of interaction between the two regions. For the ML method, the DFE is cal culated by the combination of equations ( 5 1) and ( 5 2). Table 51 shows the com parison of DFEs of intrinsic defects between supercell and ML methods. One interstitial site, the empty center of oxygen octahedron, is used for formation energy calculation of interstitial defects. As we can see from T able 51, two different methods yield small difference in the calculated DFEs. The difference in DFEs between our calculation and literature for ML method results from the difference in the size of region 1 and 2.
102 Table 5 1. DFEs calculated by various methods. The interstitials are places at the center of an empty oxygen cage. Supercell (This study) Mott Littleton (This study) Mott Littleton  Lii 7.64 7.09 7.08 Nbi 107.61 106. 51 104.12 Oi 12.64 11.41 9.47 VLi 9.73 9.85 9.81 VNb 122.40 125.32 127.56 VO 18.01 19.26 18.98 NbLi 102.69 102.45 98.37 5 .2.2 DFT with T hermodynamic F rame W orks The DFE of a defect or defects, denoted as with charge state q is defined as  ) ( ) ( ) ( ) ( V E q n perfect E q E P T q Ev F i i i total total f ( 5 3 ) where Etotal( ) is the total energy obtained from DFT calculation of a supercell with the defect(s); Etotal( perfect ) is the total energy of the supercell without any defects. ni is the number of atoms of species i that have been added to ( ni > 0) or removed from ( ni < 0) the supercell when the defects are created; i is the chemical potential of element i T is temperature and P is the oxygen partial pressure. F is the Fermi energy with respect to the valence band maximum in the bulk single crystal. Ev is the valance band maximum of the bulk single crystal system. Therefore, as discussed by Van de Walle and Neugebauer  a correction term V which is the difference in the electrostatic potentials of the defected and un defected systems, is needed to align the band structures.
103 Strictly speaking, the free energy rather than the total energy should be used in e q uation ( 5 3 ) for the calculation of the DFEs. The total internal energies of a supercell obtained from DFT calculations correspond to the Helmholtz free energy at zero temperature, neglecting zeropoint vibrations.  These calculations thus neglect the contributions from the vibrational entropy; fortunately, experimental and theoretical results for entropies of point defects typically fall between 0 to 10k (0 to 0.26 eV at 300 K), where k is the Boltzmann constant  As a consequence, we do not expect this neglect of the entropy term to qualitatively change our conclusions. Detailed analyses by He et al.  and Kohan et al.  also concluded that entr opic effects can be neglected. 5 .2.3 Comparison of DFEs of N eutral D efect C omplexes Charge neutrality requires neutral defect complexes be compared. Here we consider four possible simple defect pairs which maintain the stoichiometry of the LiNbO3 materials: Li Frenkel, Nb Frenkel, O Frenkel pair, and LiNbO3 Schottky complex. However, stoichiometric defect clusters cannot account for the Li deficiency of congruent LiNbO3. Thus, four additional defect complexes which produce nonstoichiometry are also considered: Li2O pseudo Schottky, Nb2O5 pseudoschottky, NbLi +4 VLi ( AntisiteI) and 5 NbLi +4 VN b ( AntisiteII). The purpose of this comparison is to determine the most energetically favorable defect structures. The defect reactions for these complexes are: ( 5 4) ( 5 5) ( 5 6) ( 5 7)
104 For further calculations, three different oxide states, LiNbO3, Li2O and Nb2O5 are considered for the reactions. Nb2O5 has several different structures including T Nb2O5 (Pbam) and B Nb2O5 (C2/c) at highpressure and H Nb2O5 (P2/m) at high temperatures (T > 1100oC)  metastable R Nb2O5 (C2/m), P Nb2O5 (I4122), N Nb2O5 ( C2/m) and M Nb2O5 (I4/mmm). For this study, M Nb2O5 phase is chosen because of its simplicity on crystal structure. Table 5 2 shows structures and the lattice energies of these reference oxides. The s mall differences between our calculations and previously published results  arise from minor differences in the method of calculating the electrostatic energy. Table 5 2 Lattice energy (eV) per unit formula of oxides for empirical study and DFT Space group Empirical DFT This work Reference  Li 2 O 33.15 33.16 14.461 M Nb 2 O 5 314.15 314.37 60.319 LiNbO 3 174.95 174.57 39.552 figure 5 2 compares the DFEs of stoi chiometric and nonstoichiometric defect complexes. The Li Frenkel and AntisiteI pair are the most energetically favorable for both the Li rich condition and Li deficient conditions. For the empirical potentials, it is possible to determine which are the energetically favorable defect complexes by comparing formation energies. A previous study using the ML method concluded that Li Frenkel and AntisiteI clusters are the most energetically favored among the stoichiometric and non stoichiometric defects. This is consistent with the DFT result.
105 Figure 5 2 DFE comparison of (a) stoichiometric and (b) nonstoichiometric defect pairs in LiNbO3 using DFT and empirical method. For nonstoichiometric defect pairs, two different reference states, Nb2O5 and Li2O ca n lead two variants of DFT results. For empirical study, it is not possible to identify different reference state. However, different empirical studies consistently yield same order in DFE with DFT results of Li2O reference state. The supercell method al so leads to similar results, see figure 5 2 Although the precise DFEs calculated by the supercell and ML methods do not match exaclty, both yield the same trends for DFEs. In this study, we will mainly focus on Li Frenkel and AntisiteI defect pairs using the supercell method for DFE calculations. 5 .3 Configurations of T wo D ifferent I ntrinsic D efect P airs 5 .3 .1 Li Frenkel P air ( Lii + VLi ) Figures 5 3 (a) and (b) show two possible structural arrangements of Li Frenkel defect pair. The t wo structures show similar stacking sequences as either two vacancies sit below and two Li ions sit above the center Nb ion ( figure 5 3 (a)) or vice versa ( figure 5 3 (b)). Thus in both cases the cation stacking sequence is Nb V V Nb Li Li rather than NbV Li Nb V Li as in perfect crystal T he differences between the polarization generated by the defect pair (Pd) and the spontaneous polarization (Ps) leads to the difference in energy. The configuration in figure 5 3 (a) creates Pd parallel to Ps, LiF Para while the configuration in figure 5 3 (b) generate Pd antiparallel to Ps, LiF Anti. (a) (b)
106 We also considered the unbounded state at which each Lii and VLi defect exists as an individual point defect without an association effect. The empirical study predicts that the parallel configuration of Li Frenkel ( LiF Para) is energetically more stable state than antiparallel ( LiF Anti) or unbound configuration of the Li Frenkel pair. The association energy is obtained by the difference in energy of arranged defect pair s from nonpaired (unbounded) defects. The association energy of the parallel Frenkel pair (EAsso, Li F Para ~ 0.24 eV) is lower than t he association energy of the antiparallel ( EAsso,Li F Anti ~ 0.02 eV) or unbounded Frenkel pair (EAsso,unbound ~ 0 eV). Figure 5 3(b) shows the comparison of calculated association energies of different Li Frenkel configurations using e mpirical methods. Figure 5 3 Two different configurations of Li Frenkel pair: The polarization created by defect pair (Pd) is (a) parallel and (b) antiparallel to spontaneous polarization (Ps). (c) The comparison on defect formation energies of different Li Frenkel con figurations using empirical method. Supercell methods expect parallel configuration is energetically favorable t han antiparallel configuration. 5 .3 .2 Nb antisite P air(NbLi +4VLi ) The Nb antisite cluster( NbLi +4 VLi ) is another d efect pair which can be found under Li poor conditions. The configuration of NbLi +4 VLi defect pair has previously been proposed by Kim et al  They suggest ed that the lowest energy config uration of the (a) (b) (c)
107 defect pair will show a tetrahedral shape (Fig s. 5 4 a and 5 4 b), since non uniaxial dipole moments created by defect pair may generate non trivial energy (Fig s. 5 4 c though 5 4 e). They also distinguished two different tetrahedral configurati ons depending on the direction of tetrahedron from spontaneous polarization (Ps). They conclude d that the upside down tetrahedral configuration which yields a dipole moment to the bulk polarization ( figure 5 4 a) is energetically more favorable than the ups ide configuration yielding antiparallel dipole moment ( figure 5 4 b), because a dipole moment of Nb antisite pair is created from tetrahedron top to basal plane. Fig ure 5 4 Different structural arrangements of Nb antisite defect pair. Net dipole m oment of defect pair(Pd) are (a) parallel to bulk polarizataion (b) antiparallel to bulk polarization (c) none (d) perpendicular to bulk polarization (e) leaned parallel to bulk polarization. The linear lines are the guidance for the centroid of 4 lithium vacancies. (a) (b) (c) (d) (e)
108 However, a recent DFT study  shows that the formation energy of Nb antisite pair can be lower when all four lithium vacancies are sitting at the nearest Li ion site from the antisite (See, Fig s. 5 4 c through 53e). Also, the DFT study shows that the difference in DFEs among different configurations of four Li vacancies is only 0.11 eV, thus different confi gurations do not lead major difference in defect energetics of Nb antisite pair. Finally, the study concludes that most of the configurations are energetically favorable, so that several arrangements will coexist inside the congruent LiNbO3 material. An e mpirical study with different arrangements of the Nb antisite, shown in figure 5 5 was performed to verify the capability of describing DFE of different configurations of Nb antisite pair. The DFE of the Nb antisite was calculated from: ( 5 8 ) Fig ure 5 5 The comparison of DFE and association energy of different configurations of AntisiteI defect pair using DFT and empirical methods.
109 The lattice energies of antifluorite Li2O ( 33.15 eV) and tetragonal Nb2O5 ( 314.15 eV) was determined, s ee T able 5 2. Figure 5 5 compares the DFEs determined from the DFT and empirical calculations. Consistent with DFT study, the empirical study also didnt show major differences in DFE between different arrangements of Nb antisite pair. Also, Nb antisite pair does not show any preference on the direction of dipole moment with respect to the direction of bulk polarization. The range of association energy between different configurations is from 0.127 eV to 0.165 eV. However, the empirical study predicts pos itive DFEs for all configurations while DFT predicted negative values. The positive DFEs of em pirical study represent that no defect configurations are favorable to form the Nb antisites in nature. This difference originate s from the different reference st ates ; DFT consider the reference state as Li deficient condition while the empirical study just predicts them at stoichiometric condition. Since we are considering Nb antisites at congruent condition, w e focus on the relative differences in DFEs of empiric al study to reduce the error come from the different reference state Although the DFT and empirical calculations are not directly comparable due to the very different reference states, both methods conclude that most of the configurations are equally like ly because of the small differences in formation energy. 5 4 Domain W allintrinsic D efect I nteractions We now turn to the main thrust of this study : the interaction of domain walls with intrinsic defects. We consider two domain wall structures: parallel t o c glide plane ( Y wall) and perpendicular to c glide plane (X wall). We place individual defects or defect clusters at various distances from domain wall structures, which have been previously equilibrated at constant volume. The energetics of the defect/ domain wall interactions are then determined. The interaction between the defects and the domain wall means
110 that their energies cannot be uniquely separated from each other. Here we ascribe the change in the energy of the system to a change in the defect f ormation energy; we could, with equal merit, ascribe this energy change to a change in the domain wall energy. 5 4 .1 Y wall with I ndividual P oint D efects We have determined the energetics of VLi Lii and NbLi as a function of distance from the domain walls. As figure 5 6 (a) shows, the energy of each of these defects decreases as they approach the Y wall. In particular, the energies of the VLi Lii and NbLi defects are 0.23 eV, 0.53 eV and 0.66 eV lower at the domain wall than in the bulk. The effect of the Y wall decreases rapidly as the separation increase, essentially disappearing about 6.5 away. Chapter 4 showed that the Y domain wall can break the symmetry of cations and leads to two inequivalent cation positions. Therefore, we performed DFE calculation of Lii in both sites I and II, see figure 5 6 (b). As figure 5 6 (a) shows there is only a small energy differenc e, 0.01 eV, for the DFE of Lii at the two nonsymmetric positions. Thus, although Sites I and II are structurally distinct, they show very similar energetics. As fiducial systems to validate the calculations with the empirical potentials, we hav e also performed a DFT study on a smaller system size. The DFT calculations involve neutral defects; thus these values cannot be compared directly to the values for charged defects from empirical study. The variation of formation energies determined by DFT are VLi x = 0.11 eV, Lii x = 0.12 eV and NbLi x = 0.14 eV. Both methods predict that the interaction between a point defect and Y wall can lower DFEs of point defects. A similar effect has previous been seen in a study on interaction between an oxygen vacancy and a domain wall in PbTiO3. 
111 Fig ure 5 6 The change in DFE of single point defects near Y wall. (a) The variation in DFE of the most dominant single defects in LNO near Y wall is plotted as a function of distance from the domain walls. (b) Two different interstitial sites caused by broken symmetry of oxygen near the Y wall. For Li interstitial defect, two different symmetric sites, site I and site II are considered for the site dependence of DFE near Y wall. The qualitative agreement between the DFT and atomistic calculati ons provides confidence to move forward with the determination of the energetics of other defects structures. 5 4 .2 Y wall with Li Frenkel P air Fig ure 5 7 Comparison on DEFs of Li Frenkel pair in two different configurations; parallel to bulk polarizat ion(LiF Para), antiparallel to bulk polarization(LiF Anti). Unbound state of Li Frenkel pair is also considered for the comparison. (a) (b)
112 Figure 5 7 shows the variation of the formation energy of two different configurations of Li Frenkel pairs and unpaired stat e of two defects (unbound) near the Y wall using empirical methods. Two different configurations, Pd parallel to PS and Pd antiparallel to PS are considered in order to characterize the effect of structural arrangement of Li Frenkel pairs near Y wall. As w e can see from the figure 5 7 the variation of formation energies between unbounded and bounded Li Frenkel pair shows a similar pattern near the Y wall by lowering the formation energy (0.35 eV for LiF Para, 0.32 eV for LiF Anti and 0.36 eV for unbound) f rom the bulk state. The difference in DFEs between LiF Para and LiF Anti increases from 0.22 eV in the bulk to 0.25 eV near the Y wall. The interaction between the Y wall and Li Frenkel pair disappears at 6.46 away from the Y wall. 5 4 .3 Y wall with Ant isiteI (NbLi + 4VLi ) Because the various configurations do not show much difference in energetics near the Y wall, we only consider two tetrahedr al arrangements Figs. 5 4 (a) and 5 4 (b) among the various configurations of the Anti siteI cluster Fig ure 5 8 Formation energies of AntisiteI clusters as a distance from a Y domain wall.
113 In the first, one of the three sides of basal triangle is parallel to Y wall, which creates a uniaxial component of the dipole moment. Figure 58 sh ows the variation of DFE as a function of distance from Y wall. As in the case of the Li Frenkel pair, the AntisiteI cluster has a lower energy at the Y wall by 0.30eV for unbounded and 0.24 eV for parallel and 0.23 eV for antiparallel. This does not chang e their relative order, with parallel configuration having the lowest energy. 5 4 .4 X walls with P oint D efects Chapter 4 indicates that the X walls perpendicular to c glide plane, display two different variants, XIwall and XIIwall, depending on the inp lane polarization states near the domain walls. In this section therefore, the two X wall variants are considered for the interaction with defect cluster s. As in the case of the Y wall, we initially perform the DFE calculation s with three point defects near two X walls using DFT and empirical method. Both methods predict DFEs are lowered near both XIwall and XIIwalls. However, each X wall show different pattern for interacting with single defects. Again, DFT calculation predicts DFEs of neutral charged defects, so the value cannot be compared directly with DFEs in the empirical study. When determined with the empirical potential, the largest change in the DFE is observed for the NbLi which lowers ~1.45 eV for the XIwall and ~0.77eV for XI Iwall f r o m the bulk value of 102.69 eV. L ii shows smaller changes: ~ 0.71 eV for XIwall and ~0.63 eV for XIIwall from the bulk value of 7.64 eV The smallest domain wall/defect interaction is VLi as ~0.16 eV for both XIwall and XIIwall from the bulk value of 9.73 eV We also perform DFT study to see the variation in DFEs near the X walls. For DFT study, neutral defect s are considered by adding or removing neutral atom for the simplicity of calculation. Thus the change in DFE s near the XIwall and XIIwall are
114 smaller than in the empirical study. The largest decrease is for Lii x with 0.27 eV for XIwall and 0.32 eV for XIIwall from the bulk value of 2.76 eV VLi x shows a smaller variation of 0.14 eV for the XIwall and 0.18 eV for the XIIwall from the bulk value of 2.14 eV NbLi x shows the lowest variation, 0.11 eV for the XIwall and 0.14 eV for the XIIwall from the bulk value of 6.06 eV In contrast to the empirical study, no significant diff erence in the variation of DFEs at both the XIwall and the XIIwall are observed from the DFT study with neutral defects A lthough t he DFEs of two different methods cannot be compared directly, the domain wall defect interactions lead to decrease in the f ormation energy of defects regardless of their charge state. Also, we can predict that higher charged state of NbLi may cause the higher interaction with domain walls. Comparing these values with those from the Y wall, the X walls show a st ronger interaction with defects. Fig ure 5 9 Variation of DFE of point defects near (a) XIwall (b) XIIwall 5 4 .5 X walls with Li F renkel P air We calculate the DFEs of two different Li Frenkel configurations, parallel and antiparallel near the two X w all variants The variations of DFE of Li Frenkel pairs as a function of the distance from X walls are shown on figure 5 10. The DFE of the XIwall (a) (b)
115 decreases by 0.30 eV for LiF Para and 0.32 eV of DFE of LiF Anti. The XIIwall shows a slightly higher inter action, with the DFE of LiF Para lowered by 0.34 eV and the DFE of LiF Anti lowered by 0.36 eV. Thus, both of X walls decrease the formation energies of the Li Frenkels by almost one third from the bulk values. As figure 5 10 shows, the consistent variati on between different configurations explains that LiF Para is always the lowest energy configuration regardless of domain wall defect interactions. Fig ure 5 10. DFE of Li Frenkel defect pair near (a) XIwall (b) XIIwall 5 4 .6 Two X wall s with Antisi teI P air The two X walls also lead to two variants of DFE for the AntisiteI clusters, see Figs. 5 4 (a) and 5 4 (b). For the parallel state, the DFEs of AntisiteI defects are lowered by 0.36 eV for XIwall and by 0.29 eV for XIIwall at the nearest plane, 1.3 away from the X walls Similarly, the DFE of the antiparallel state are lowered by 0.32 eV for XIwall and 0.25 eV for XIIwall. The parallel configuration shows a slightly lower energy (0.07 eV for near XIwall and 0.08 eV for XIIwall) than antiparallel configuration at the nearest plane An interesting phenomenon is observed at 4.5 away from the XIIwall. The DFEs between parallel and antiparallel are switched near XIIwall. This result indicates the possible rearrangement of the defect cluster f rom a frustrate d mode to stable (a) (b)
116 configuration during the domain wall motion, an effect previously discussed by Gopalan et al  Fig ure 5 1 1 DFE of Nb antisite defect pair near (a) XIwall (b) XIIwall 5 5 E xtrinsic D efect Mg 5 .5.1 Reaction S chemes with Mg D opants Donnerberg et al.  demonstrated that cation sites are energetically more favorable than interstitial sites for dopants A recent DFT study  also con sidered Er dopants as substitutional forms A p revious experimental study  on Mg dopants in LiNbO3 show s that most of Mg dopants stays as a substitutional form instead of interstitial. Therefore, we consider five different reactions, each of which maintains the ch arge neutrality of the system. ( 5 9 ) ( 5 10) ( 5 11) ( 5 12) ( 5 13) (a) (b)
117 The energy associated with each reaction is obtained from the formation energy of the products less the formation energies of reactants. Therefore, th e formation energy for each reaction can be shown as below: ( 5 1 4 ) ( 5 1 5 ) ( 5 1 6 ) ( 5 1 7 ) ( 5 1 8 ) Before presenting the results for defect/domai n wall interactions, we perform a study to reproduce the results for bulk defects using the potentials. Table 5 3 shows that the lattice energies of LiNbO3, Li2O and MgO : the very small difference between the results of our calculations using an inhouse c ode and the published results with GULP arise from minor differences in the method of calculating the electrostatic energy. Table 5 3 Lattice energy of oxides Space group Lattice energy per unit formula (eV) This work Reference 41.03 41.04 33.15 33.16 174.95 174.57 The previous calculations used the Mott Littleton approach for calculating the DFEs. Because we focus on defect/domain wall interactions, it is more convenient for us to use a supercell method. We initially perform the calculation of DFE of a single dopant defect in order to compare our resul t with previous calculation by Araujo et al Our l attice energy calculation predict s the formation energy of Mg dopant at different site as 15.47 eV for
118 MgLi and 93.89 eV for MgNb which is the similar value ( 15.36 eV for MgLi and 95.18 eV for MgNb ) from previous Mott Littleton calculation. Table 5 4 compares the DFEs of single dopant calculated with the two methods. Table 5 4 DFE of isolated defect s: comparison of bulk DFEs calculated with a supercell method and the Mott Littleton method. Supercell (This study) Mott Littleton (This study) Mott Littleton 15.47 15. 41 15.36 9 3.89 94.45 95.18 Table 5 5 compares the formation energies of five different reactions. The differences between published and this study come from the different sizes for each region during the ML method. This study considers a large radius, 12 for reg ion 1 and 2a for the convergence of the ML method. As we can see from the Table 5 5, the differences between supercell and ML method decrease with the large radius of each region. Table 5 5 Energies for five reactions, calculated from the DFEs of isolated defects and the bulk lattice energies in Tables 5 2, 5 3 and 5 4 Supercell (This study) Mott Littleton (This study) Mott Littleton Reaction 1 2.17 2.32 2.33 Reaction 2 0.88 1.06 1.38 Reaction 3 5.97 6.62 10.51 Reaction 4 8.96 9.73 14.50 Reaction 5 1.79 2.43 3.01 Reactions 3 and 4 show relatively higher formation energy per dopant than other reactions. A previous experimental study suggests that the dominant defect configuration wil l change from reaction 1 at low concentration of dopant to reaction 1 and 5 at the intermediate dopant level t o reaction 2 at the high defect concentration of
119 Mg dopant. Therefore, we only focus on reactions 1, 2 and 5 for the rest of the calculations. 5 5 2 Domain W all Mg D opant I nteractions Now, we consider the interaction between a single Mg dopant and a Y wall. To develop a basic understanding of defect doma i n wall interactions, we focus on the Y wall. The variation of formation energy of single Mg dopant occupying a Li site or a Nb site near Y wall is shown at figure 5 1 2 The reference states are chosen as the DFE of each dopant at bulk state. As we can see from the empirical study, the formation energy of Mg dopant occupying the Li site increases in the vicinity of Y wall; thus a single Mg dopant will not be trapped by domain wall at Li site. However, a Mg dopant occupying a Nb site has a lower formation energy near the Y wall, and will thus be trapped by the domain wall. Charge neutrality requires t hat there to be defects with compensating charges. The above results suggest that a low energy configuration would consist of MgNb at the domain wall, with the MgLi standing off from the domain wall by ~5. However, the association effect between two Mg dopants also needs to be considered. If the association energy is larger than the energy caused by domainwall/defe ct interaction, MgLi can also sit near the Y wall. Therefore, we performed the simulation, MgNb at the domain wall, with MgLi at varying distance in order to check how the association energy affect the configuration of Mg dopants near the domain wall. The green line of the figure 5 1 2 shows the variation of DFE of the system in which MgLi only varies its position with MgNb at the Y wall. The DFE at the reference state is chosen as the addition of DFE of bulk MgLi and DFE of MgNb at the Y wall. L owest defect formation energy is observed when MgLi sit near the Y wall. Thus, the association
120 energy between MgLi and MgNb is larger than the domain wall /Mg dopant interaction energy 0.05 eV for MgLi and 0.15 eV for MgNb As a result Mg dopants can be exist at both cation sites near the domain walls Fig ure 5 1 2 The variation of DFE of Mg dopant occupying cation sites near Y wall. Positively charged Mg in the Li site show increase in DFE near Y wall which indicate Mg at Li site will not be trapped by domain wall and less effect on domain wall motion and pinning effect. However, negatively charged Mg in the Nb site show the decrease in DFE as other intrinsic defect and defect pair and causing the pinning effect on domain wall motion. For each reaction, we consider the defect as a different cluster shape in the LiNbO3 system. The investigation of the associati on effect caused by different configuration of defect clusters show s the stable configuration of each defect pair near and far away from the Y wall. For comparison, the unbounded state, i n which each defect locates far apart as a single defect without any connectivity between components of defect pair, is also considered. Reaction 1 : Two different arrangements of MgLi compensated by VLi can be distinguished depending on the direction between the dipole moment of the defect pair and the bulk polarization (Ps). Figure 5 1 3 (a) shows the dipole moment of defect pair
121 parallel to bulk Ps; figure 5 1 3 (b) shows the case of the dipole moment a ntiparallel to bulk polarization. Fig ure 5 1 3 Possible arrangement for r eaction 1. (a) dipole moment of defect pair parallel to bulk polarization (b) dipole moment of defect pair antiparallel to bulk polarization Fig ure 5 1 4 The variation of DFE of defect pair using reaction 1 near Y wall. DFE is considered as formation energy per dopant. Although single MgLi shows slightly higher formation energy near the Y wall, the energy compensation by VLi and association effect between two defects cause the Mg dopant to be trapped at the Li site near the Y wall. The association energy for two different configurations are 0.03 eV for antiparallel and 0.04 eV for parallel. Thus, the parallel configuration is slightly more stable than antiparallel configuration. Also, the unbound defect state shows the lowest energy configuration around 4.5 away from (a) (b)
122 the Y wall, which indicate the possibility of the rearrangement of defect pair. The rearrangement of defe ct configuration from frustrated mode to stable configuration during the domain wall motion has been discussed by Gopalan et al.  Reaction 2 : The arrangement of defect pair, MgNb compensated by 3 MgLi can have multiple different arrangements. As we only consider the nearest Li sites from MgNb six possible MgLi sites are shown in figure 5 1 5 (a). Three Li sites among six different sites wil l be replaced by MgLi while the rest will remaine as Li ions. Although we consider the symmetry of the system, this will lead to six distinct configurations. Moreover, the broken symmetry near the domain wall will bring more complexity. In order to simplify the structural configurations, only two configurations : dipole moment parallel or antiparallel to bulk polarization are consi dered for this study, see figure 5 1 5 (b) and (c). Fig ure 5 1 5 Possible arrangement of Mg dopants in LiNbO3 for reaction 2 (a) the position of MgNb and six possible nearest Li sites for the arrangement (b) possible arrangement for parallel dipole moment (c) possible arrangement for antiparallel dipole moment to bulk polarization The variation of format ion energy for Reaction 2 near Y wall is shown on figure 5 1 6 The DFE per dopant predict s that antiparallel configuration i s the most ene rgetically stable arrangement in the bulk state. The association energies of different configurations are 0.27 eV for parallel and 0.42 eV for antiparallel configuration. (a) (b) (c)
123 However, neither the unbound defect pair nor the antiparallel defect pair show s much interaction with the Y wall In particular, the formation energy for antiparallel configuration increases 0.04 eV near the Y wall. By contrast, parallel configuration lowers by 0.13 eV of its formation energy. However, as we can see from the Fig 5 16, the most energetically favorable position for the defect pair, MgNb compensated by 3 MgLi is at t he bulk state with antiparallel configuration. Thus, we can predict Y wall will show less behavior to trap the defect pair, MgNb compensated by 3MgLi Fig ure 5 1 6 The variation of DFE of defect pair with reaction 2 near Y wall. DFE is considered as formation energy per dopant. Fig ure 5 1 7 Possible arrangement of Mg dopants in LiNbO3 for reaction scheme (5). Five Li sites out of six nearest Li sites from Nb vacancy will be occupied by Mg ion. (a) three MgLi generat e parallel dipole moment to bulk polarization and two MgLi create antiparallel dipole moment to bulk polarization (b) two Mg ions at Li site generate antiparallel dipole moment and three Mg ions create parallel dipole moment. (a) (b)
124 Reaction 5: As we can see from the previous figure 5 1 5 (a), the six nearest Li sites are available from one Nb site. By assuming Mg dopants occupy only the nearest Li sites of VNb two different arrangements of five MgLi are possible in the bulk s tate. The two different arrangements can be understood in terms of their relati onship to the bulk polarization: 3 MgLi sit above and 2MgLi sit below the VNb or vice versa. Thus, for the case, in which the 3MgLi is above the VNb the overall parallel dipole moments is parallel to the bulk polarization. figure 5 1 7 shows schematic view of two possible configurations of reaction 5. Figure 5 1 8 shows the comparison of DFE of different configurations for reaction 5. Although the p arallel configuration shows a slightly lower DFE~ 0.06 eV in the bulk region, the interaction with Y wall almost completely removes the preference on configuration for reaction 5. The association energ ies, 0. 74 eV for parall el and 0.68 eV for antiparallel at bulk region change to 0. 6 1 eV for both parallel and antiparallel in the vicinity of Y wall Fig ure 5 1 8 The variation of DFE of defect pair with reaction 5 near Y wall 5 .6 Conclusion Various point defects and defect pairs interact with domain walls in LiNbO3 and lower the energy of the defective system The decrease in the DFE near the domain
125 wall demonstrates that defects and defect pairs can be easily trapped by domain wall. The t rapped defects will generate a loca l cluster field near the domain wall which can be explained to change the ferroelectric and optical properties. A decrease in energy will make the domain wall more stable state and create a higher energy barrier for domain wall motion because the activation energy for domain wall motion is increased by the migration energy for defect pairs. This is domain wall pining. In particular the trapped defects and defect pairs at domain walls will increase the coercive energy and threshold field, which is the pinni ng effect, making the d omain wall motion difficult. A pinning effect affects the ferroelectric and optical properties and causes the ferroelectric degradation known as fatigue.  S uccessful growth of stoichiometric LiNbO3 can prevent the pinning effect by defects, seen in congruent LiNbO3, and improve the stability of LiNbO3. The comparison on the energetics between intrinsic defects and Mg dopants explains the optical stability and low coercive field of Mg doped LiNbO3. When we compare these values ( 0.05 for MgLi and 0.15 for MgNb ) with th at of intrinsic defects (0.23 eV for VLi 0.57 eV for Lii and 0.66 for NbLi ) Mg dopants show much less interaction with the Y wall and are less trapped by Y wall. The decrease in the trapped defects near the domain walls improves the optical birefringence and the domain wall motion. Thus, Mg doped LiNbO3 yields better optical properties and a reduction in coercive field after doping with Mg.
126 CHAPTER 6 DOMAIN WALLS IN LITHIUM TANTALATE 6.1 Introduction L iTaO3 is a ferroelectric material with trigonal structure, similar to LiNbO3. It has a high Curie temperature (~940 K)  and spontaneous polarization (~50 2)  Because of its unique optical properties, LiTaO3 is a candidate for non linear frequency convertor s  electro optics  and holography  The control of shape of the domain is a key for optical applications. D ifferent domain shapes can focus or defocus light in the plane and convert the input wave length to different wave lengths. Congruent LiTaO3 sh ows lithium(Li) deficient composition with CLi/[CLi+CTa]~0.485  Thus congruent LiTaO3 shows a large difference in Curie temperature from stoichiometric LiTaO3  This non stoichiometry of congruent LiTaO3 can cause drastic changes in optical, electric and elastic properties near domain walls and affect the domain dynamics.  Recent studies [166, 167] show that the threshold field for domain reversal c an be changed depending on the stoichiometry. Moreover, the crystallographic plane and shape of stable domain wall s are changed from Y wall, parallel to c glide plane, to X wall, perpendicular to c glide plane. Figure 61 shows different domain wall struct ures in LiNbO3 and LiTaO3. At stoichiometric conditions, b oth form hexagonal shape domain wall structures of Y wall, see figure 6 1(a) and (b). As their compositions are change to Li deficient condition for congruent system, however, they show different do main wall structures, see figure 6 1(c) and (d). LiTaO3 changes its stable domain wall structure to X wall with a triangular domain shape, but LiNbO3 maintains the hexagonal shape domain with Y walls Gopalan et al. explained the hexagonal domain shapes us ing anisotropy of piezoelectric effect s  The change
127 of domain wall shape from hexagonal to triangular under different electromechanical coupling  have been studied by Shur. However, no systematic understanding of the role of defects for the change of the domain wall shape in congruent LiTaO3 has developed. In this study, therefore, we focus on understanding the variation of the domain shape fro m the study of the energetics of the domain walls and intrinsic defects in LiTaO3. Fi gure 61. Different shapes of domains in LiNbO3 and LiTaO3: (a) stoichiometric LiNbO3, (b) stoichiometric LiTaO3,  (c) congruent LiNbO3, (d) congruent LiTaO3.  Optical images: (a) domain revealed by etching, (b) phasecontrast microscopy, (c) and (d) piezoelectric force microscopy phase contrast images. (a) (b) (c) (d)
128 6.2 Potential D evelopment As discussed for LiNbO3, f or systematic study such of domain wall s, we use empirical method because of the size limitation on density functional study (DFT) methods Since no potential is available for LiTaO3 system, we developed a new potential by simultaneously fitting to the structure and elastic properties  of LiTaO3. The General utility lattice program (GULP) was used for the fitting. In fitting we retain the Li O and O O potential parameters developed by Jackson et al.  This allows additional study for composite m aterials between LiNbO3 and LiTaO3. The full details for potential parameters are discussed in Chapter 2. Table 6 1 shows comparison between the potentials of LiNbO3 and LiTaO3 Table 61 The comparison o f potential parameters between LiNbO3 and LiTaO3 L iNbO 3 LiTaO 3 Interaction A (eV) () C (eV/ 6 ) A (eV) () C (eV/ 6 ) Li O 950.0 0.2610 0.0 0 950.0 0 0.2610 0.0 0 Nb(Ta) O 1425.0 0.3650 0.0 0 1131.56 0.3845 0.0 0 O O 22764.0 0.1490 27.88 22764.0 0 0.1490 27.88 LiNbO 3 LiTaO 3 Shell parameter Shell charge, Y Spring constant, K 2 (eV/ 2 ) Shell charge, Y Spring constant, K 2 (eV/ 2 ) O2 2.9 70.0 2.9 70.0 LiNbO 3 LiTaO 3 Three body parameter Force constant, k (eV/rad) Equilibrium angle, 0 (deg) Force constant, k (eV/rad) Equilibrium angle, 0 ( deg) O Nb(Ta) O 0.5776 90 0.1405 90 The two body interaction parameters, A, and C between Nb and O are changed for Ta O. Similarly, the three body interaction parameter, k for O Nb O are changed for the O Ta O interaction. 6.3 Validation of LiTaO3 P otential The six oxygen octahedral cages in LiTaO3 are continuously stacked with shared close packed planes as shown in figure 6 2 Similar to LiNbO3, the octahedral is filled in
129 the sequence of Li, Ta, vacant, Li, Ta, vacant a long the z axis of trigonal structure. As the polarization state switches the sequence can be change to Li, vacant, Ta, Li, vacant, Ta. Figure 62. Crystal structure of LiTaO3 for (a) +Ps and (b) Ps in z direction. The trigonal structure with rhombohedral axes can be described by hexagonal lattice parameter s a=5.154 and c=13.783 in the six formula unit hexagonal cell [92, 165] The structural coordinates were obtained by least squares fitting method of neutron scattering data The elastic constants were measured by Takanaga and coworker s We have compar ed of our simulation results using new LiTaO3 potential with structural coordinates and elastic constants, see T able 62 and figure 6 3. Our new inter atomic potential reproduces the structural parameters extremely well, the l arge st error bar only 0.4% for the lattice parameter (a) (b)
130 Table 62 Comparison of structural parameters and elastic properties of LiTaO3 between experiments and simulations Parameters Simulations Experiments % Error Lattice parameter a 5.171499 5.154 0.3 4 c 13.781961 13.783 0.008 Fractional atomic coordinate z (Li) 0 .276429 0.279 0.003 z (Ta) 0 0 0 x (O) 0.050146 0.0501 0.00005 y(O) 0.352944 0.3436 0.009 z(O) 0.068474 0.0687 0.0002 Figure 63. Comparison on elastic constants between experiment and simulation. The elastic constants are important material parameters for mechanical applications. The simulations reproduce elastic properties of LiTaO3 within acceptable limits. The C11 and C33 are the most important elastic constants for reproducing pure
131 longitudinal mechanical properties of trigonal materials. The errors from our simulation are only 22.3% for C11 and 17.75% for C33. The elastic constants, C44 and C66 reproduce pure shear mechanical properties. Our simulation study also predicts them well with the 34.64% error for C44 and the 8.70% error for C66. The elastic constants, C12, C13 and C14 are relatively less important for mechanical application. Figure 64 The comparison of normalized polarization between simulation and experiment  The temperature dependence of the polarization was also determined using molecular dynamics simulations with the new inter atomic poten tial. Figure 6 4 shows the variation of normalized polarization as a function of temperature. Although a difference in polarization profile at high temperature is observed, both experiment  and simulation results show gradual changes in polarization at low temperature and a relatively abrupt change of polarization near Curie temperature. Our simulation yield the
132 Curie temperature as 1000K which matches quite well with ex perimental value, of ~974K  for stoichiometric LiTaO3. Therefore, we have confidence that the new potential will be able to discuss the energetics of domain wall and defects. 6.4 Energetics of D omain W alls in L ithium Tantalate Similar to LiNbO3, LiTaO3 also has two different domain walls; Y walls, parallel to cglide plane and X walls, perpendicular to c glide plane. Simulation studies are performed in order to determine the energetics of the two domain walls us ing both DFT calculations and empirical method with our potential. Both methods consistently show that the Y wall sits between oxygen planes and X wall sits between atom planes, see figure 6 5. Because the DFT method is known to be more accurate than the e mpirical method, the comparison between two methods is a very good way to validate our potential. DFT calculation yields a Y wall energy of 60 mJ/m2 when the center of Y wall is sitting between two anion planes. It is also possible to determine the Y wall energy at transition state by fixing the wall at the atom plane. For the Y wall at the cation plane, the maximum domain wall energy of 78 mJ/m2 is observed. Assuming that this is the maximum in energy, these calculations yield a Peierls barrier to motion of 18 mJ/m2  The corresponding analysis using the empirical potential yield an minimum energy of 70 mJ/m2 and maximum energy of 113 mJ/m2 with two different locat ions of the Y wall. Thus the empirical study with our potential yields a Peierls barrier of 43 mJ/m2, which is tw ice larger value than DFT study. The Peierls barrier is typically large, one. We shall see that the domain wall widths are quite narrow, is consistent with the Peierls potential.
133 The two methods also agree on the equilibrium position of the X wall. The lowest energy for the X wall is observed at the center between the ion planes. DFT calculation yields an X wall energy of 63 mJ/m2 and empirical method yields 80 mJ/m2. We also perform the calculation of energy by fixing the X wall position to atom plane. The DFT calculation yields a X wall energy of 75 mJ/m2 at atom planes, and empirical method yields 97 mJ/m2. Thus the DFT method yields a Peierls barrier of 12 mJ/m2 and empirical study with our potential yields a Peierls barrier of 17 mJ/m2. Figure 65 Domain wall energy comparison in LiTaO3. (a) Y wall at cation plane and between anion planes (b) X wall at cationanion mixed plane and between mixed planes. Although, as is to be expected, the DFT and empirical values do not precisely agree, both methods do predict the Y wall shows slightly lower energy t han X wall, 3 mJ/m2 for DFT method and 10 mJ/m2 for empirical study. This conclusion is consistent (a) (b)
134 with the conclusions of Scrymgeour and his coworkers  based on calculations using Ginzburg LandauDevonshire theory. The comparison on the energetics between LiNbO3 and LiTaO3 will improve our understanding on the properties of the domain walls. Figure 66 compares the domain wall energetics between LiNbO3 and LiTaO3. The domain wall energies of LiTaO3 are only one third of that of LiNbO3. Higher domain wall energy of LiNbO3 can be understood by difference in the ionic radii of a Ta ion and a Nb ion. The smaller ionic r adius of the Nb ion allows it to displace a larger distance and leads to a larger lattice mismatch with higher piezoelectric and electrostrictive constants.  Thus, LiNbO3 shows higher domain wall energies with larger piezoelectric effect. The lower domain wall energy can also explain the relativel y lower coercive field for LiTaO3 (15~17 kV/cm)  than that of LiNbO3 (33.14~40.54 KV/cm)  Figure 66. Domain wall energy comparison between LiNbO3 (LNO) and LiTaO3 (LTO).
135 As discussed i n Sec. 6.1, the domain wall structures between LiTaO3 and LiNbO3 are significantly different at Li deficient condition. This can be understood by the difference in domain wall energetics between LiTaO3 and LiNbO3. For LiNbO3, significant energy differences between Y wall and X wall are observed from both of MD and DFT methods; Y wall is energetically lower by 40 mJ/m2 fo r DFT and 30 mJ/m2 for empirical study than X wall. For LiTaO3, however, the only minor differences in domain wall energies are expected between Y wall and X wall from both methods; the difference is only 3 mJ/m2 for DFT and 10 mJ/m2 for empirical study. A s we discussed on Chapter 5, the X wall has a relative higher tendency to low er its energy during interaction with antisite defects that that of Y wall. As we will discuss on Sec 6 6, TaLi compensated by four VLi defects are the major defect pair in congruent LiTaO3. With such small differences in domain wall energy of LiTaO3, the preferred domain walls are easily switched by defects domain walls interactions. Thus, the shift of preferred domain wall orientation from Y wall to X wall is observed from LiTaO3. The change in the shape of domain walls can thus be explained. Figure 67 shows two possible shapes of X walls: triangle and hexagon. The white region represents the up polarization and shaded region represents down polarization region. The arrow represents the orientation of displacements of three different oxygens near the X walls. For triangular domain shape with X walls, all three X walls, aligned along three edges of the triangle have the same origin of oxygen displacements, see figure 6 7(a). In this particular case, XIwalls are created along each edge of the triangle. For the opposite case, three XIIwalls are expected. For the hexagonal shape, however, additional edges, which are aligned with oppositely orientated X walls, are expected, see figure 6 7(b). Because the
136 antisite defect show high er in teraction with especially XIwall to lower its formation energy, additional three XIIwall edges are less preferred. Therefore, the hexagonal shapes of th e domain walls are less likely for congruent LiTaO3. Th us XIwalls with the triangular shapes are observed for congruent LiTaO3. This is also consistent with experimental preference of XIwall for congruent LiTaO3.  Figure 67. The schematic view of two possible shapes of X walls: (a) triangle and (b)hexagon. 6.5 Domain W all S tructure at T ransition S tate In Sec 6 4, we discussed the Peierl potentials and the two structurally distinct positions of domain walls. The Y wall is sit s on either the cation plane or between two anion planes. Similarly, the X wall sits on either cationanion mixed plane or between mixed p lanes In Chapter 4, we discussed the structural change near the domain walls at optimized states. Intuitively, we can expect that the atomic structure at the domain wall plane will be the average of the structures on the opposite sides. In this sec tion, w e discuss the structural differences in the Y wall and X wall at their transition states. (a) (b)
137 Figure 68. Structural arrangement of cations on (a) the Y wall plane and (b) the X wall plane at the transition state. Figure 68 shows the domain wall planes o f LiTaO3 at the transition state s. For both domain walls, Ta atom sit at the center of the oxygen cage, which is their position in paraelectric state. However, the walls show two different behaviors for the Li atoms on the atom plane. For the Y wall, the Li atom moves above and below from the oxygen plane in hexagonal Z direction, to creates a stack of local dipole moment on the domain wall. In other words, the structure of Li atoms on Y wall plane looks like a mixture of the up and down polarization. Howev er, the local dipole moments are cancelled out by equal and opposite motion of Li atoms, so the Y wall leads net zero polarization on the domain wall plane. For the X wall, however, Li atom sits at the exact center of the oxygen planes which can be underst ood as the paraelectric structural position. The positional change of Li atoms is more continuous through the X wall. Thus, any of the (a) (b)
138 local z directional dipole moments can be observed from the X w all planes. Although neither domain wall shows a net polar ization in domain wall plane, Y wall plane shows a mixture of the up and down ferroelectric structures while the X wall has the paraelectric structure. We also observed similar patterns for LiNbO3. 6.6 Defects in Lithium Tantalate Defect formation energy ( DFE) calculation s are also performed using the Mott Littleton (ML) method with our developed potential. In Table 63, the DFEs are reported for possible intrinsic defects in the structure. The interstitial site is chosen at fractional coordinate (0.0, 0.0, 0.13936) which corresponds to the center of the empty oxygen octahedron of the structure. To understand the effect of temperature for DFEs, DFE calculations are performed at two different temperatures (0K and 300K). The simulation results show very similar values with minor difference in DFEs at two different temperatures. Table 63 Intrinsic defect formation energies. Formation energy per defect (eV) 0K 300K Li i 7.48 7.56 Ta i 104.82 105.50 O i 11.39 11.75 V Li 9.50 9.26 V Ta 120.87 120.00 V O 19.44 19.28 Ta Li 100.94 101.42 We should note that the individual point DFEs do not have any direct significance unless they are combined in charge neutral defect clusters. Therefore, the combination
139 of those defects to form the Schottky, Frenkel, or pseudoSchottky defect pairs are considered in Table 63. Pseudo Schottky defect pair consists of two vacancies among three components of LiTa O3. Two different oxide states, Li2O and Ta2O5, are considered as the reference state for nonstoichiometric defect pair, see Table 65. As we can see from the Table 64, the Li Frenkel pair is the most stable defect pair among different configurations. Table 6 4 Frenkel, Schottky and pseudoSchottky defect formation energies (eV/defect) (association effects are not considered between defects) Formation energy per defect (eV/defect) 0K 300K Li Frenkel 1.01 0.85 Ta Frenkel 8.02 7.25 O Frenkel 4.02 3.76 Schottky LiTaO 3 4.81 4.36 Pseudo Schottky Li 2 O 1.76 1.64 Pseudo Schottky Ta 2 O 5 3.84 3.52 Table 65 Lattice energies of different oxide states Space group Lattce energies (eV) 0K 300K Li2O Fm 3 m 33.15 32.89 Ta 2 O 5 P6/mmm 312.07 311.77 LiTaO 3 R3c 174.26 174.02 Similar to LiNbO3 study, complex defect models are also considered in order to explain the experimentally observed lithium deficiency condition for congruent LiTaO3. Three different models which can explain the congruent growth condition of LiTaO3 are considered; Antisite Ta compensated by Li vacancies, Antisite Ta compensated by Ta
140 vacancies, and Interstitial Ta compensated by Li vacancies. Table 66 shows the energy comparison between different defect models. Our cal culation predicts Antisite Ta compensated by four Li vacancies to be the most energetically favorable defect pairs among three different models. Table 66 Energies of possible reactions that lead to Li2O deficiency Formation energy per Li 2 O unit (eV) 0K 300K Ta Li +4 V Li 2.05 1.86 5Ta Li +4 V Ta 5.78 5.17 Ta i +5 V Li 2.64 2.41 6 .7 C onclusion A new LiTaO3 potential is developed and shows very good agreement with study experimental properties on LiTaO3. Also, th e potential gives very good agreement with DFT results on domain wall energetics and defect formation energies. The comparison of domain wall energetics and structures between LiNbO3 and LiTaO3 shows the two systems can lead very different domain wall ener getics in spite of the similarity in their structures and properties. LiTaO3 shows a domain wall energy of only onethird of that in LiNbO3 which explains the difference between two systems. Preferred orientation and shape of domain walls in LiTaO3 and LiN bO3 are explained by the differences in the energetics of domain walls and defects domain wall interactions. For LiNbO3, the domain wall orientation remains oriented along crystallographic y axes in congruent as well as stoichiometric LiNbO3 crystals, since X wall yields much higher domain wall energy than Y wall. However, small energy difference between two domain walls in LiTaO3 can be overcome by defect domain wall interactions. The stable domain wall
141 structure of LiTaO3 is shifted from Y wall with hexag onal shape to X wall with triangular shape during the change in stoichiometry.
142 CHAPTER 7 SURFACE DEPOSITION O N POLYTHIOPHENE MOLECULES 7 .1 Introduction Surface polymerization by ionassist ed deposition (SPIAD) enables the control of thin film chemistry a nd morphology on the nanos cale during growth of conducting polymer thin films. Because SPIAD can be performed by a mass selected beam of thiophene, it can simplify the growth mechanism by eliminating all nonmass selected sources such as fragment ions, radicals and protons. The SPIAD of thiophene incidents has been carried out by both experimental and in a density functional theory molecular dynamics (DFT MD) study  Previous experimental work on SPIAD found that the polymerization proces s can only occur at specific conditions such as specific ion/neutral ratios and energies of incident s [105, 115] T he previous DFT MD study considered t he deposition of mass selected beams of thiophene at 100 t o 500 eV/molecule energy range with thin film substrate of terth iophene (3T) oligomers The study predicted that the terthiophenes on the substrate are mainly damaged by incident thiophenes at the energy range  Thus, further study is necessary for optimizing the SPIAD process for polythiophene. As a n extension of the previous study, th e SPIAD of thiophene incidents is performed with relatively lower incident energy, 25 eV. We also perform a study of the surface chemistry between a nanocrystalline quantum dot and t hiophene polymer chains. PbS and PbSe nanocrystal photovoltaics devices are receiving concentrated att ention for their potential usage as a solar cell [176, 177] Especially, the combination between PbS nanocrystalline and organic materials, so called organic photovoltaics has been widely studied because of the low production cost and easy roll to roll process ing [178 182] M ultijunction photovoltaics made by a
143 pair of near infrared and visible active layers could theoretically improve the solar to electrical conversion efficiency by up to 50% High electron affinity PbS nanocrystals mainly act as an infrared active layer for multijunction photovoltaics while the polymer layer acts as a visual active layer, leading to enhanc ed charge transfer. However, most multijunction photovoltaics are limited in their efficiency to 10 ~ 20% which is much below their theoretical limit [183, 184] Recently, a nonheterojunction high efficiency ( 50%) photovoltaic device has been suggested using the property of the multiple excition generation of PbS nanocryst als [185, 186] Unlike traditional photovoltaic devices which generate one excition for each observed photon, multiple excitions generated by one photon can improve the effi ciency of photovoltaic device to up to 50% of solar energy. However, solar to electrical conversion efficiency is still below 10% because of surface chemistry, charge transfer efficiency, lifetimes of the excitons and other phenomena. [187 189] In this study, therefore, we focus on elucidating the detailed mechanism for the growth and polymerization of polythiophene and the charge transfer mechanism between PbS nanocrystalline and polythiophene layer to improve transferability of excitons gener ated by photons. 7.2 SPIAD The mass selected SPIAD technique is modeled by using a computer simulation tool, Spanish Initiative for Electronic Simulations with Thousands of Atoms (SIESTA) As we discussed on Chapter 2, SIESTA allows us to examine both the electronic structure and the kinetics of molecules using the DFT MD methodology. Also, the linear scaling with system size of SIESTA allows us to simulate larger supercell system s by impro ving the high computational cost of plane wave DFT MD simulations. I nitially, the basis set for Gaussian type electron orbitals is optimized for carbon, hydrogen and
144 sulfur system In collaborati o n with Prof. Julian Gale of the Curtin University of Australia, one of the coauthor s of SIESTA, we successfully optimized the n ecessary basis sets. Figure 71 Comparison between hybridized HF calculation and linear scaled DFT calculation. Two different methods yield slight difference in bonding energies and structures of 3T polymer chain. Figure 71 shows the comparison betw een quantum chemical calculation and linearly scaled DFT study for the structure of terthiophene. The B3LYP method [190 193] with correlationconsistent polarized valence quadruple zeta (cc pVQZ) basis sets is used for hybridized quantum chemical calculation and PBE pseudopotential with DZP basis sets is used for DFT calculation. The energy on the left side of figure 7 1 represents the bonding energy of 3T. The DFT study with optimized basis sets estimate approximately 4 eV stronger bonding than the B3LYP calculation. The difference corresponds to a 3.3% difference in bonding energy of 3T polymer chain. The optimized structures are compared on right side. DFT overestimate the bonding lengths of C C and C H, but underestimate the bonding lengths of C S. Also, the two methods show (a) (b)
145 small differences in the optimized bonding angle. However, the quantum chemical method requires much larger memory with a longer time scale; the quantum chemical method is approximately 200 tim es slower than that of DFT calculation. Taking into account the computational costs for the calculation, however, the linear scaled DFT efficiently predict reasonable results within a relatively small sacrifice on accuracy. Figure 7 2 The schematics o f two structures with different system sizes. (a) Two layered structure with four terthiophenes (3T); Two 3Ts are parallel to each other in each layer. The backbone chains of 3Ts in first layer are perpendicular to that of second layer. Each 3T is separated by 3.5 (b) Five layered structure of thiophene oligomers; Three 3Ts and two dithiophenes are combined in a layer. The backbone chain of each layer is angled 45 upper layer. The distance between each layer is about 3.5 To model the SPIAD proc ess, we create substrate structures of two different system sizes: a small simulation cell with only 3T molecules and a large simulation cell with a combination of dithiophene(2T) and terthiophene (3T), see figure 7 2 The large simulation cell includes five layer s of polythiophene with a total of 505 atoms. Thus, each layer has 101 atoms combination of three 3T and two 2T polythiophenes. The small cell includes 92 atoms with two layer s with two 3T in each layer. Each polythiophene is separated by 3.5 wit hin a layer and between the layers. For the small (a) ( b)
146 system the density of the substrates are chosen as 1.667 g/cm3 which is similar to the value, 1.503 g/cm3, examined from X ray study on 3T crystalline. For the large substrate, a relatively smaller density thiophene structure, 1.034 g/cm3 is chosen to allow the study of the effect of density on reactivity of polythiophene as well as the effect of consecutive deposition. Figure 73 The angle between incident thiophene ring and substrate 3T ring is determ ined in order to study angle dependency of reaction. Probability chart shows that breaking C S bonds are less correlated with angles of incident molecules from substrate 3T. Because incident thiophenes with kinetic energy of 35 eV or higher cause da mage to polythiophene and form debris after the collision, i nitial energy, 25 eV is used for incident thiophene molecules for deposition. 10 0 different trajectories of single thiophene deposition are examined f or the small simulation cell. About 25 different t rajectories (2 5 % among single deposition trajectories) lead to breaking carbon(C) sulfur(S) bonds of an incident thiophene or 3Ts in substrate. The C and S ions from the broken sites act as active particles for reaction and actively lead to either polymeri zation or cross linking to the incident molecule or to another 3T. Depending on the orientation and deposition sites of incident molecules, different dynamic behaviors are observed. Figure 73 shows the probability of breaking the C S bonds based on the angle between (a) (b)
147 the incident thiophene and the substrate 3T. The angle between two thiophene rings is chosen for the measuring, see figure 7 3(a). Each trajectory is distinguished depending on the angle and split into a 10 degree bin. As we can see from the fi gure 7 3(b), no general rule between breaking C S bonds and the angles increase of the thiophenes is observed; Only one trajectory among 10 attempts successfully lead to the breaking of the C S bonding at the angle between 0 reactions between incident and substrate 3T are strongly correlated with the collisio n sites. Figure 74 through 77 show f our different trajectories observed from the study. The first trajectory, in which an inc ident thiophene is deposited into the side of a 3T molecule, forms a polymerization between a 3T and an imag e 3T in substrate, se e figure 7 4 With breaking C S bonds from the incident thiophene, two 3Ts create consecutive polymer chain with a bridge of the incident thiophene. This reaction generates a polymer fragment with a molecular weight, 576 g/mol. Another trajectory shows cro sslinking between two 3Ts at different depths in the substrate, see figure 7 5 The incident thiophene breaks a C S bond by pushing away the sulfur in 3T polymer chain. The carbon ion from the broken C S bonding actively creates the bonding with carbon fr om another 3T molecule. This reaction leads to a final product with molecular weight, 492 g/mol which is equivalent to the molecular weight of two 3Ts. The reaction between an incident thiophene and a 3T molecule is also observed from figure 7 6 This reac tion is obtained, when the C S bonds with the incident thiophene is broken during the impact with the substrate. The incident molecule with a broken chain immediately reacts with a substrate 3T to generate a branching polymer fragment The molecular
148 weight of final product is 330 g/mol. The f ourth reaction shown in figure 7 7 breaks the C S bon ds in the center of 3T molecule. In this case, however, the broken 3T polymer chain fails to create a stable C C bonding with either incident thiophene or other 3Ts. Therefore, the broken 3T polymer c hain remains an active state. The molecular weights of the final products remain unchanged Therefore, polymerization or crosslinking under 25 eV deposition of a thiophene can be understood as two step processes; breaking C S bonding and creating stable C C bonding between polythiophenes. Figure 74. Snapshots at different time from DFT MD simulations of neutral system. Surface polymerization and cross linkings are observed from depositing of a thiophene molecules with 25 eV on a terthiophene oligomer thin film; polymerization between two 3Ts and an incident thiophene. Figure 78 shows the distribution of molecular weight among successful trajectories leaded to broken C S bonds. 13 successful trajectories, which successfully lead to new molecular weight products, are only observed among 100 trials. Even though it is difficult to determine probabilities with small number of successful trajectories, we observe 10 successful trajectories for [ 3 T]C4H4S resulting from third reaction, two trajectories for [2T]C4H4S[3T] and one case for [3T]C4H4S[3T]. Our study successfully
149 reproduces different molecular weight products from surface polymerization and cross linking processes. Figure 75. Snapshots at different time from DFT MD simulations of neutral system. Surface polymerization and cross linkings are observed from depositing of a thiophene molecules with 25 eV on a terthiophene oligomer thin film; crosslinking between two 3Ts. Figure 76. Snapshots at different time from DFT MD simulations of neutral system. Surface polymerization and cross linkings are observed from depositing of a thiophene molecules with 25 eV on a terthiophene oligomer thin film; polymerization between a 3T and an incident thiophene.
150 Figure 77 Sna pshots at different time from DFT MD simulations of neutral system. Surface polymerization and cross linkings are observed from depositing of a thiophene molecules with 25 eV on a terthiophene oligomer thi n film; broken 3T chain and an incident thiophene. Figure 7 8 The distribution of molecular weight for four different trajectories which lead successful reaction with 3Ts in the substrate
151 Because of the very much larger computational cost, only 5 trajectories of single thiophene deposition are perf ormed with 25 eV of incident thiophene in large unitcell s. As we observed from small system size results, surface polymerization and cross linking of 3T with incident thiophene of 25 eV is a relatively low probability event. Moreover, the lower density of the large system also causes a lower reaction rate between polythiophenes. As a result no polymerization or cross linking reaction is observed from single deposition trajectories in a large simulation cell. We are currently running consecutive deposition in large unitcell system to characterize accumulated damage. 7.3 C2H+ ion deposition on 3T polymer with PbS nanocrystalline As we address ed in Sec 7 1, the combinations of PbS nanocrystalline with polythiophenes has a potential applications for organic photovoltaics. The low electric conductivity between PbS nanoparticle and polythiophenes limits its efficiency for an organic solar cell. This limitation can be overcome with a change in surface morphology between thiophene and PbS nanocystalline. We perfor m the simulation study by depositing C2H+ ion on the substrate, in which a PbS nanocluster is surround by 2Ts and 3Ts. The positive charge was place on the system as a whole and was determined to be located on the substrate An incident energy of 100 eV ki netic energy is considered for the incident C2H+ ion. Our simulation study shows that a bond between PbS nanocluster and a polythiophene is created by a process similar to the growth mechanism, discussed in previous section. Figure 79 shows the iso surfac e of electronic charge density of DFT MD. I nitially, no bo nding between terthiophene and PbS nanocrystalline is observed, as indicated in figure 7 9 (a). As C2H+ molecules are deposited into the mixed substrate of polythiophenes and PbS nanocrystalline, ind uced
152 bonding between terthiophene and PbS nanocrystalline is predicted to occur In particular t he incident molecule breaks the C S bonding of 2T and generates the C2TCin bonding with it. The iso surface of C2TCin bonding is clearly connected to the isosurface of PbS nanocrystalline, see figure 7 9 (b). Obviously, 2T creates electronic connection to the PbS nanocrystlline. Therefore, the active polymer chain generated from broken C S bond of polythiophene can introduce electronic connectivity between thio phene oligomer and PbS quantum dots. Figure 7 9 The iso surface of electron density. The C2TCin bonding between incident particle, C2H+ and 2T at the site of the broken C S bonding leads the connection to PbS nanocrystalline. C2T and Cin represent the carbon from 2T and incident particle. 7.4 Conclusion s The growth and cross linking of polythiophene is simulated using SPIAD method. M any different factors including the kinds of incident particles, collision angles and positions, and kinetic energy o f i ncident particle affect the reactivity during SPIAD. In this study, the possible reaction mechanisms of thiophene incidents with polythiophenes at relatively low kinetic energy collision (25 eV) are demonstrated. The growth and (a) (b)
153 polymerization mechanism are observed as two step processes; breaking C S bonding of thiophene and creating C C bonding between broken bonds. We categorize four different reactions for the growth of thiophene chains. The study on surface reaction between PbS quantum dots and polythiophenes also confirms that broken C S bonds are a key for further reaction of polythiophenes. We are currently performing simulations with larger cell in order to understand the effect on polymerization and cross linking processes of polythiophenes from consecutive depositions of thiophenes.
154 CHAPTER 8 ADATIVE KINETIC MONTE CARLO (AKMC) 8 .1 Introduction Classical MD simulations are limited to short periods of time almost always less than a microsecond This is because t he time step for MD simulations is in the order of a femtosecond. Thus a simulation of 1 microsecond require approximately 1 billion steps, which is computationally very expensive. S everal mechanisms have been proposed to bridg e the time scale gap ; parallel replica dynamics  t emperature a ccelerated dynamics  (TAD) and k inetic Monte Carlo [196, 197] ( kM C ). Parallel replica dynamics provides a linear acceleration by running a replica of a state in each computing node and accepting the first escape.  TAD accelerates dynamics by performing MD simulation at elevated temperature. It obtains the trajectory of dynamics by neglect ing nonphysical behavior at an interesting temperature. TAD can reproduce true dynamics well, since the approximations are well controlled in. However, it still requires huge amount s of computing time to reproduce rare events KMC is an alternative method. The kMC does not attempt to reproduce the dynamic trajectory. Instead, it calculates a list of possible transition events, and then chooses to execute one transition at random The rate of each event will determine the probability of each event in the kMC The list of transition events is then updated according to the current configuration. A trajectory for the system evolution is obtained by repeating this process. One problem with kMC is t hat it usually requires approximations to generate a list of possible escape paths.  A second problem is that the possible escape mechanism has to be known in advance. If the specified mechanism is not intuitiv e it will be very difficult to obtain all possible transition states. The adaptive kinetic monte c arlo ( akMC )
155 method can re solv e this issue. In akMC the rate for each of the events can be calculated accurately from transition state theory   The list of events for the system at specific condition can be obtained using either classical interatomic potentials or even first principles calculations   I n th is chapter, I will discuss the mechanism of akMC and how this concept is implemented into a code. I will test the code for the dynamic behavior of aluminum atom on an aluminum (111) surface. 8.2 Mechanism 8.2 .1 Kinetic Monte Carlo ( KMC ) The KMC algorithm simulates the time evolution of a system when processes can occur with known rates, r. At t=0, the system is in state i. Then, all kinetic events that are available at the state i are listed with the rate for the each event. We now calculate the total rate, Rtot for all events at state i using the cumulative function. N j j totr R1 for N j 2 1 (8 1) where N is the total number of transitions at the state i. Then, one process among the possible kinetic events is chosen randomly with the probability proportional to its rate. The event will be chosen by j tot jR R u R 1 (8 2) where u is a random number between 0 to 1. Figure 81 shows a schematic illustration of the c hosen process. One of processes will be chosen from the distribution of rates proportional to the rate constant. The chosen event will be performed and the state i will move from another state k. The amount of time spent during this process is now calculat ed as:
156 totR u t log (8 3) where u is the random number between 0 to 1. Figure 81. Schematic of the procedure for choosing the reaction pathway to advance the system to the next state in the standard KMC algorithm. The size of box is proportional to the rate constant for each event. Figure 82 Flow chart illustrating the basic steps in the KMC algorithm.
157 Therefore, each step includes the probability of all events, even though that do not occurred. Figure 82 shows the flow chart for the KMC procedure. 8.2.2 Adaptive K inetic Mont e Carlo ( AKMC ) Because unidentified reaction pathways are almost invariably missing from the list of rates, KMC simulations cannot reproduce the real dynamics For the diffusion process, atomic movements can be different from those we might give by intuiti on alone. Although we can get more insight into the mechanism this approach also cannot guarantee to find all the reactive events that could occur during the evolution of the system. This is the reason AKMC method d eveloped. The AKMC method can reach long time scales while maintaining the accuracy of MD simulation Henkelman and Jonsson  have proposed a variation on the KMC method, in which the list of rates can be obtained at each state. The main idea to this approach is to search for saddle points that are connected to the current state of the syst em. The dimer method  is typically used for searching saddle points. . In the dimer algorithm, atoms climb uphill along the lowest eigenvector of the Hessian matrix until they reach the saddle point at the top. Figure 8 3 shows (a) a schematic view and (b) flow chart of the dimer searching met hod. The dimer method uses effective force s along the lowest curvature direction in order to find the saddle point. Thus, it is computationally efficient since the first derivatives of the potential is only required.  I f all saddle points are searched, the list of rates can be obtained from the activation energy for each saddle point Although the dimer method does not guarantee to find all saddles, most of the low lying barriers can be found by a large number of randomly initiated searches. 
158 Figure 8 3. (a) The schematic view of dimer saddle point searching mechanism.  (b) Flow chart of dimer method. 8.3 Dimer M ethod with T wo D imensional P otential I have impl emented t he dimer method into a new code. A two dimensional test problem case  has been analyzed in order to test the new developed code. The dimer method successfully identifies three different saddle points. Figure 8 4 shows the comparison between our test using the new code and literature resul ts  The comparison shows that our code suc cessfully reproduces the proper dimer behavior with this implemented two dimensional potential. (a) (b)
159 Figure 8 4 The comparison on searching behavior of dimer method between new developed code(a) and literature(b)  8. 4 Diffusion P rocess on (111) A luminum(Al) S urface KMC method is now also implement ed into the new code. In order to test the code, we simulate the diffusion process of an Al atom adsorbed on (111) Al surface. We used the simple LennardJones (LJ) potential for simplicity of the model problem. The standard from for LJ potential is ] ) / ( ) / [( 4 ) (6 12r r r VLJ LJ LJ LJ (8 4 ) w here, LJ =0.392 eV and LJ =2.64 are obtained from the literature.  A six layered Al structure which has a (111) planar surface is prepared and single Al atom is positioned on the top of t he surface. Figure 85 shows a schematic view of init ial Al structure with (111) surface. The color of the top Al atom is changed for only clari ty The (a) (b)
160 line s represent the border of the PBC cell box; therefore, atoms are repeated in x and y direction and only z direction remains as surface state with vacuum space. We performed AKMC simulation with the dimer method. Because the FCC structure has ABCABC staking sequence along the (111) planar direction, two adsorption s ite s, A and B are available, see figure 8 6. Thus, we can expect four different diffusion be havior s from two states occupying two s it es; A to A, A to B, B to A, B to B. All four processes are obs e r v ed in our AKMC study. Figure 87 shows four possible process obtained from our study. The study shows that the s t at e A has a 0.002 eV lower energy than that of st at e B. However, only one saddle point with a 0.102 eV higher from the s t at e A, is observed from all different processes Thus, two different activation proces ses are possible with activation energies, 0. 102 eV from s t at e A and 0. 100 eV from s t a t e B. We consider the energy at state A as a ground energy. Therefore the energy at state A is zero With the dimer method, the adsorbed atom moved to the transition state along the lowest curvature direction. The variation of substrate structure is the r esult of shifting of atoms along the PBC. The di mer searching steps, identifying saddle point s using dimer method, are repeated in order to find out almost all possible transition states. For this study, 12 dimer searching steps are used for each transitio n event. After searching 12 events of transition, the total rate for a transition event is calculated. Using the random number, one transition event among 12 different trajectories is chosen and accepted. With a different random number, time for the transi tion event is calculated. In this case, the calculated time for the first transition event is 2.11 ps. Figure 87 represents the time
161 evolution of kinetic process until the AKMC simulation reaches to 0. 1 5 ns. 130 transition events occur for this run. Figu re 8 5. The schematic view of six layered intial structure on (111) Al surface. The color of top Al atom is changed to brown for the clarification. The blue line represents the PBC cell box; therefore, atoms are repeated using the PBC in x and y direction and only z direction remains as surface state with vacuum space. Figure 86 Two possible adsorption sites for an Al adatom.
162 Figure 87 Four different diffusion processes of an Al adatom on (111) Al surface; (a) state A to state A (b) state A to state B (c) state B to state A (d) state B to state B. (b) (a) (c) (d)
163 Figure 88 The evolution of time for AKMC simulation. 75 transition events occur until it reaches the 0. 1 5 ns. 8. 5 Conclusion The AKMC code was developed and tested for a simple case. The AKMC code does not require the list of rates for the system, because it identifies for the transition state s using the dimer method. The implemented dimer method was test ed with a two dimensional potential case. The dimer method successfully identifies possible transition states in the two dimension potential. The AKMC code was tested with an Al adatom on the (111) Al surface. The AKMC code successfully searched the transition state and performed the trajectories of kinetic events. This work is an initial stage for studying the catalytic behavi or on ZnO surface. Further work will combine the AKMC method with other complex potentials, so as to be able to describe other material system s.
164 CHAPTER 9 GENERAL CONCLUSIONS Different level s of atomistic simulation methods have been successfully used to understand ferroelectric dom ain walls, organic gas surfaces and organic semiconductor interfaces. Also, an a daptive kinetic Monte Carlo ( AKMC ) code using dimer transition searching mechanism has been successfully developed in order to study catalytic activities on oxide surfaces. 9 .1 Domain W alls In this thesis, the structure and energetics of 180 domain walls of LiNbO3 were analyzed using electronic and atomic level simulations T he energetics of two structurally distinct doma in walls showed that Y walls parallel to c glide plane, were more favorable than X walls, perpendicular to c glide plane. The comparison o f activation energies of the two domain walls show ed that the Y wall is energetically less mobile than the X wall. Th is energetics show ed good agreement with the literature. In addition to the energy comparisons, the structural distinctions of two domain wall s were addressed by evaluating atomic displacements near domain walls. Because of the different structural orientations, one Y wa ll and two different X walls were available. Both Y wall and X walls cause d the atomic displacements along normal and transverse direction as well as longitudinal direction s. These displacements resulted in additional polarization components normal to the domain wall ( Pn) and transverse to the domain wall ( Pt) along with the primary ferroelectric polarization ( Pz). Thus, additional variations of polarization components are studied along with a kink type variation of Pz. For the Y wall, both o f Pn and Pt were observed near the domain wall center. The maximum
165 magnitude of Pn and Pt was calculated to be ~1.66 C/cm2 and 0.60 2 at the 2.58 away from the Y wall center. For X walls, only normal polarization components were observed in additio n to the primary ferroelectric polarization, because transverse components are cancelled by equal and opposite displacements of ions near X walls The maximum magnitude of Pn was observed to be 0.61 C/cm2 and 0.71 C/cm2 for two X walls at 1.49 away from the X walls. The presence of Pn and Pt along with Pz required a simultaneous fitting to calculate the domain wall width. However, none of previously defi n ed models were able to fit both components simultaneously. Thus, fits of Pz to single hyperbolic func tion yield ed sharp domain wall width ( ~ 2 .1 ) for both domain walls. A s imilar effect has recently observed in recent theoretical calculation in PbTiO3. D efects play an important role for the properties of bulk ferroelectric materials. The structure and energetics of intrinsic   and extrinsic defect  clusters in LiNbO3 had been studied by density functional theory (DFT). The pin ning of domain wall motion was observed by DFT study on oxygen vacancy domain wall interactions in PbTiO3.  However, no study has been made for defects domain wall interactions for LiNbO3. Thus, we perform ed a simulation study o f defects / domain wall interactions using DFT and empirical metho ds. Our study showed the most of defects / domain wall interactions lowered the defect formation energy significantly. Particular defect structures Li Frenkel and NbLi +4 VLi were considered. For the Li Frenkel pair, the p arallel configuration was energetically more favorable than antiparallel configuration. The interaction between a Li Frenkel pair and domain walls lowered the formation energy by almost onethird for both a Y wall and two X walls For NbLi +4 VLi however,
166 different energetics between Y wall and two X walls were observed. The formation energy of the parallel configurations wa s lowered by 0.24 eV/defect near a Y wall ; thereby, the parallel configuration was energetically favorable within entire regions across Y wall. However, possible rearrangement of defect pair between parallel and antiparallel configurations was observed near X walls. Also, X walls show ed higher defects domain wall interactions lowering defect formation energy of NbLi +4 VLi by 0.32 eV/defect than Y wall s Therefore, the NbLi +4 VLi defect cluster shows relative strong interaction with X walls. When we consider ed the similarity in DFE between LiNbO3 and LiTaO3, we could expect that TaLi +4 VLi defects will also interact stronger with an X wall than a Y wall. This is the major reason for switching of stable domain wall direction for LiTaO3. As the number of TaLi +4 VLi clusters increase, X wall could be energetically favorable in LiTaO3, which has only a small difference in domain wall energies The interactions of Y walls with extrinsic Mg dopants were also studied. Several different reaction mechanisms were considered for the study. Consistent with intrinsic defect clusters an Mg dopant pair also has a lower formation energy in the vicinity of Y wall. However, the interactions were relatively weaker than for intrinsic defect pairs. Since defect domain walls interactions cause d the pinning effect of domain wall motion, weaker defect domain wall interactions could improve the ferroelectric properties of LiNbO3. Therefore, our study successfully supported better optical resistance and relatively lower cohesive and threshold field of Mg doped LiNbO3 obs erved from experiments. Similar studies were performed on LiTaO3. Because an empirical potential for LiTaO3 system was not available, we developed an empirical potential for LiTaO3. The
167 new potential successfully captures the effect of temperature on ferroelectricity in particular the ferroelectric to paraelectric transition as well as structural parameter and elastic properties of LiTaO3. The energetics of domain walls in LiTaO3 show ed the Y wall wa s still energetically favorable for LiTaO3. However, the difference in energy between Y wall and X walls at the minimum energy state is small (3 mJ/m2 for DFT and 10 mJ/m2 for empirical) compared to LiNbO3 (40 mJ/m2 for DFT and 30 mJ/m2 for empirical). In addition, the energetics on defects domai n wall interacti ons demonstrated that defects are more strongly interacting with X walls than Y wall. Therefore, the shift of preferred domain wall direction from Y wall to X wall in Li deficient congruent LiTaO3 could be explained by defects / X wall interactions, which ov ercome the energy difference between Y wall and X wall. 9 .2 Organic Gas I nterface C ollision induced chemical reactions at organic gas interface were also investigated. U sing DFT MD methods, detailed processes for four different reaction mechanisms were r evealed. The results provide insights regarding reactions b etween organic molecules and predict how surface polymerization by ionassisted deposition (SPIAD) influences the polymerization and crosslinking process of polythiophene. The mechanism of surface polymerization of terthiophene oligomers was simulated using thiophene incident molecules. The DFT MD simulations predicted incident thiophenes having 25 eV kinetic energy could pos sibly yield surface reactions with a low probability (less than 30%). The i ncrease in the kinetic energy of incident thiophene to 35 eV could break the C C bonds of polythiophene and le a d ed to different surface reactions for polythiophenes. Further increases in the kinetic energy caused the damage to polythiophene and formed debr is from the collision. The mechanism for
168 polymerization and crosslinking of polythiophene was understood as a two step processes; breaking C S bonds of polythiophenes and creating stable C C bonds. Based on this mechanism four different reaction processes were observed; p olymerization, crosslinking, branched polymer and active state with brok en C S bond. Each reaction showed different probability. In a similar manner, bonding chemistry at or ganic semiconductor interface was also investigated. The incident C2H+ molec ule with 100 eV kinetic energy was deposited into the PbS nanocrystalline embedded in polythiophenes. The study on charge density distribution across the interface between thiophene and PbS quantum dot showed possible mechanism for creation of bonds. After the collision, a C S bond from dithiophene were broken and created a new bond with incident particle. This new bond produced the electronic connection between PbS nanocrystalline and polythiophene. Therefore breaking C S bonds of polythiophenes wa s a key mechanism for reaction of polythiophenes. 9.3 Metal gas I nterface Surface diffusion of Al adatom on Al (111) surface was studied with a new ly developed AKMC code. Two different sites were available for Al adatom on (111) surface. The study with LennardJones potential obtained from literature showed the energy difference between two different sites was 0.002 eV. AKMC simulations were performed to understand diffusion process of an Al adatom on Al (111) surface. Our study showed that only one trans ition state was available for surface diffusion of an Al adatom. So, the study yielded two different activation energies and four different diffusion processes. As we denoted two different positions of an adatom as A and B, the activation energy for leavin g two sites was 0.100 eV for site A and 0.102 eV for site B. Four different diffusion
169 processes were A to A, A to B, B to A and B to B. This study will be expanded to understand catalytic activity on ZnO2 surface.
170 LIST OF REFERENCES  C. R. Barrett, W. D. Nix, and A. S. Tetelman, The Principles of Engineering Materials New Jersey: Prentice Hall, 1973.  R. M. Martin, Electronic Structure: Basic Theory and Practical Methods New York: Cambridge University Press, 2004.  E. Schr odinger, "An undulatory theory of the mechanics of atoms and molecules," Physical Review, vol. 28, pp. 10491070, 1926.  M. Born and R. Oppenheimer, "Quantum theory of molecules," Annalen Der Physik, vol. 84, pp. 04570484, 1927.  J. C. Slater, "Th e theory of complex spectra," Physical Review, vol. 34, pp. 1293 1322, 1929.  D. S. Scholl and J. A. Steckel, Density Functional Theory: A Practical Introduction Hoboken, New Jersey: John Wiley & Sons, Inc., 2009.  P. Hohenberg and W. Kohn, "Inhom ogeneous Electron Gas," Phys. Rev., vol. 136, pp. B864B871, 1964.  W. Kohn and L. J. Sham, "Self Consistent Equations Including Exchange And Correlation Effects," Phys. Rev., vol. 140, pp. A1133A1138, 1965.  C. J. Cramer, "Essentials of Computati onal Chemistry: Theories and Models," Second Edition ed: John Wiley & Sons, Ltd, 2004.  J. P. Perdew, "Accurate Density Functional for the Energy: Real Space Cutoff of the Gradient Expansion for the Exchange Hole," Phys. Rev. Lett., vol. 55, pp. 16651668, 1985.  N. W. Ashcroft and N. D. Mermin, Solid State Physics New York: Holt, Rinehart and Winston, 1976.  H. J. Monkhorst and J. D. Pack, "Special points for Brillouinzone integrations," Phys. Rev. B, vol. 13, pp. 51885192, 1976.  S. F. Boys, "Electronic wave functions .1. A general method of calculation for the stationary states of any molecular system," Proceedings of the Royal Society of London Series aMathematical and Physical Sciences, vol. 200, pp. 542554, 1950.  G. Kress e and J. Furthmuller, "Efficiency of abinitio total energy calculations for metals and semiconductors using a planewave basis set," Comp. Mater. Sci., vol. 6, pp. 1550, 1996.
171  G. Kresse and J. Furthmuller, "Efficient iterative schemes for ab initio total energy calculations using a planewave basis set," Phys. Rev. B, vol. 54, pp. 1116911186, 1996.  D. Vanderbilt, "Soft self consistent pseudopotentials in a generalized eigenvalue formalism," Phys. Rev. B, vol. 41, pp. 78927895, 1990.  P. E. Blochl, "PROJECTOR AUGMENTED WAVE METHOD," Phys. Rev. B, vol. 50, pp. 1795317979, 1994.  G. Kresse and J. Joubert, "From Ultrasoft Pseudopotentials to the Projector AugmentedWave Method," Phys. Rev. B, vol. 59, pp. 17581775, 1999.  D. Lee, H. X. Xu, V. Dierolf, V. Gopalan, and S. R. Phillpot, "Structure and energetics of ferroelectric domain walls in LiNbO3 from atomic level simulation," Phys. Rev. B, p. to be Submitted, 2010.  J. P. Perdew and W. Yue, "Accurate And Simple Density Func tional For The Electronic Exchange Energy Generalized Gradient Approximation," Phys. Rev. B, vol. 33, pp. 88008802, 1986.  Q. K. Li, B. Wang, C. H. Woo, H. Wang, and R. Wang, "Firstprinciples study on the formation energies of intrinsic defects in LiNbO3," J. Phys. Chem. Solids, vol. 68, pp. 13361340, 2007.  M. P. Teter, M. C. Payne, and D. C. Allan, "Solution Of Schrodinger Equation For Large Systems," Phys. Rev. B, vol. 40, pp. 1225512263, 1989.  P. Pulay, "Convergence Acceleration Of Iterative Sequences The Case Of Scf Iteration," Chem. Phys. Lett., vol. 73, pp. 393398, 1980.  H. X. Xu, D. Lee, J. He, S. B. Sinnott, V. Gopalan, V. Dierolf, and S. R. Phillpot, "Stability of intrinsic defects and defect clusters in LiNbO3 from density functional theory calculations," Phys. Rev. B, vol. 78, p. 174103, Nov 2008.  N. Troullier and J. L. Martins, "Efficient Pseudopotentials for PlaneWave Calculations," Physical Review B, vol. 43, pp. 19932006, Jan 15 1991.  J. M. Soler, E Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez Portal, "The SIESTA method for ab initio order N materials simulation," Journal of Physics: Condensed Matter, vol. 14, pp. 27452779, 2002.  G. D. Purvis and R. J. Bartlett, "A Full Coupled Cluster Singles and Doubles Model the Inclusion of Disconnected Triples," Journal of Chemical Physics, vol. 76, pp. 19101918, 1982.
172  T. H. Dunning, "GaussianBasis Sets for Use in Correlated Molecular Calculations .1. The Atoms Boron t hrough Neon and Hydrogen," Journal of Chemical Physics, vol. 90, pp. 10071023, Jan 15 1989.  J. A. Pople, M. Headgordon, and K. Raghavachari, "Quadratic ConfigurationInteraction a General Technique for Determining Electron Correlation Energies," Journal of Chemical Physics, vol. 87, pp. 59685975, Nov 15 1987.  R. S. Mulliken, "Electronic Population Analysis On LcaoMo Molecular Wave Functions .1," J. Chem. Phys., vol. 23, pp. 18331840, 1955.  J. P. Foster and F. Weinhold, "Natural Hybri d Orbitals," J. Am. Chem. Soc., vol. 102, pp. 72117218, 1980.  P. O. Lowdin, "On The NonOrthogonality Problem Connected With The Use Of Atomic Wave Functions In The Theory Of Molecules And Crystals," J. Chem. Phys., vol. 18, pp. 365375, 1950.  R. F. W. Bader and T. T. Nguyendang, "Quantum Theory Of Atoms In Molecules Dalton Revisited," Adv. Quantum Chem., vol. 14, pp. 63124, 1981.  G. Henkelman, A. Arnaldsson, and H. Jonsson, "A fast and robust algorithm for Bader decomposition of charg e density," Comp. Mater. Sci., vol. 36, pp. 354360, Jun 2006.  E. Sanville, S. D. Kenny, R. Smith, and G. Henkelman, "Improved gridbased algorithm for Bader charge allocation," J. Comput. Chem., vol. 28, pp. 899908, 2007.  W. Tang, E. Sanville and G. Henkelman, "A gridbased Bader analysis algorithm without lattice bias," J. Phys. Condens. Matter, vol. 21, p. 084204, 2009.  A. J. Cohen, P. Mori Sanchez, and W. T. Yang, "Insights into current limitations of density functional theory," Scie nce, vol. 321, pp. 792794, 2008.  Y. K. Zhang and W. T. Yang, "A challenge for density functionals: Self interaction error increases for systems with a noninteger number of electrons," J. Chem. Phys., vol. 109, pp. 26042608, 1998.  P. Mori Sanc hez, Q. Wu, and W. T. Yang, "Accurate polymer polarizabilities with exact exchange density functional theory," J. Chem. Phys., vol. 119, pp. 1100111004, 2003.
173  C. Toher, A. Filippetti, S. Sanvito, and K. Burke, "Self interaction errors in density fun ctional calculations of electronic transport," Phys. Rev. Lett., vol. 95, 2005.  S. H. Ke, H. U. Baranger, and W. T. Yang, "Role of the exchangecorrelation potential in ab initio electron transport calculations," J. Chem. Phys., vol. 126, 2007.  S. Grimme, J. Antony, T. Schwabe, and C. Muck Lichtenfeld, "Density functional theory with dispersion corrections for supramolecular structures, aggregates, and complexes of (bio)organic molecules," Organic & Biomolecular Chemistry, vol. 5, pp. 741758, 2007.  R. C. Albers, N. E. Christensen, and A. Svane, "HubbardU bandstructure methods," J. Phys. Condens. Matter, vol. 21, 2009.  A. J. Cohen, P. Mori Sanchez, and W. T. Yang, "Development of exchangecorrelation functionals with minimal many el ectron self interaction error," J. Chem. Phys., vol. 126, 2007.  M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids First ed. New York: Oxford University Press Inc., 1987.  D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r( 1) summation," J. Chem. Phys., vol. 110, pp. 82548282, May 1999.  R. A. Buckingham, "The classical equation of state of gaseous helium, neon and argon," Procee dings of the Royal Society of London Series aMathematical and Physical Sciences, vol. 168, pp. 264283, 1938.  J. E. Jones, "On the determination of molecular fields II From the equation of state of a gas," Proceedings of the Royal Society of London Series a Containing Papers of a Mathematical and Physical Character, vol. 106, pp. 463477, 1924.  P. M. Morse, "Diatomic molecules according to the wave mechanics. II. Vibrational levels," Physical Review, vol. 34, pp. 5764, 1929.  R. S. Chhikara and L. Folks, The Inverse Gaussian Distribution: Theory, Methodology, and Applications. Statistics, textbooks and monographs vol. 95. New York: M. DEKKER, 1989.  K. T. Tang and J. P. Toennies, "An Improved SimpleModel For The Vanderwaals Potenti al Based On Universal Damping Functions For The Dispersion Coefficients," J. Chem. Phys., vol. 80, pp. 37263741, 1984.
174  Y. P. Varshni, "Comparative Study Of Potential Energy Functions For Diatomic Molecules," Reviews of Modern Physics, vol. 29, pp. 664 682, 1957.  R. A. Jackson and M. E. G. Valerio, "A new interatomic potential for the ferroelectric and paraelectric phases of LiNbO3," J. Phys. Condens. Matter, vol. 17, pp. 837843, Feb 2005.  L. Verlet, "Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of LennardJones Molecules," Physical Review, vol. 159, pp. 98103, 1967.  S. Toxvaerd, "Molecular Dynamics At Constant Temperature And Pressure," Physical Revi ew E, vol. 47, pp. 343350, 1993.  D. Beeman, "Some Multistep Methods For Use In Molecular Dynamics Calculations," Journal of Computational Physics, vol. 20, pp. 130139, 1976.  C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations Englewood Cliffs, NJ: PrenticeHall, 1971.  M. E. Tuckerman and G. J. Martyna, "Understanding modern molecular dynamics: Techniques and applications," J. Phys. Chem. B, vol. 104, pp. 159178, 2000.  M. P. Allen and D. J. Tildesley, C omputer Simulation of Liquids Oxford: Oxford Science, 2004.  D. G. Luenberger and Y. Ye, Linear and Nonlinear Programming, Third ed. New York: Springer Science+Business media, LLC, 2008.  M. Parrinello and A. Rahman, "Crystal Structure and Pair Potentials: A Molecular Dynamics Study," Phys. Rev. Lett., vol. 45, pp. 11961199, 1980.  M. Parrinello and A. Rahman, "Strain fluctuations and elastic constants," J. Chem. Phys., vol. 76, pp. 26622666, 1982.  H. C. Andersen, "Molecular dynamics simulations at constant pressuure and/or temperature," J. Chem. Phys., vol. 72, pp. 23842393, 1980.  S. Nos, "A unified formulation of the constant temperature molecular dynamics methods," J. Chem. Phys., vol. 81 pp. 511519, 1984.  W. G. Hoov er, "Canonical Dynamics Equilibrium PhaseSpace Distributions," Physical Review A vol. 31, pp. 16951697, 1985.
175  H. J. C. Berendsen, J. P. M. Postma, W. F. v. Gunsteren, A. DiNola, and J. R. Haak, "Molecular dynamics with coupling to an external bat h," J. Chem. Phys., vol. 81 pp. 36843690, 1984.  S. A. Adelman and J. D. Doll, "Generalized Langevin Equation Approach for Atom/Solid Surface Scattering: Collinear Atom/Harmonic Chain Model," J. Chem. Phys., vol. 61, pp. 42424245, 1974.  S. A. Adelman and J. D. Doll, "Generalized Langevin Equation Approach for Atom Solid Surface Scattering General Formulation for Classical Scattering Off Harmonic Solids," J. Chem. Phys., vol. 64, pp. 2375 2388, 1976.  K. K. Wong, Properties of Lithium Niobate. London: INSPEC, The Institution of Electrical Engineers, United Kingdom, 2002.  T. Volk and M. Whlecke, Lithium Niobate: Defects, Photorefraction and Ferroelectric Switching Berlin: Springer, 2008.  P. F. Bordui, R. G. Norwood, D. H. Jund t, and M. M. Fejer, "Preparation And Characterization Of Off Congruent Lithium Niobate Crystals," J. Appl. Phys., vol. 71, pp. 875879, 1992.  Y. H. Han and D. M. Smyth, "Nonstoichiometry And Defects In Linbo3," American Ceramic Society Bulletin, vol. 62, pp. 852852, 1983.  H. M. Obryan, P. K. Gallagher, and C. D. Brandle, "Congruent Composition And Li Rich PhaseBoundary Of Linbo3," Journal of the American Ceramic Society, vol. 68, pp. 493496, 1985.  N. Iyi, K. Kitamura, F. Izumi, J. K. Ya mamoto, T. Hayashi, H. Asano, and S. Kimura, "ComparativeStudy Of Defect Structures In Lithium Niobate With Different Compositions," Journal of Solid State Chemistry, vol. 101, pp. 340352, 1992.  W. Bollmann, "Stoichiometry And Point Defects In Lith ium Niobate Crystals," Crystal Research and Technology, vol. 18, pp. 11471149, 1983.  O. F. Schirmer, O. Thiemann, and M. Wohlecke, "Defects in linbo3 .1. Experimental aspects," J. Phys. Chem. Solids, vol. 52, pp. 185200, 1991.  H. Donnerberg, S. M. Tomlinson, C. R. A. Catlow, and O. F. Schirmer, "Computer Simulation Studies Of Intrinsic Defects In Linbo3 Crystals," Phys. Rev. B, vol. 40, pp. 1190911916, Dec 1989.
176  H. X. Xu, D. Lee, S. B. Sinnott, V. Dierolf, V. Gopalan, and S. R. Phillpot "Structure and diffusion of intrinsic defect complexes in LiNbO3 from density functional theory calculations," J. Phys. Condens. Matter, vol. 22, 2010.  R. M. Araujo, K. Lengyel, R. A. Jackson, L. Kovacs, and M. E. G. Valerio, "A computational study of intrinsic and extrinsic defects in LiNbO3," J. Phys. Condens. Matter, vol. 19, 2007.  V. Gopalan, V. Dierolf, and D. A. Scrymgeour, "Defect domain wall interactions in trigonal ferroelectrics," Ann. Rev. Mater. Res., vol. 37, pp. 449489, 2007. [ 81] P. Paufler, "Complete online set of International tables for crystallography," Acta Crystallogr., Sect. A: Found. Crystallogr., vol. 63(6), p. 483, 2007.  H. Boysen and F. Altorfer, "A Neutron Powder Investigation Of The HighTemperature Structure And PhaseTransition In Linbo3," Acta Crystallogr., Sect. B: Struct. Sci, vol. 50, pp. 405414, 1994.  K. Parlinski and Z. Q. Li, "Ab initio calculations of phonons in LiNbO3," Phys. Rev. B, vol. 61, pp. 272278, 2000.  A. V. Postnikov, V. Caciu c, and G. Borstel, "Structure optimization and frozen phonons in LiNbO3," J. Phys. Chem. Solids, vol. 61, pp. 295299, 2000.  M. Veithen and P. Ghosez, "First principles study of the dielectric and dynamical properties of lithium niobate," Phys. Rev. B, vol. 65, p. 214302, 2002.  I. U. o. Crystallography, International tables for crystallography New York: Springer, 2006.  S. C. Abrahams, J. M. Reddy, and Bernstei.Jl, "Ferroelectric lithium niobate .3. Single crystal x ray diffraction study at 24 degrees c," J. Phys. Chem. Solids, vol. 27, pp. 997&, 1966.  H. D. Megaw, "A Note On Structure Of Lithium Niobate Linbo3," Acta Crystallogr., sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., vol. A 24, pp. 583&, 1968.  S. C. Abraha ms, Levinste.Hj, and J. M. Reddy, "Ferroelectric lithium niobate .5. Polycrystal x ray diffraction study between 24 degrees and 1200 degrees c," J. Phys. Chem. Solids, vol. 27, pp. 1019&, 1966.  A. M. Glazer, "Simple Ways Of Determining Perovskite St ructures," Acta Crystallogr. Sect. A, vol. 31, pp. 756762, 1975.  H. D. Megaw, Crystal structures a working approach. Philadelphia: W. B. Saunders, 1973.
177  H. Landolt and R. Brnstein, Numerical Data and Functional Relations in Science and Techn ology vol. 15a. Berlin: Springer, 1981.  S. R. Phillpot and V. Gopalan, "Coupled displacive and order disorder dynamics in LiNbO3 by molecular dynamics simulation," Appl. Phys. Lett., vol. 84, pp. 19161918, Mar 2004.  O. Auciello, J. F. Scott, a nd R. Ramesh, "The physics of ferroelectric memories," Physics Today, vol. 51, pp. 2227, Jul 1998.  A. J. Moulson and J. M. Herbert, Electroceramics : materials, properties, applications Second ed. New York: Wiley, J., 2003.  D. A. Scrymgeour, V. Gopalan, A. Itagi, A. Saxena, and P. J. Swart, "Phenomenological theory of a single domain wall in uniaxial trigonal ferroelectrics: Lithium niobate and lithium tantalate," Phys. Rev. B, vol. 71, p. 184110, May 2005.  H. S. Nalwa, Handbook of Organic Conductive Molecules and Polymers Chichester: Wiley, 1997.  T. A. Skotheim, R. L. Elsenbaumer, and J. R. Reynolds, Handbook of Conducting Polymers New York: M. Dekker, 1998.  F. Jonas and G. Heywang, "Technical applications for conductive polymers," Electrochim. Acta, vol. 39, pp. 13451347, 1994.  S. Ghosh and O. Ingan, "Conducting Polymer Hydrogels as 3D Electrodes: Applications for Supercapacitors," Adv. Mater., vol. 11, pp. 12141218, 1999.  M. Gerard, A. Chaubey, and B. D. Malhotra, "Application of conducting polymers to biosensors," Biosens. Bioelectron., vol. 17, pp. 345359, 2002.  C. Pratt, "Conducting Polymers," http://homepage.ntlworld.com/colin.pratt/cpoly.pdf 1996.  D. Fichou, "Handbook of Oligoand Polythiophenes," Wiley VCH: Weinheim, 1999.  A. J. Heeger, "Semiconducting and metallic polymers: The fourth generation of polymeric materials," J. Phys. Chem. B, vol. 105, pp. 84758491, 2001.  S. Tepavcevic, Y. Choi, and L. Hanley, "Surface polymerization by ionassisted deposition for polythiophene film growth," J. Am. Chem. Soc., vol. 125, pp. 23962397, 2003.
178  G. Tourillon, "Conjugate polymers: processing and applications, in Handbook of Conducting Polymers Third ed. vol. 1, T. A. Skotheim and J. R. Reynolds, Eds.: CRC Press, Taylor & Francis Group, LLC, 2007, p. 293.  A. J. Heeger, D. Moses, and M. Sinclair, "Nonlinear Excitations And Nonlinear Phenomena In Conduct ive Polymers," ACS Symp. Ser., vol. 346, pp. 372380, 1987.  F. Vanbolhuis, H. Wynberg, E. E. Havinga, E. W. Meijer, and E. G. J. Staring, "The X Ray Structure And Mndo Calculations Of AlphaTerthienyl A Model For Polythiophenes," Synth. Met., vol. 30, pp. 381389, 1989.  G. J. Ashwell, Molecular Electronics Taunton, Somerset, England: Research Studies Press, 1992.  L Alccer Conducting Polymers; Special Applications : Proceedings of the Workshop Held at Sintra Portugal, July 2831, 1986. Dordrecht: D. Reidel, 1987.  G. A. Inzelt, Conducting Polymers New York: Springer, 2008.  B. J. Schwartz, "Conjugated polymers as molecular materials : How chain conformation and film morphology influence energy transfer and interchain i nteractions," Annu. Rev. Phys. Chem., vol. 54, pp. 141172, 2003.  J. R. Sheats, "Manufacturing and commercialization issues in organic electronics," J. Mater. Res., vol. 19, pp. 1974 1989, 2004.  Y. S. Choi, S. Tepavcevic, Z. Xu, and L. Hanley "Optical and chemical properties of polythiophene films produced via surface polymerization by ionassisted deposition," Chem. Mater., vol. 16, pp. 19241931, 2004.  S. Tepavcevic, Y. Choi, and L. Hanley, "Growth of novel polythiophene and polyphenyl films via surface polymerization by ionassisted deposition," Langmuir, vol. 20, pp. 87548761, 2004.  Y. Choi, A. Zachary, S. Tepavcevic, C. P. Wu, and L. Hanley, "Polyatomicity and kinetic energy effects on surface polymerization by ionassisted deposition," Int. J. Mass spectrom., vol. 241, pp. 139147, 2005.  S. Tepavcevic, A. T. Wroble, M. Bissen, D. J. Wallace, Y. Choi, and L. Hanley, "Photoemission studies of polythiophene and polyphenyl films produced via surface polymerization by ionassisted deposition," J. Phys. Chem. B, vol. 109, pp. 71347140, 2005.
179  J. W. Shen, C. Evans, N. Wade, and R. G. Cooks, "Ion ion collisions leading to formation of C C bonds at surfaces: An interfacial Kolbe reaction," J. Am. Chem. Soc., vol. 121, pp. 9762 9763, 1999.  N. Wade, C. Evans, S. C. Jo, and R. G. Cooks, "Silylation of an OH terminated self assembled monolayer surface through low energy collisions of ions: a novel route to synthesis and patterning of surfaces," J. Mass Spectrom., vol. 37, pp. 591602, 2002.  H. Usui, "Polymeric film deposition by ionizationassisted method for optical and optoelectronic applications," Thin Solid Films, vol. 365, pp. 2229, 2000.  J. Y. Kim, E. S. Kim, and J. H. Choi, "Poly[2 (N carbazolyl) 5 (2 ethylhexyloxy) 1,4 phenylenevinylene]/tris(8hydroxyquinoline) aluminum heterojunction electroluminescent devices produced by cluster beam deposition methods," J. Appl. Phys., vol. 91, pp. 19441951, 2002.  L. Hanley, Y. Choi, and S. Tepavcevic Materials Research Society Symposium Proceeding, pp. 804, JJ5.3.1JJ5.3.6, 2004.  I. Jang, B. Ni, and S. B. Sinnott, "Study of angular influence of C3H5+ ion deposition on polystyrene surfaces using molecular dynamics simulations," Journal of Vacuum Science & Technology aVacuum Surfaces and Films, vol. 20, pp. 564568, 2002.  D. G. Lim, B. S. Jang, S. I. Moon, C. Y. Won, and J. Yi, "Characteristics of LiNbO3 memory capacitors fabricated using a low thermal budget process," Solid State Electronics, vol. 45, pp. 11591163, 2001.  J. Macario, P. Yao, R. Shireen, C. A. Schuetz, S. Y. Shi, and D. W. Prather, "Development of ElectroOptic Phase Modulator for 94 GHz Imaging System," J. Lightwave Technol., vol. 27, pp. 56985703, 2009.  A. Tehranchi and R. Kashyap, "Efficient wavelength converters with flattop responses based on counterpropagating cascaded SFG and DFG in low loss QPM LiNbO3 waveguides," Optics Express, vol. 17, pp. 1911319119, 2009.  A. M. Prokhorov and Y. S. Kuz'minov, Physics and Chemistry of Crystalline Lithium Niobate Taylor & Francis, 1990.  L. Arizmendi, "Photonic applications of lithium niobate crystals," Phys. Status Solidi A Appl. Res., vol. 201, pp. 253283, Jan 2004.
180  V. Gopalan, M. J. Kawas, M. C. Gupta, T. E. Schlesinger, and D. D. Stancil, "Integrated quasi phasematched secondharmonic generator and electrooptic scanner on LiTaO3 single crystals," IEEE Photonics Technol. Lett., vol. 8, pp. 17041706, Dec 1996.  V. Soluch and M. Lysakowska, "Surface acoustic waves on X cut LiNbO3," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 52, pp. 145147, Jan 2005.  J. Padilla, W. Zhong, and D. Vanderbilt, "First principles investigation of 180 degrees domain walls in BaTiO3," Phys. Rev. B, vol. 53, pp. R5969R5973, Mar 1996.  A. Goldman, The Handbook of Modern Ferromagnetic Materials 1st edition ed. New York: Springer, 1999.  B. Meyer and D. Vanderbilt, "Ab initio study of ferroelectric domain walls in PbTiO3," Phys. Rev. B, vol. 65, p. 104111, 2002.  K. C. Kao, Dielectric Phenomena in Solids 1st edition ed. New York: Academic Press, 2004.  D. Lee, R. K. Behera, P. Wu, H. X. Xu, S. B. Sinnott, S. R. Phillpot, L. Q. Chen, and V. Gopalan, "Mixed BlochNel Ising character of 180 ferroelectric domain walls," Phys. Rev. B, vol. 80, p. 060102(R), August 2009.  A. F. Devonshire, "Theory Of Barium Titanate .2," Philos. Mag., vol. 42, pp. 10651079, 1951.  V. A. Zhirnov, "A Contribution To The Theory Of Dom ain Walls In Ferroelectrics," Zh. Eksp. Teor. Fiz., vol. 35, p. 1175, 1958.  V. L. Ginzburg, "Some Remarks On Phase Transitions Of The 2nd Kind And The Microscopic Theory Of Ferroelectric Materials," Fiz. Tverd. Tela, vol. 2, p. 2031, 1960.  F. Jona and G. Shirane, Ferroelectric crystals Oxford: Pergamon Press, 1962.  L. N. Bulaevskii, "Thermodynamic Theory Of Domain Walls In Ferroelectric Materials With The Perovskite Structure," Fiz. Tverd. Tela (St. Petersburg), vol. 5, p. 3183, 1963.  Y. Ishibashi, "Computational Method Of ActivationEnergy Of Thick DomainWalls," J. Phys. Soc. Jpn., vol. 46, pp. 12541257, 1979.
181  Y. Ishibashi and I. Suzuki, "The Structure, The Formation And The ActivationEnergies Of A Domain Wall In Cryst als Undergoing The 1st Order Transition," J. Phys. Soc. Jpn., vol. 53, pp. 1093 1096, 1984.  S. Choudhury, Y. L. Li, N. Odagawa, A. Vasudevarao, L. Tian, P. Capek, V. Dierolf, A. N. Morozovska, E. A. Eliseev, S. Kalinin, Y. S. Cho, L. Q. Chen, and V. Gopalan, "The influence of 180 degrees ferroelectric domain wall width on the threshold field for wall motion," J. Appl. Phys., vol. 104, p. 084107, 2008.  M. Leslie and M. J. Gillan, "The Energy And Elastic Dipole Tensor Of Defects In Ionic Crystal s Calculated By The Supercell Method," Journal of Physics C Solid State Physics, vol. 18, pp. 973982, 1985.  M. F. Mott and M. J. Littleton, "Conduction in polar crystals. I. Electrolytic conduction in solid salts," Transactions of the Faraday Society, vol. 34, pp. 0485 0499, 1938.  C. G. Van de Walle and J. Neugebauer, "First principles calculations for defects and impurities: Applications to III nitrides," Journal of Applied Physics, vol. 95, pp. 38513879, Apr 15 2004.  K. Reuter and M Scheffler, "Composition, structure, and stability of RuO2(110) as a function of oxygen pressure," Physical Review B, vol. 65, pp. Jan 15 2002.  J. He, R. K. Behera, M. W. Finnis, X. Li, E. C. Dickey, S. R. Phillpot, and S. B. Sinnott, "Predictio n of high temperature point defect formation in TiO2 from combined ab initio and thermodynamic calculations," Acta Materialia, vol. 55, pp. 43254337, Aug 2007.  A. F. Kohan, G. Ceder, D. Morgan, and C. G. Van de Walle, "First principles study of nat ive point defects in ZnO," Physical Review B, vol. 61, pp. 1501915027, Jun 1 2000.  I. P. Zibrov, V. P. Filonenko, P. E. Werner, B. O. Marinder, and M. Sundberg, "A new highpressure modification of Nb2O5," Journal of Solid State Chemistry, vol. 141, pp. 205 211, Nov 15 1998.  J. B. Goodenough and A. Hamnett, Physics of Non Tetrahedrally Bonded Binary Compounds III Berlin: Springer Verlag, 1984.  S. Kim, V. Gopalan, K. Kitamura, and Y. Furukawa, "Domain reversal and nonstoichiometry in l ithium tantalate," J. Appl. Phys., vol. 90, pp. 29492963, 2001.  L. X. He and D. Vanderbilt, "First principles study of oxygenvacancy pinning of domain walls in PbTiO3," Phys. Rev. B, vol. 68, 2003.
182  H. Donnerberg, S. M. Tomlinson, C. R. A. Ca tlow, and O. F. Schirmer, "Computer Simulation Studies Of Extrinsic Defects In Linbo3 Crystals," Phys. Rev. B, vol. 44, pp. 48774883, 1991.  H. X. Xu, D. Lee, S. B. Sinnott, V. Gopalan, V. Dierolf, and S. R. Phillpot, "Structure and energetics of Er defects in LiNbO3 from first principles and thermodynamic calculations," Phys. Rev. B, vol. 80, 2009.  F. Abdi, M. Aillerie, P. Bourson, and M. D. Fontana, "Defect structure in Mg doped LiNbO3: Revisited study," J. Appl. Phys., vol. 106, 2009. [157 ] H. Z. Jin and J. Zhu, "Size effect and fatigue mechanism in ferroelectric thin films," J. Appl. Phys., vol. 92, pp. 45944598, Oct 2002.  B. S. Razbirin, "Curie Temperature Of The Ferroelectric Litao3 Edge Luminescence Of Crystals," Soviet PhysicsSolid State, vol. 6, pp. 254255, 1964.  S. H. Wemple, Didomeni.M, and I. Camlibel, "Relationship Between Linear And Quadratic ElectroOptic Coefficients In Linbo3 Litao3 And Other OxygenOctahedra Ferroelectrics Based On Direct Measurement Of Sponta neous Polarization," Appl. Phys. Lett., vol. 12, pp. 209&, 1968.  R. G. Batchko, V. Y. Shur, M. M. Fejer, and R. L. Byer, "Backswitch poling in lithium niobate for high fidelity domain patterning and efficient blue light generation," Appl. Phys. Let t., vol. 75, pp. 16731675, 1999.  K. T. Gahagan, V. Gopalan, J. M. Robinson, Q. Z. X. Jia, T. E. Mitchell, M. J. Kawas, T. E. Schlesinger, and D. D. Stancil, "Integrated electrooptic lens/scanner in a LiTaO3 single crystal," Applied Optics, vol. 38, pp. 11861190, 1999.  D. Psaltis and G. W. Burr, "Holographic data storage," Computer, vol. 31, pp. 52 +, 1998.  K. Kitamura, J. K. Yamamoto, N. Iyi, S. Kimura, and T. Hayashi, "Stoichiometric Linbo3 SingleCrystal Growth By Double Crucible C zochralski Method Using Automatic Powder Supply System," J. Cryst. Growth, vol. 116, pp. 327332, 1992.  L. Tian, V. Gopalan, and L. Galambos, "Domain reversal in stoichiometric LiTaO3 prepared by vapor transport equilibration," Appl. Phys. Lett., vo l. 85, pp. 44454447, 2004.
183  V. Gopalan, N. A. Sanford, J. A. Aust, K. Kitamura, Y. Furukawa, N. Hari Singh, M.Sc, and Ph.D, "Crystal growth, characterization, and domain studies in lithium niobate and lithium tantalate ferroelectrics," in Handbook of Advanced Electronic and Photonic Materials and Devices Burlington: Academic Press, 2001, pp. 57114.  V. Gopalan, T. E. Mitchell, Y. Furukawa, and K. Kitamura, "The role of nonstoichiometry in 180 degrees domain switching of LiNbO3 crystals," Ap pl. Phys. Lett., vol. 72, pp. 19811983, 1998.  K. Kitamura, Y. Furukawa, K. Niwa, V. Gopalan, and T. E. Mitchell, "Crystal growth and low coercive field 180 degrees domain switching characteristics of stoichiometric LiTaO3," Appl. Phys. Lett., vol. 73, pp. 30733075, 1998.  V. Y. Shur, "Kinetics of ferroelectric domains: Application of general approach to LiNbO3 and LiTaO3," Journal of Materials Science, vol. 41, pp. 199210, 2006.  S. C. Abrahams, E. Buehler, W. C. Hamilton, and S. J. La placa, "Ferroelectric lithium tantalate .3. Temperaturedependence of structure in ferroelectric phase and paraelectric structure at 940 degrees k," J. Phys. Chem. Solids, vol. 34, pp. 521532, 1973.  I. Takanaga and J. Kushibiki, "Elastic constants of multidomain LiTaO3 crystal," J. Appl. Phys., vol. 86, pp. 33423346, 1999.  A. M. Glass, "Dielectric Thermal And Pyroelectric Properties Of Ferroelectric Litao3," Physical Review, vol. 172, pp. 564&, 1968.  R. Peierls, "The size of a disloc ation," Proc. Phys. Soc, vol. 52, pp. 3437, Jan 1940.  R. T. Smith and F. S. Welsh, "Temperature Dependence Of Elastic, Piezoelectric, And Dielectric Constants Of Lithium Tantalate And Lithium Niobate," J. Appl. Phys., vol. 42, pp. 2219&, 1971. [1 74] V. Gopalan and M. C. Gupta, "Observation of internal field in LiTaO3 single crystals: Its origin and timetemperature dependence," Appl. Phys. Lett., vol. 68, pp. 888890, 1996.  W. D. Hsu, S. Tepavcevic, L. Hanley, and S. B. Sinnott, "Mechanisti c studies of surface polymerization by ionassisted deposition," Journal of Physical Chemistry C, vol. 111, pp. 41994208, 2007.  E. H. Sargent, "Solar Cells, Photodetectors, and Optical Sources from Infrared Colloidal Quantum Dots," Adv. Mater., vol 20, pp. 39583964, 2008.
184  A. L. Rogach, A. Eychmuller, S. G. Hickey, and S. V. Kershaw, "Infraredemitting colloidal nanocrystals: Synthesis, assembly, spectroscopy, and applications," Small, vol. 3, pp. 536557, Apr 2007.  A. A. R. Watt, D. Blake, J. H. Warner, E. A. Thomsen, E. L. Tavenner, H. RubinszteinDunlop, and P. Meredith, "Lead sulfide nanocrystal: conducting polymer solar cells," J. Phys. D Appl. Phys., vol. 38, pp. 20062012, Jun 2005.  D. H. Cui, J. Xu, T. Zhu, G. Paradee, S Ashok, and M. Gerhold, "Harvest of near infrared light in PbSe nanocrystal polymer hybrid photovoltaic cells," Appl. Phys. Lett., vol. 88, p. 3, May 2006.  K. P. Fritz, S. Guenes, J. Luther, S. Kumar, N. S. Saricifitci, and G. D. Scholes, "IV VI na nocrystal polymer solar cells," J. Photochem. Photobiol. A Chem., vol. 195, pp. 3946, Mar 2008.  S. J. Kim, W. J. Kim, A. N. Cartwright, and P. N. Prasad, "Carrier multiplication in a PbSe nanocrystal and P3HT/PCBM tandem cell," Appl. Phys. Lett., v ol. 92, p. 3, May 2008.  A. M. Zachary, I. L. Bolotin, D. J. Asunskis, A. T. Wroble, and L. Hanley, "Cluster Beam Deposition of Lead Sulfide Nanocrystals into Organic Matrices," Acs Applied Materials & Interfaces, vol. 1, pp. 1770 1777, 2009.  M. A. Green, "Photovoltaics: technology overview," Energy Policy, vol. 28, pp. 989998, 2000.  A. Goetzberger, C. Hebling, and H. W. Schock, "Photovoltaic materials, history, status and outlook," Materials Science & Engineering R Reports, vol. 40, pp 1 46, 2003.  V. I. Klimov, "Detailed balance power conversion limits of nanocrystal quantum dot solar cells in the presence of carrier multiplication," Appl. Phys. Lett., vol. 89, p. 3, Sep 2006.  M. C. Hanna and A. J. Nozik, "Solar conversion efficiency of photovoltaic and photoelectrolysis cells with carrier multiplication absorbers," J. Appl. Phys., vol. 100, p. 8, Oct 2006.  M. T. Trinh, A. J. Houtepen, J. M. Schins, T. Hanrath, J. Piris, W. Knulst, A. Goossens, and L. D. A. Siebbeles, "In spite of recent doubts carrier multiplication does occur in PbSe nanocrystals," Nano Lett., vol. 8, pp. 17131718, Jun 2008.  G. I. Koleilat, L. Levina, H. Shukla, S. H. Myrskog, S. Hinds, A. G. Pattantyus Abraham, and E. H. Sargent, "Efficien t, stable infrared photovoltaics based on solutioncast colloidal quantum dots," ACS Nano, vol. 2, pp. 833840, May 2008.
185  J. M. Luther, M. Law, M. C. Beard, Q. Song, M. O. Reese, R. J. Ellingson, and A. J. Nozik, "Schottky Solar Cells Based on Colloi dal Nanocrystal Films," Nano Lett., vol. 8, pp. 34883492, Oct 2008.  A. D. Becke, "A New Mixing Of HartreeFock And Local Density Functional Theories," J. Chem. Phys., vol. 98, pp. 1372 1377, 1993.  C. T. Lee, W. T. Yang, and R. G. Parr, "Development Of The ColleSalvetti CorrelationEnergy Formula Into A Functional Of The ElectronDensity," Phys. Rev. B, vol. 37, pp. 785789, 1988.  S. H. Vosko, L. Wilk, and M. Nusair, "Accurate SpinDependent Electron Liquid Correlation Energies For Loc al Spin Density Calculations A Critical Analysis," Can. J. Phys., vol. 58, pp. 12001211, 1980.  P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, "Ab Initio Calculation Of Vibrational Absorption And Circular Dichroism Spectra Usin g Density Functional ForceFields," J. Phys. Chem., vol. 98, pp. 11623 11627, 1994.  A. F. Voter, "Parallel replica method for dynamics of infrequent events," Phys. Rev. B, vol. 57, pp. 1398513988, 1998.  M. R. Sorensen and A. F. Voter, "Temperatureaccelerated dynamics for simulation of infrequent events," J. Chem. Phys., vol. 112, pp. 95999606, Jun 2000.  A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, "New Algorithm For Monte Carlo Simulation Of Ising Spin Systems," Journal of Computati onal Physics, vol. 17, pp. 1018, 1975.  D. T. Gillespie, "General Method For Numerically Simulating Stochastic Time Evolution Of Coupled Chemical Reactions," Journal of Computational Physics, vol. 22, pp. 403434, 1976.  A. F. Voter, F. Montal enti, and T. C. Germann, "Extending the time scale in atomistic simulation of materials," Ann. Rev. Mater. Res., vol. 32, pp. 321346, 2002.  G. Henkelman and H. Jonsson, "A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives," J. Chem. Phys., vol. 111, pp. 70107022, 1999.
186  D. Sheppard, R. Terrell, and G. Henkelman, "Optimization methods for finding minimum energy paths," J. Chem. Phys., vol. 128, 2008.  J. Jacobsen, K. W. Jacobsen, P. Stoltze, and J. K. Norskov, "Island ShapeInduced Transition From 2d To 3d Growth For Pt/Pt(111)," Phys. Rev. Lett., vol. 74, pp. 22952298, Mar 1995.  L. J. Xu and G. Henkelman, "Adaptive kinetic Monte Carlo for first principles accelerated dynamics," J. Chem. Phys., vol. 129, 2008.  G. Henkelman and H. Jonsson, "Long time scale kinetic Monte Carlo simulations without lattice approximation and predefined event table," J. Chem. Phys., vol. 115, pp. 96579666, 2001.  A. F. Voter, "Hyper dynamics: Accelerated molecular dynamics of infrequent events," Phys. Rev. Lett., vol. 78, pp. 3908 3911, 1997.  G. Henkelman and H. Jonsson, "Multiple time scale simulations of metal crystal growth reveal the importance of multiatom surface processes," Phys. Rev. Lett., vol. 90, 2003.  P. Puri and V. Yang, "Effect of particle size on melting of aluminum at nano scales," Journal of Physical Chemistry C, vol. 111, pp. 1177611783, 2007.
187 BIOGRAPHICAL SKETCH Dong h wa Lee was born in South Korea. He earned his bachelors degree in materials science and engineering in 2004 from Korea Universi ty, South Korea. He joined the Department of Materials Science and E ngineering at the University of Florida in 2005. He joined the doctoral degree at the Uni versity of Florida (UF), Gainesville. He started working with Prof. Hans J. Seifert for a year. After Prof. Seifert left, Lee joined the computational science focus group under Prof. Simon R. Phillpot. He received his M.S in materials science and engineeri ng from the University of Florida in fall of 2007. He is complet ing his Ph. D. degree in August 2010. His research focused on using computational simulation methods to characterize different material system.